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?
 A: There are lots of problems with these kind of energy arguments. Here is one: 
Let  ${\bf B}_{\rm ext}$ be a uniform magnetic field in which is immersed  a compactly supported  current distribution ${\bf J}$, which is itself the source of a  field ${\bf B}_{\rm int}$.
We can write ${\bf B}_{\rm ext}= \nabla\times {\bf A}_{\rm ext}$ where ${\bf A}_{\rm ext} = ({\bf B}_{\rm ext}\times {\bf r})/2$ and set out to compute the cross term
$$
 \int  \frac 1 {\mu_0} {\bf B}_{\rm int}
\cdot {\bf  B}_{\rm ext}\, d^3x= \int \frac 1{\mu_0}  {\bf B}_{\rm int}
\cdot (\nabla\times {\bf A}_{\rm ext})\, d^3 x
$$
in the total field energy.
By using the identity
$$
\nabla \cdot ({\bf B}_{\rm int} \times {\bf A}_{\rm ext})= (\nabla\times {\bf  B}_{\rm int})\cdot {\bf A}_{\rm ext} - {\bf B}_{\rm int}\cdot(\nabla \times {\bf A}_{\rm ext})
$$ 
and discarding the total divergence, we  get 
$$
 \left\{\frac 1 {\mu_0} \int  {\bf B}_{\rm int}
\, d^3x\right\}\cdot {\bf  B}_{\rm ext}  =  \frac 1 {\mu_0}  \int {\bf B}_{\rm int}\cdot {\bf  B}_{\rm ext} 
\, d^3x\\
=
\frac 1 {\mu_0}  \int  {\bf B}_{\rm int}\cdot (\nabla\times {\bf A}_{\rm ext})\, d^3x\\
\stackrel{?}{=}\int \frac 1 {\mu_0} (\nabla\times {\bf B}_{\rm int})\cdot {\bf A}_{\rm ext}\, d^3x\\
=\int {\bf J}\cdot {\bf A}_{\rm ext}\, d^3x\nonumber\\
= \int {\bf J} \cdot ({\bf B}_{\rm ext} \times {\bf r})/2\, d^3x\nonumber\\
= \left\{   \int \frac 12 ({\bf r}\times {\bf J})\, d^3x \right\} \cdot {\bf B}_{\rm ext} \nonumber\\
={\bf \mu}\cdot {\bf B}_{\rm ext},{\phantom \int} \nonumber
$$
In the last line 
$$
\mu\equiv   \int \frac 12 ({\bf r}\times {\bf J})\, d^3x
$$
is the standard definition of te magnetic moment of a currebnt distribution.
From this manipulation we might try to deduce that he magnetic moment is given by 
$$
{\mu}\stackrel{?}{=}  \frac 1 {\mu_0} \int_{{\mathbb R}^3}    {\bf B}_{\rm int}\, d^3x=  \frac 1 {\mu_0} \int_{{\mathbb R}^3}  ( {\bf B}-  {\bf B}_{\rm ext})\, d^3x,
$$
where $ {\bf B}$ is the total magnetic field.
The  reason for the "?" over some  of the equals signs  above  is that $|B_{\rm int}|= O(R^{-3})$ and $|A_{\rm ext}|=O(R)$ near infinity, so  the integrated out term 
$$
\int  ({\bf B}_{\rm int} \times {\bf A}_{\rm ext})\cdot {\bf n} \,d|S|
$$ is $O(1)$ and may not be neglected.  Indeed the integral
$$
\int_{{\mathbb R}^3}   {\bf B}_{\rm int}\, d^3x
$$
is only conditionally convergent. It  depends on  the aspect ratio of the region in which it is evaluated.   The formula is correct for an the case when ${\bf B}_{\rm ext}$ is generated by an infinitely long solenoid of finite radius. 
However, even when  correctly evaluated as ${\mu}\cdot {\bf B}_{\rm ext}$, this cross term  is equal to minus the dipole/field interaction energy, which  is correctly given by $-{\mu}\cdot {\bf B}_{\rm ext}$. The sign changes  because we have  ignored the work done by various  induced EMF's when we move the currents.
A  good reference for this issue is ``Magnetic dipole orientation energetics'', G. H. Goedecke, Roy C. Wood, and Paul Nachman. Am. J. Phys. 67 45 (1999)
A: You get symmetric equations for both fields if your magnetic dipole is built from north and south poles, instead of a current loop. Check Wikipedia's Internal magnetic field of a dipole.
Intuitively, with magnetic poles, the field inside the dipole is parallel to the (nearby) outside field, instead of anti-parallel; and that's the same situation as for the electric dipole formed by a positive and a negative charges.
As for the paradoxical case $-$ the spheres of current loops $-$ I believe that you have to take into account the energy of each electrical current $\mathrm{d}I$ under the magnetic field generated by the remaining of the sphere, $\mathrm{d}U = \mathrm{d}I\int B\, \mathrm{d}A$, in line with, e.g., Feynman's Lectures.
A: This is not an answer, but to long for a comment. The only thing I want to show is why the inside directions are different for electric and magnetic charges.
An electric dipole is made of separated electric charges. Having two conducting bodies and moving some amount of electrons from the one to the other body we get something like this (from Wikipedia):

Since the electric field of a isolated subatomic particles is equally distributed in space around the charge the sketched result is highly intuitive for conglomerates of charges.

A macroscopic magnetic dipole is made of aligned elementary magnetic dipoles, mostly this dipoles are electrons. This is the principle behind  permanent magnets.
What happens with the elementary magnetic dipoles you can repeat on the macroscopic level. Imagine you have a lot of bar magnets. You stick them together and what you get is a good model of what happens inside any permanent magnet.
