I'm not sure, what you are really considering. I guess it's magnetostatics. Then the multipole expansion for the magnetic field, given a current density ##\vec{j}(\vec{x})## which is ##0## outside and at the boundary of a sphere of radius ##a## can be derived as follows:
It all starts, of course, with Maxwell's equations in differential form, and I'm assuming the fundamental vacuum form, i.e., ##\vec{j}## is the total current density. Then you have (in SI units)
$$\vec{\nabla} \times \vec{B}=\mu_0 \vec{j}, \quad \vec{\nabla} \cdot \vec{B}=0.$$
From the 2nd equation you know that there's a vector potential ##\vec{A}## such that
$$\vec{B}=\vec{\nabla} \times \vec{A}.$$
Plugging this into the 1st Maxwell equation (Ampere's Law) you get
$$\vec{\nabla} \times (\vec{\nabla} \times \vec{A})=\mu \vec{j}.$$
That can be rewritten as
$$\vec{\nabla} (\vec{\nabla} \cdot \vec{A})-\Delta \vec{A}=\mu_0 \vec{j}. \qquad (*)$$
Now the vector potential is determined from the physics only up to a "gauge transformation", i.e., you can add any gradient field to ##\vec{A}## without changing the physics at all. This means you can impose some convenient constraint on ##\vec{A}##. In magnetostatics the most convenient "choice of gauge" is the Coulomb gauge, which demands
$$\vec{\nabla} \cdot \vec{A}=0.$$
Then (*) simplifies to
$$-\Delta \vec{A}=\mu_0 \vec{j}.$$
For the Maxwell equations of magnetostatics to be solvable at all you must have
$$\vec{\nabla} \cdot \vec{j}=\frac{1}{\mu_0} \vec{\nabla} \cdot (\vec{\nabla} \times \vec{B})=0,$$
which just describes the conservation of charge for a static configuration.
The solution for the potential is now given by the Green's function of ##-\Delta##, which is well-known from electrostatics:
$$\vec{A}(\vec{x})=\mu_0 \int_{B_a} \mathrm{d}^3 x' \frac{\vec{j}(\vec{x}')}{4 \pi |\vec{x}-\vec{x}'|}.$$
We also know from electrostatics that the Green's function can be expanded in terms of spherical harmonics. Using the notation ##r_<=\text{min}(r,r')## and ##r_>=\text{max}(r,r')## it reads
$$G(\vec{x},\vec{x}')=\frac{1}{4 \pi |\vec{x}-\vec{x}'|} = \sum_{\ell=0}^{\infty} \frac{1}{2 \ell +1} \frac{r_{<}^{\ell}}{r_{>}^{\ell+1}} \sum_{m=-\ell}^{\ell} \text{Y}_{\ell m}(\vartheta,\varphi) \text{Y}_{\ell m}^*(\vartheta',\varphi').$$
Now for ##r=|\vec{x}|>a## (i.e., for a point outside of the current distribution) you thus get
$$\vec{A}(\vec{x}) = \sum_{\ell=0}^{\infty} \frac{\mu_0}{(2 \ell +1) r^{\ell+1}} \sum_{m=-\ell}^{\ell} \text{Y}_{\ell m}(\vartheta,\varphi) \int_{B_a} \mathrm{d^3 x'} \mathrm{d} \varphi r^{\prime \ell} \text{Y}_{\ell m}^*(\vartheta',\varphi') \vec{j}(\vec{x}').$$
Defining the "vector multipole moments" (I'm not sure, where to find this in the literature in this way, so maybe my convention concerning coefficients is differen from the literature) such that
$$\vec{A}(\vec{x}) = \frac{\mu_0}{4 \pi} \sum_{\ell=0}^{\infty} \sum_{m=-\ell}^{\ell} \frac{\vec{M}_{\ell m}}{r^{\ell+1}},$$
you have
$$\vec{M}_{\ell m} = \frac{4 \pi}{2 \ell+1} \int_{B_a} \mathrm{d}^3 x' r^{\prime \ell} \text{Y}_{\ell m}^*(\vartheta',\varphi') \vec{j}(\vec{x}').$$
Now let's look at the "monopole term", i.e., ##\ell=m=0##. Since ##Y_{00}=1/\sqrt{4 \pi}## you simply need the integral
$$\vec{M}_{00} = \sqrt{4 \pi} \int_{B_a} \mathrm{d}^3 x' r^{\prime \ell} \vec{j}(\vec{x}').$$
To show that this vanishes, given that ##\vec{\nabla} \cdot \vec{j}=0##. Just consider
$$\partial_k (x_j j_k)=\delta_{jk} j_k + x_j \partial_k j_k=j_j$$
and use Gauss's integral theorem
$$\int_{B_a} \mathrm{d}^3 x \partial_k (x_j j_k) = \int_{B_a} \mathrm{d}^3 x j_j = \int_{\partial B_a} \mathrm{d}^2 f_k (x_j j_k)=0,$$
because by assumption ##\vec{j}(\vec{x})=0## for ##r=|\vec{x}| \geq a##. Thus you indeed have ##M_{00}=0##, and the expansion starts with the dipole contribution.
