electromagnetism - Electromagnetic field energy "paradox"

The starting point is the energy density in electromagnetic fields:

$$ u = \frac{1}{2} \left( \epsilon_0 |\mathbf{E}|^2 + \frac{1}{\mu_0} |\mathbf{B}|^2 \right) \tag{1}$$

The "paradox" is if we use this to test whether it is lower energy for electric dipoles to align (or anti-align) with each other's field it gives the correct answer, but for magnetic dipoles it gives the opposite answer.

---

I ran across this "paradox" several years ago, and calculating everything out is messy and distracting, so if you will forgive me I will do my best to lay this out in such a way that none of the integrals need to be worked out explicitly. Obviously there is no real paradox, the question is what is wrong with the reasoning or what physics is being overlooked that leads to a contradiction.

If we consider two dipoles, they should have the lowest energy when the dipoles are pointing in opposite directions, so that each dipole is aligned with the field of the other:

The goal will be to compare the energy in the fields for the aligned and anti-aligned case of uniformly polarized (or magnetized) spheres, to see how the energy depends on orientation of the dipoles. First we'll start with electric dipoles, but the concept is the same for considering magnetic dipoles.

$$ U = \frac{\epsilon_0}{2} \int |\mathbf{E}_{\text{dipole1}} + \mathbf{E}_{\text{dipole2}}|^2 d^3r = \frac{\epsilon_0}{2} \int \left( |\mathbf{E}_1|^2 + 2 \mathbf{E}_1 \cdot \mathbf{E}_2 + |\mathbf{E}_2|^2 \right) d^3r \tag{2}$$

Only the $\mathbf{E}_1 \cdot \mathbf{E}_2$ term will depend on the orientation of the dipoles.

Since we need to integrate over all space, we cannot neglect the field inside the spheres. This actually turns out to be important, and is where the electric and magnetic dipoles differ. The electric field lines switch directions in the polarized sphere, while the magnetic field lines do not.

$$\mathbf{E}_\text{dipole} = \begin{cases} \frac{1}{4\pi\epsilon_0 r^3} ( 3(\mathbf{p}\cdot\mathbf{\hat{r}})\mathbf{\hat{r}} - \mathbf{p}) &, r > R \\ -\frac{1}{3\epsilon_0}\frac{\mathbf{p}}{\frac{4}{3}\pi R^3} = -\frac{1}{3\epsilon_0} \mathbf{P} &, r < R \end{cases} \tag{3}$$

$$\mathbf{B}_\text{dipole} = \begin{cases} \frac{\mu_0}{4\pi r^3} ( 3(\mathbf{m}\cdot\mathbf{\hat{r}})\mathbf{\hat{r}} - \mathbf{m}) &, r > R \\ +\frac{2\mu_0}{3}\frac{\mathbf{m}}{\frac{4}{3}\pi R^3} = +\frac{2\mu_0}{3}\mathbf{M} &, r < R \end{cases} \tag{4}$$

So when evaluating the field energy, $\mathbf{E}_1 \cdot \mathbf{E}_2$ breaks into three pieces. One is the region outside both spheres, and the other two are the regions in either sphere (which due to symmetry must give the same value). When the dipoles are aligned with each other's fields, we can see that all regions contribute negatively. And when the dipoles are anti-aligned with each other's fields, all regions will contribute positively.

So, as expected, the lower energy orientation is when each dipole is aligned with the other's field.

If this is not clear, it is easiest to analyze in the limit of a perfect dipole (or approximating the case where $R \ll d$ where $R$ is the radius of the sphere and $d$ is the distance between them). Then the inside one sphere becomes just $$\int \left(\frac{-1}{3}\mathbf{p}_1 \ \delta^3(\mathbf{r_1}) \right) \cdot \mathbf{E}_2(\mathbf{r}) d^3r \tag{5}$$ The outside of the spheres case can be analyzed noting that we can choose the origin at one dipole, and expand the field of the other in multipole moments. Since the other dipole is offset, the dipole moment will no longer be the only non-zero moment in this expansion, but since we are eventually integrating $\mathbf{E}_1 \cdot \mathbf{E}_2$ and the spherical harmonic components are orthogonal, only the dipole-dipole term will contribute. In the far field, the dipole moment in this offset expansion will clearly still be in the same direction as the old moment, just a change in magnitude.

If we actually worked out the integral, we should find the energy of orientation is:

$$U = - \mathbf{p}_1 \cdot \mathbf{E}_2(\mathbf{r}_1) = - \mathbf{p}_2 \cdot \mathbf{E}_1(\mathbf{r}_2) \tag{6}$$

Now, even without working out the integral to verify that, we can see something is wrong. If we look at magnetic dipole interaction, the fields are equivalent to the electric dipoles situation when outside the spheres, but are in the opposite direction inside the spheres. This means the same calculation for the magnetic dipoles cannot work. When the magnetic dipoles are aligned with each other's field, the region outside the spheres will still contribute negatively (all is good there), while the region inside the spheres will now contribute positively. Something is wrong.

Indeed, if one does the calculation for the magnetic dipoles, the result will even have the correct magnitude ... but the wrong sign! This is the paradox.

The calculation works fine for electric dipoles, but claims the energy is lowest when the magnetic dipoles are parallel with each other (anti-aligned to each other's fields).

What is going on?