Magnetic Multipole Tensor When the electric scalar potential is expanded into spherical coordinates, one gets
\begin{align}
\phi (\vec r) =  \frac{1}{4\pi\varepsilon_0} \sum_{l=0}^{\infty} \sum_{m=-l}^l \sqrt{\frac{4\pi}{2l+1}} Y_{lm} (\theta, \varphi) \frac{Q_{lm}}{r^{l+1}} 
\end{align}
with the electric multipole tensor
\begin{align}
Q_{lm} = \sqrt{\frac{4\pi}{2l+1}}\int r'^{l} Y^*_{lm} (\theta', \varphi') \rho(\vec r\,^\prime) d^3 r' \,.
\end{align}
This is simple and clear. But when the same is done for the magnetic vector potential, one gets
\begin{align}
\vec A(\vec r) = \mu_0 \sum_{l=0}^{\infty} \frac{1}{2l+1} \frac{1}{r^{l+1}}\sum_{m=-l}^l Y_{lm} (\theta, \varphi) \int  r'^l  Y^*_{lm} (\theta', \varphi') \vec j (\vec r\,^\prime)\, d^3 r' \, ,
\end{align}
and I don't see a way to define a magnetic multipole tensor like in the electric case, because now, in the integral, there is a vector, what makes the definition of a tensor of arbitrary rank quite awkward. So how can a magnetic multipole tensor be defined in this case?
 A: One simple solution is to use the scalar magnetic potential. This method works for magnetostatical problems with permanent magnets (that is, without free currents).
In macroscopic electrodynamics and in the absence of free currents we have
$$\nabla \times \vec H = 0.$$
Therefore, given a simply connected domain, there exists a scalar potential $\psi$, such that:
$$ \vec H = - \nabla \psi. $$
(Formally, this is completely analogous to the existence of the scalar potential for the electrostatic field.)
From the relation $\vec B = \mu_0 (\vec H + \vec M)$ we get (note that $\nabla \cdot \vec B = 0$ does always hold, even in macroscopic electrodynamics):
$$ \nabla \cdot \vec H = -\nabla \cdot \vec M. $$
Inserting the above relations one gets:
$$ \nabla \cdot \vec H = \nabla \cdot (-\nabla \psi) = -\nabla \cdot \vec M. $$
Rearranging gives:
$$ \Delta \psi = \nabla \cdot \vec M =: -\rho_m. $$
As this is formally equal to the Poisson equation from electrodynamics one can easily do a multipole expansion of $\psi$ in terms of multipoles of $\rho_m$ (the quantity $\rho_m$ is called bound magnetic charge). Outside of the permanent magnets the magnetization is zero, and therefore $\vec B$ can easily be computed from $\psi$ as:
$$ \vec B = -\mu_0 \nabla \psi. $$
Note, that all this does not imply there are magnetic monopoles! $\nabla \cdot \vec B = 0$ does hold all the time! $\rho_m$ is not a magnetic charge in the sense, that it is a source of $\vec B$, it is an artefact of handling some currents as bound and separating the flux density in magnetization and field strength.
A: There's no reason why you can't use the same method you used for the electric multipole expansion for magnetic multipole expansion.  
The indices $l,m$ indicate the order of multipole expansion used in the approximation in those equations.  
Since the $\vec A(\vec r)$ field you provided implies $\nabla \cdot \vec B(\vec r)=0$, i.e. there is no magnetic monopole, I would exploit the symmetry of magnetic dipole and use a power series of closed circuit loops
$$\vec A(\vec r)=\sum_{n=0}^{\infty}\vec A_n(\vec r)$$
where 
$$\vec A_n(\vec r)=\frac {\mu_0}{4\pi} \frac{1}{r^{n+1}} \oint_{v^{'}} {(r^{'})}^n P_n(cos\theta^{'})\vec J(\vec r^{'})d\tau^{'} $$
and $$cos\theta{'}= \hat{r} \cdot \hat r^{'}$$ $$r >> r^{'}$$
Since there is no known magnetic monopole in general - at least everywhere except possibly Stanford Univerity, then $\vec A_0(\vec r)=0$. 
And $\vec A_1(\vec r)$ is the dipole used in magnetostatics, $\vec A_2(\vec r)$ is the quadrupole, etc.
