I've never worked with magnetic quadrupoles. Jackson doesn't appear to give any explicit expressions for magnetic quadrupole moments. I did find the following expression in Morse and Feshbach's
Methods of Mathematical Physics for a general distribution of current density ##\vec J##:
##\mathbf{M} = \frac{1}{2}\int\left[ \vec r (\vec r \times \vec J) + (\vec r \times \vec J) \vec r \right] d^3r##
This is dyadic notation. Thus to get ##{M}_{xy}##, you dot the integrand from the left with a unit vector in the x-direction and dot the integrand from the right with a unit vector in the y-direction. So,
##M_{xy} = \frac{1}{2} \int \left[x(zJ_x-xJ_z) + (yJ_z-zJ_y)y\right]d^3r##
## \,\,\,\,\,\,\,\,\,\,\,\,= \frac{1}{2} \int \left[ (y^2-x^2)J_z + xzJ_x - yzJ_y \right] d^3x##
##M_{zz} = \int \left( zxJ_y - zyJ_x \right) d^3r##
etc.
The dimensions of ##M## are current times distance cubed; whereas, the dimensions of your expression for ##M_{ij}## are current times distance to the fourth.
I applied this dyadic definition to the case of a magnetic quadrupole formed from two current loops with opposite currents.
View attachment 257972
I got the same result for the magnetic quadrupole moment components as given in
this link. See pages 9 through 11. The final result is at the bottom of page 11.
The vector potential ##\mathbf A## produced by this quadrupole is given by the last term shown in the equation at the top of page 11. But this is only for the far field. I would imagine that the near field would be very tedious to work out.