quadrupole moment and higher for simple current loop I'm working on a field configuration that needs to die off rapidly and I got to a $1/r^5$ dependence with canceling the dipole moment of the system cleverly, but to go get better arrangements I need to figure out the higher order moments of a simple loop as well. 
Griffiths (page 243, 3rd ed) refers to this integral $\oint (r')^{n}P_n(cos \theta ') d L'$, with $P_n$ as the Legendre polynomials but I am unclear on how to go about evaluating it for higher moments because for the dipole term Griffiths does a clever trick that I'm not super sure extrapolates. I feel like the current loop is such a simple, tidy system it's higher moments would be tabulated somewhere just as an example of how higher orders are calculated, but I'm unable to find them after a few a days of googling and librarying.
 A: Whenever $|\mathbf{r}|$ lies far outside the loop, we can expand the vector potential $\mathbf{A}(\mathbf{r})$ in an irregular solid harmonic series
$$\mathbf{A}(\mathbf{r})=\frac{\mu_0}{4\pi}\int\frac{\mathbf{J}(\mathbf{r}')}{|\mathbf{r}-\mathbf{r}
'|}dV=\frac{\mu_0}{4\pi}\sum_{L=0}^\infty\sum_{m=-L}^L I_{Lm}(\mathbf{r})\langle R_{Lm},\mathbf{J}\rangle$$
where $R_{Lm},I_{Lm}$ are the regular and irregular real solid harmonic functions and $\langle R_{Lm},\mathbf{J}\rangle$ is the integral inner product between the two (aka, the $Lm$ multipole moment). 
Letting the loop with radius $\rho$ lie in the $xy$-plane with $z$-axis height $h$, the interior moments $\langle R_{Lm},\mathbf{J}\rangle$ can be computed as
$$\langle R_{Lm},\mathbf{J}\rangle=\int_0^{2\pi}R_{Lm}\left(\sqrt{\rho^2+h^2},\arctan(h,\rho),\phi\right)\mathbf{J}(\phi)\,\rho d\phi.$$
Letting the loop current flow counterclockwise, in Cartesian coordinates we have 
$$\mathbf{J}(\phi)=
\left(
\begin{array}{cc}
-J_0\sin(\phi)\\
J_0\cos(\phi)\\0
\end{array}
\right)
$$
and we obtain the following moments up to octupole order (copy paste into Mathematica and hit Cell > Convert To > StandardForm to properly format it):
SphericalHarmonicYr[L_, m_, \[Theta]_, \[Phi]_] := (-1)^m*
   Sqrt[((2*L + 1)*(L - Abs[m])!)/((2*Pi)*(L + Abs[m])!)]*
   LegendreP[L, Abs[m], Cos[\[Theta]]]*
       Piecewise[{{Sin[Abs[m]*\[Phi]], m < 0}, {Sqrt[1/2], 
      m == 0}, {Cos[Abs[m]*\[Phi]], m > 0}}]; 
SolidHarmonicIr[l_, m_, 
   r_, \[Theta]_, \[Phi]_] := (Sqrt[(4*Pi)/(2*l + 1)]*
     SphericalHarmonicYr[l, m, \[Theta], \[Phi]])/r^(l + 1); 
SolidHarmonicRr[l_, m_, r_, \[Theta]_, \[Phi]_] := 
  Sqrt[(4*Pi)/(2*l + 1)]*r^l*
   SphericalHarmonicYr[l, m, \[Theta], \[Phi]]; 
J[\[Phi]_] := J0*{-Sin[\[Phi]], Cos[\[Phi]], 0}; 
Grid[Simplify[
  Table[MatrixForm[
    Integrate[\[Rho]*J[\[Phi]]*
      SolidHarmonicRr[L, m, Sqrt[\[Rho]^2 + h^2], 
       ArcTan[h, \[Rho]], \[Phi]], {\[Phi], 0, 2*Pi}]], 
       {L, 0, 3}, {m, -L, L}], 
  Assumptions -> {\[Rho] > 0, Element[h, Reals]}], Frame -> All]




The code above computes the moments out to $L=3$, but you can change 3 to whatever multipole order is necessary. Here is a numerical verification that the irregular solid harmonic series $\frac{\mu_0}{4\pi}\sum_{L=0}^\infty\sum_{m=-L}^L I_{Lm}(\mathbf{r})\langle R_{Lm},\mathbf{J}\rangle$ and the integral $\frac{\mu_0}{4\pi}\int\frac{\mathbf{J}(\mathbf{r}')}{|\mathbf{r}-\mathbf{r}
'|}dV$ are actually the same thing (again, hit Cell > Convert To > StandardForm to properly format the appearance):
Subscript[Q, L_, m_] := 
  Integrate[\[Rho]*J[\[Phi]]*
    SolidHarmonicRr[L, m, Sqrt[\[Rho]^2 + h^2], 
     ArcTan[h, \[Rho]], \[Phi]], {\[Phi], 0, 2*Pi}]; 
r = N[{5, 6, 7}]; 
rSph = CoordinateTransform["Cartesian" -> "Spherical", r]; 
h = 2; 
\[Rho] = 1; 
J0 = 2; 
NIntegrate[(J[\[Phi]]/
    Norm[r - {\[Rho]*Cos[\[Phi]], \[Rho]*Sin[\[Phi]], 
       h}])*\[Rho], {\[Phi], 0, 2*Pi}]
Sum[(Subscript[Q, L, m]*SolidHarmonicIr[L, m, ##1] & ) @@ rSph, {L, 
  0, 5}, {m, -L, L}]
Clear[h, \[Rho], J0]


{-0.0471743, 0.039312, 0.}
{-0.0471859, 0.0393216, 0.}

Alternately, as an additional simplification, note that $\mathbf{J}$ can be expressed in spherical vector coordinates as
$$\mathbf{J}(\phi)=\left(
\begin{array}{cc}
0\\
0\\
J_0
\end{array}
\right)$$
giving (this time up to $L=4$):
J[\[Phi]_] := {0, 0, J0}; 
Grid[Simplify[
  Table[MatrixForm[
    Integrate[\[Rho]*J[\[Phi]]*
      SolidHarmonicRr[L, m, Sqrt[\[Rho]^2 + h^2], 
       ArcTan[h, \[Rho]], \[Phi]], {\[Phi], 0, 2*Pi}]], 
       {L, 0, 4}, {m, -L, L}], 
  Assumptions -> {\[Rho] > 0, Element[h, Reals]}], Frame -> All]




although for this coordinate choice I'm a bit confused as to why the $L=0$ coefficient is not zero (will have to check that later).
The magnetic field can then be obtained from the above coefficients $\langle R_{Lm},\mathbf{J}\rangle$ by
$$\mathbf{B}(\mathbf{r})=\nabla\times\mathbf{A}(\mathbf{r})=\frac{\mu_0}{4\pi}\sum_{L=0}^\infty\sum_{m=-L}^L \nabla\times I_{Lm}(\mathbf{r})\langle R_{Lm},\mathbf{J}\rangle.$$
