ClearAll["Global`*"]

(* Smooth load used in [de Goes and James 2017] *)
g[r_] := 15 ϵ^4/(8 Pi) Power[r^2+ϵ^2, -7/2];
Sphere2Euc[r_, θ_, ϕ_] := {{r Sin[θ]Cos[ϕ]}, {r Sin[θ]Sin[ϕ]}, {r Cos[θ]}};

(* Voigt index mapping *)
VoigtIndex[i_, j_]:= Which[i==j,i,(i==2 && j==3 )|| (i==3 && j==2),4,(i==1&&j==3) || (i==3 &&j==1),5,(i==1&&j==2 ||i==2&&j==1),6];
(* Convert reduced 6x6 matrix to a full fourth-order tensor *)
MatTensor[M_] := (
rtnC=ConstantArray[0, {3,3,3,3}];
For[i=1, i<=3, i++, 
For[j=1,j<=3,j++, 
For[k=1,k<=3,k++, 
For[l=1,l<=3,l++, 
rtnC[[i,j,k,l]]=M[[VoigtIndex[i,j],VoigtIndex[k,l]]]
]
]
]
];
Return[rtnC];
);
ElasSymbol[C_, ξ_] := (
rtn=ConstantArray[0,{3,3}]; 
For[i=1,i<=3,i++,
For[j=1,j<=3,j++,
For[k=1,k<=3,k++,
For[l=1,l<=3,l++,
rtn[[i,k]]+=C[[i,j,k,l]]ξ[[j,1]]ξ[[l,1]]
]
]
]
];
Return[rtn];
);

(* Transform smooth load by Hankel. [isReg] True for g(r); False for δ(r) *)
hatg[ζ_,isReg_] := If[isReg, Power[2 Pi,3/2]/Sqrt[ζ]Integrate[Simplify[BesselJ[1/2,r ζ]  g[r]Power[r,3/2]], {r, 0, Infinity}, Assumptions->ζ≥0&& ϵ>0],1] ;
Print["Fourier g(ζ)=",hatg[ζ,True]];

(* Radial term R_l *)
R[r_, l_]:=Integrate[hatg[ζ, True]*SphericalBesselJ[l, r ζ],{ζ,0,Infinity}, Assumptions->ϵ>0 && l≥0 && r>0];
For [l=0, l≤2, l++,
Print["l=",l,", R=",R[r, l]];
];

(* Material term P_l^m, integrating on a unit sphere *)
Plm[l_, m_, invS_] :=  Integrate[invS*Conjugate[SphericalHarmonicY[l,m,  Θ, Φ]]Sin[Θ], {Θ, 0, Pi}, {Φ, 0, 2*Pi}];

(* Assemble SH expansion *)
GeneralGF[r_, θ_, ϕ_, invS_] :=(
result=ConstantArray[0, {3, 3}];
For[l=0, l≤2, l++, 
For[m=-l, m ≤ l, m++, 
(*Print[l," ", m, " ", MatrixForm[Plm[l, m, invS]]];*)
result+= 1/(2 Pi^2)I^l SphericalHarmonicY[l,m, θ,ϕ] R[r,l]Plm[l, m, invS];
]
];
Return[Simplify[ExpToTrig[result]]];
);

(* Apply to isotropic material *)
isoM = {
{λ+2μ,λ,λ,0,0,0},
{λ,λ+2μ,λ,0, 0, 0},
{λ,λ,λ+2μ,0, 0, 0}, 
{0, 0, 0, μ, 0, 0}, 
{0, 0, 0, 0, μ, 0},
{0, 0, 0, 0, 0, μ}
};
ξ= Sphere2Euc[1, Θ, Φ]; invS =Simplify[Inverse[ElasSymbol[MatTensor[isoM], ξ]]];
MatrixForm[Simplify[ GeneralGF[r, θ, ϕ, invS], Assumptions->ϵ>0]]

(* Formula for singular and regularized Kelvinlet [de Goes and James 2017] *)
a = 1/(4Pi μ);
v=λ/(2(μ+λ));
b = Simplify[a/(4(1-v))];
u=Sphere2Euc[r, θ, ϕ];
kel=MatrixForm[Simplify[(a-b)IdentityMatrix[3]/r + b/r^3 u .Transpose[u]]]
regkel = MatrixForm[Simplify[(a-b)IdentityMatrix[3]/Sqrt[r^2+ϵ^2]+b/Sqrt[r^2+ϵ^2]^3 u . Transpose[u]+a ϵ^2/(2 Sqrt[r^2+ϵ^2]^3) IdentityMatrix[3]]]