Solution (source code)

= Solution

Use source coordinates $(x',y',0)$ and observation point $(0,0,z)$. The axial numerator in the <Biot-Savart law> is
$$
\left[d\boldsymbol\ell'\times(-x',-y',z)\right]_z
=x'\,dy'-y'\,dx'.
$$
For $|z|$ much larger than the loop dimensions, replace the denominator $[z^2+x'^2+y'^2]^{3/2}$ by $|z|^3$ to leading order. Transverse components cancel between opposite sides. The loop integral is twice its oriented area, so
$$
B_z\simeq\frac{\mu_0I}{4\pi|z|^3}
\oint(x'\,dy'-y'\,dx')
=\frac{\mu_0I}{4\pi|z|^3}(2\Delta x\Delta y).
$$
Therefore
$$
\boxed{\mathbf B(0,0,z)\simeq\frac{\mu_0\mathbf m}{2\pi|z|^3}.}
$$
It points along $\mathbf m$ on both halves of the dipole axis, not along the signed radial unit vector.