Compute the sum of all the quadropole potential terms in two ways (following Jackson, chapter 4).

ClearAll[y2m, q2m, m, tp, p, tp, pp, rp, phim, r, Pre, Pim, Ps]
$Assumptions = And@@(#>0&)/@{tp, pp, rp, t, p, r};

y2m[m_, r_, t_, p_] := r^2 SphericalHarmonicY[2, m, t, p]

phim[m_] := ((4 Pi)/5) Conjugate[y2m[m, rp, tp, pp]] y2m[m, r, t, p]

Ps = Sum[ phim[m]  , {m,-2,2}] // FullSimplify ;
Pim = Ps // Im // FullSimplify
Pre = Ps // Re // FullSimplify
Assert[ Pim == 0]
P= Pre // TrigExpand

(1/16) (r^2) (rp^2) (1+3 (Cos[2 tp]+Cos[2 θ] (1+3 Cos[2 tp])+4 Cos[2 ϕ-2 ϕp] ((Sin[θ])^2) ((Sin[tp])^2)+4 Cos[ϕ-ϕp] Sin[2 θ] Sin[2 tp]))

ClearAll[i, j]
spherical[r_, t_, p_] := CoordinateTransformData["Spherical"->"Cartesian", "Mapping", {r,t, p}] ;
(*r{ Sin[t]Cos[p], Sin[t]Sin[p], Cos[t]} ;*)

Qij[i_, j_, r_, t_, p_] := 3spherical[r,t,p][[i]] spherical[r,t,p][[j]] -KroneckerDelta[i,j] r^2 ;

Q =Sum[Qij[i,j,rp, tp, pp] spherical[r,t, p][[i]] spherical[r,t, p][[i]] ,{i,1,3}, {j,1,3}]/2  //FullSimplify

(*P
qq = Q // TrigExpand*)

P

P /. {Cos[pp] -> xp/Sin[tp]/rp, Sin[pp] -> yp/Sin[tp]/rp , Cos[tp] -> zp/rp, Cos[p] -> x/Sin[t]/r, Sin[p] -> y/Sin[t]/r, Cos[t]-> z/r} 

(*P /. Cos[t] -> z/r*)

c = {x, y, z}
cp = {xp, yp, zp}
qij[i_, j_] := 3cp[[i]] cp[[j]] -KroneckerDelta[i,j] (xp^2 + yp^2 + zp^2);

q =Sum[qij[i,j] c[[i]] c[[i]] ,{i,1,3}, {j,1,3}]/2  //FullSimplify

SphericalHarmonicY[2, #, θ, ϕ]&/@ {0, 1,2}