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]
![enter image description here](https://i.sstatic.net/wNgCL.png)
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]
![enter image description here](https://i.sstatic.net/FQHx3.png)
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.$$