As per my own answer to my question about the tensor Laplacian, the formula for the tensor Laplacian is
$$(\boldsymbol{\triangle}\mathbf{T})^{i_1\dots i_r}{}_{j_1\dots j_s}=g^{kl}(\nabla(\nabla\mathbf{T}))^{i_1\dots i_r}{}_{j_1\dots j_s~kl}$$
By second order tensor I assume you mean a $(0,2)$ tensor in which case we get
$$(\boldsymbol{\triangle}\mathbf{T})_{ij}=g^{kl}(\nabla(\nabla\mathbf{T}))_{ij~kl}$$
NB: In the following I am using the so called mathematics convention for spherical coordinates, i.e $x=r\cos\theta\sin\phi$, etc.
There are two relevent formulae you'll need to compute the tensor Laplacian. The first of course is the inverse metric (written in the proper $(2,0)$ contravariant form)
$$\mathbf{g}^{-1} =\begin{bmatrix}
\begin{bmatrix}
1\\
0\\
0
\end{bmatrix}\\
\begin{bmatrix}
0\\
(r^2\sin^2\phi)^{-1}\\
0
\end{bmatrix}\\
\begin{bmatrix}
0\\
0\\
1/r^2
\end{bmatrix}
\end{bmatrix}$$
And the covariant derivative, which is going to be applied twice, first for a $(0,2)$ tensor
$$(\nabla \mathbf{T})_{ij~k}=\partial_kT_{ij}-\Gamma^a_{ik}T_{a j}-\Gamma^{a}_{jk}T_{ia}$$
And then for a $(0,3)$ tensor:
$$(\nabla(\nabla\mathbf{T}))_{ijk~l}=\partial_lT_{ijk}-\Gamma^a_{il}T_{ajk}-\Gamma^{a}_{jl}T_{iak}-\Gamma^{a}_{kl}T_{ija}$$
The computations are pretty heavy, so I'm going to use Mathematica to assist me.
I'll paste in the Mathematica code below.
coords = {r, \[Theta], \[Phi]};
n = Length[coords];
p = {r*Cos[\[Theta]] Sin[\[Phi]], r*Sin[\[Theta]] Sin[\[Phi]],
r*Cos[\[Phi]]};
g = Simplify[
Table[D[p, coords[[i]]].D[p, coords[[j]]], {i, 1, n}, {j, 1, n}]];
inverseg = Simplify[Inverse[g]];
\[CapitalGamma] =
Simplify[Table[
inverseg[[k, k]]*
D[p, coords[[k]]].D[D[p, coords[[i]]], coords[[j]]], {k, 1,
n}, {i, 1, n}, {j, 1, n}]];
T = {{Trr[r, \[Theta], \[Phi]], Tr\[Theta][r, \[Theta], \[Phi]],
Tr\[Phi][r, \[Theta], \[Phi]]}, {T\[Theta]r[r, \[Theta], \[Phi]],
T\[Theta]\[Theta][r, \[Theta], \[Phi]],
T\[Theta]\[Phi][r, \[Theta], \[Phi]]}, {T\[Phi]r[
r, \[Theta], \[Phi]], T\[Phi]\[Theta][r, \[Theta], \[Phi]],
T\[Phi]\[Phi][r, \[Theta], \[Phi]]}};
Clear[CovariantDerivative]
CovariantDerivativeSecondOrder[T_] :=
Table[D[T[[i, j]], coords[[k]]] -
Sum[\[CapitalGamma][[a, i, k]] T[[a, j]], {a, 1, n}] -
Sum[\[CapitalGamma][[a, j, k]] T[[i, a]], {a, 1, n}], {i, 1,
n}, {j, 1, n}, {k, 1, n}]
CovariantDerivativeThirdOrder[T_] :=
Table[D[T[[i, j, k]], coords[[l]]] -
Sum[\[CapitalGamma][[a, i, l]] T[[a, j, k]], {a, 1, n}] -
Sum[\[CapitalGamma][[a, j, l]] T[[i, a, k]], {a, 1, n}] -
Sum[\[CapitalGamma][[a, k, l]] T[[i, j, a]], {a, 1, n}], {i, 1,
n}, {j, 1, n}, {k, 1, n}, {l, 1, n}]
Simplify[Table[
Sum[inverseg[[k, l]]*
CovariantDerivativeThirdOrder[
CovariantDerivativeSecondOrder[T]][[i, j, k, l]], {k, 1, n}, {l,
1, n}], {i, 1, n}, {j, 1, n}]]
Here is a screenshot. I've removed the redundant [r,θ,ϕ] all over the place. Note that in Mathematica, the (a,b,c) superscript means an a-fold partial derivative in the first argument of the function and so on.
This would take me quite a while to typeset, and I am unfortunately very busy at the moment so I really don't have time. I will get to it hopefully in the next couple of days.
