\begin{eqnarray} \nabla \cdot \boldsymbol \tau &=& 2 \mu \nabla \cdot \boldsymbol \varepsilon\\ &=& \mu \nabla \cdot \left( \nabla\mathbf{u} + (\nabla\mathbf{u}) ^\mathrm{T} \right)\\ &=& \mu \, \nabla^2 \mathbf{u} \end{eqnarray}
Given $$ \nabla \cdot \boldsymbol u =0 \ $$ How do you go from step 3 to step 4 and get rid of the transpose? Isn't the divergence operator distributive
The divergence of the transpose of a vector gradient is equal to the gradient of the divergence of the vector. You can prove it quickly using some index notation:
$$\nabla \cdot (\nabla \vec{u})^T = \frac{\partial^2 u_j}{\partial x_i \partial x_j} = \frac{\partial}{\partial x_i}\frac{\partial u_j}{\partial x_j} = \nabla(\nabla \cdot \vec{u})$$
I assumed everything is in orthonormal coordinates so no need to raise or lower indices. Since the flow is incompressible by your assumption, this term cancels out, and you're only left with the divergence of the vector gradient (which gives you the vector Laplacian).
I'm going to guess that the RHS really means $$ \nabla^j ( \nabla_i u_j + \nabla_j u_i ) = \nabla_j \nabla_i u^j + \nabla^2 u_i . $$ Now, $$ \nabla_j \nabla_i u^j = [ \nabla_j , \nabla_i ] u^j + \nabla_i (\nabla_j u^j) = R_{ij} u^j + \nabla_i (\nabla_j u^j) . $$ You're probably working in a Ricci flat background and the velocity vector is divergenceless so $\nabla_j \nabla_i u^j=0$ and therefore $$ \nabla^j ( \nabla_i u_j + \nabla_j u_i ) = \nabla^2 u_i . $$
