When mutual inductance are occurring between two coils, why the two coils have the same mutual inductance? I learned that 
$$\epsilon_1 = -M_{12}\frac{di_2}{dt}$$
$$\epsilon_2 = -M_{21}\frac{di_1}{dt}$$
And the book tells us directly that $M_{12} = M_{21}$ without a reason. Is there a mathematical proof for this?   
 A: An elegant and elementary derivation is given by Crawford:"Mutual inductance $M_{12}=M_{21}$", American Journal of Physics, vol 60 , February 1962, p186.
The idea is the following: the total stored energy rate, "power" is $$\frac{du}{dt} = L_1i_1\frac{di_1}{dt} + M_{12}i_1\frac{di_2}{dt} + L_2i_2\frac{di_2}{dt} + M_{21}i_2\frac{di_1}{dt}$$
This can be written as the differential 
$${du} = L_1d(i_1^2/2) + M_{12}i_1di_2 + L_2d(i_2^2/2) + M_{21}i_2di_1$$
Now integrate from $i_1(t=0)=i_2(t=0)=0$ to $I_1,I_2$ and get $$U=\frac{1}{2}L_1I_1^2 +M_{12}I_1I_2+\frac{1}{2}L_2I_2^2+(M_{21}-M_{12})\int {i_2di_1}$$.
In general, the integral $\int {i_2di_1}$ can be anything and will depend on the history of the currents, and thus $U$ is not necessarily single valued function of the present currents $I_1,I_2$ unless $M_{21}=M_{12}$! Now you can see that for a simple medium we must have equality of the mutual inductances.
A: The answer by hyportnex is an elegant solution that does not even explicitly use the definition of mutual inductance. I also learned this in the textbook but had been wondering about whether one can prove this through brutal force. This is, therefore, a complementary proof that gives a straightforward approach by performing integration using the definition. Sometimes it is worthy to work the non-elegant part out explicitly.
$$M_{12}=M_{2\to 1}=\frac{1}{I_2}\phi_{2\to 1}
=\frac{1}{I_2}\int_{S_2}\mathbf{B}_{2\to 1}\cdot d\mathbf{S}_1$$
$$=\frac{\mu_0}{4\pi I_2}\int_{S_2}\oint_{L_1}\frac{I_2d\mathbf{r}_2\times \mathbf{r_{12}}}{r_{12}^3}\cdot d\mathbf{S}_1
=\frac{\mu_0}{4\pi}\int_{S_2}\oint_{L_1}\frac{d\mathbf{r}_2\times \mathbf{r_{12}}}{r_{12}^3}\cdot d\mathbf{S}_1 ,$$
where $\mathbf{S_1}$ represents the surface of the coil $\mathbf{L_2}$, subscript $1\to 2$ indicates a quantity generated at the location of the second coil by the current of first coil, for instance $\mathbf{B}_{1\to 2}$ means the magnetic induction at position $\mathbf{r}_2$ generated by an infinitesimal piece of coil $d\mathbf{r}_1$ with current $I_1$. Also, $\mathbf{r}_{12}=\mathbf{r}_1-\mathbf{r}_2$, $r_{12}=|\mathbf{r}_{12}|$, $\nabla_1\equiv \nabla_{\mathbf{r}_1}$.
We observe that the resultant expression is entirely geometrical.
To be specific, it only depends on the geometrical layouts of the two coils but not the currents. One can similarly write down the corresponding expression for $M_{12}$, and try to show that these two expressions are indeed the same by some manipulations, which are the focus of this write-up.
To proceed further, we note the following mathematical relations:
$$\frac{\mathbf{r}}{r^3}=-\nabla\frac{1}{r} ,$$
and in particular, one has
$$\frac{\mathbf{r_{12}}}{r_{12}^3}=-\nabla_1\frac{1}{r_{12}}=+\nabla_2\frac{1}{r_{12}}.$$
This is because for any scalar function $f=f(r_{12})$, we have, in general,
$$\nabla_1 f(r_{12})=-\nabla_2f(r_{12})=-\nabla_2f(r_{21}) ,$$
and therefore
$$\nabla_1\nabla_1 f(r_{12})=\nabla_2\nabla_2f(r_{12})=\nabla_2\nabla_2f(r_{21}) .$$
We further introduce the vectorial part. To start with, we have
$$\mathbf{A}\times(\mathbf{B}\times \mathbf{C})=(\mathbf{A}\cdot \mathbf{C})\mathbf{B}-(\mathbf{A}\cdot \mathbf{B})\mathbf{C} .$$
we can further introduce the vectorial operator $\nabla$ into the above formula
$$\nabla\times(\mathbf{B}\times \nabla \phi)=(\nabla\cdot \nabla\phi)\mathbf{B}-(\nabla\cdot \mathbf{B})\nabla\phi= -(\nabla\cdot \mathbf{B})\nabla\phi ,$$
where $\phi$ is an arbitrary scalar field that satisfies
$$\nabla\cdot \nabla\phi =0.$$
Now we have all the pieces to continue the calculation, using the Stokes theorem
$$\frac{\mu_0}{4\pi}\int_{S_2}\oint_{L_1}\frac{d\mathbf{r}_2\times \mathbf{r_{12}}}{r_{12}^3}\cdot d\mathbf{S}_1
=\frac{\mu_0}{4\pi}\int_{S_2}\oint_{L_1}\left(d\mathbf{r}_2\times \frac{\mathbf{r_{12}}}{r_{12}^3}\right)\cdot d\mathbf{S}_1$$
$$=\frac{\mu_0}{4\pi}\oint_{L_2}\oint_{L_1}\left[\nabla_1\times\left(d\mathbf{r}_2\times \frac{\mathbf{r_{12}}}{r_{12}^3}\right)\right]\cdot d\mathbf{r}_1
=\frac{\mu_0}{4\pi}\oint_{L_2}\oint_{L_1}\left[\nabla_1\times\left(d\mathbf{r}_2\times (-1)\nabla_1\frac{1}{r_{12}}\right)\right]\cdot d\mathbf{r}_1$$
$$=\frac{\mu_0}{4\pi}\oint_{L_2}\oint_{L_1}\left[\nabla_1\cdot d\mathbf{r}_2\nabla_1\frac{1}{r_{12}}\right]\cdot d\mathbf{r}_1
=\frac{\mu_0}{4\pi}\oint_{L_2}\oint_{L_1}\left[\nabla_1\nabla_1 \left(\frac{1}{r_{12}}\right)\right]: d\mathbf{r}_1d\mathbf{r}_2
=\frac{\mu_0}{4\pi}\oint_{L_2}\oint_{L_1}\left[\nabla_2\nabla_2 \left(\frac{1}{r_{21}}\right)\right]: d\mathbf{r}_1d\mathbf{r}_2 .$$
The last step is readily shown identical to the result obtained the other way around, in other words, one finds
$$M_{12}=M_{2\to 1}=M_{1\to 2}=M_{21} .$$
