Radial dependence

The form of the solution for \(R(r)\) is easy to guess since the derivatives are balanced by the powers of \(r\), and therefore, a function of the form \(r^l\) will preserve its form up to a coefficient.

\[\begin{eqnarray} \frac{d}{d r}\left(r^2\frac{d r^l }{d r}\right)-c r^l=\left[ l(l+1)-c\right]r^l=0 \implies c=l(l+1). \tag{19} \end{eqnarray}\] However, notice the unexpected symmetry of \(l(l+1)\) under \(l\to-l-1\). This means, if \(r^l\) is a solution, so is \(r^{-l-1}\). This suggests the following form of solution for \(R\): \[\begin{eqnarray} R(r)=a_l r^l +\frac{b_l}{r^{l+1}}. \tag{20} \end{eqnarray}\]

Angular part

Putting this back in Eq. (18) yields \[\begin{eqnarray} \frac{1}{\sin\theta}\frac{d }{d\theta}\left(\sin\theta \frac{ d \Theta}{d\theta}\right) +\left[l(l+1)-\frac{m^2}{\sin^2\theta} \right]\Theta =0. \tag{21} \end{eqnarray}\] Now define \(\cos\theta=x\), which gives \(\frac{d}{d\theta}=\frac{dx}{d\theta}\frac{d}{dx}=-\sin\theta \frac{d}{dx}\) and insert this back in Eq.(21): \[\begin{eqnarray} (1-x^2)\frac{d^2 \Theta}{dx^2}-2x\frac{d\Theta}{dx} +\left[l(l+1)-\frac{m^2}{1-x^2} \right]\Theta =0. \tag{22} \end{eqnarray}\] Let’s first attempt to solve this for \(m=0\) using power series expansion[4]:

\[\begin{eqnarray} \Theta(x)=\sum_{k=0}^\infty c_k x^k. \tag{23} \end{eqnarray}\] Inserting this back to Eq. (22) we get:

\[\begin{eqnarray} 0&=&\sum_{k=0}^\infty c_k x^{k-2}k(k-1)-\sum_{k=0}^\infty c_k x^{k}k(k-1)+\sum_{k=0}^\infty c_k x^k \left[l(l+1)-2k\right]\nonumber\\ &=&\sum_{k=2}^\infty c_k x^{k-2}k(k-1)-\sum_{k=0}^\infty c_k x^{k}k(k-1)+\sum_{k=0}^\infty c_k x^k \left[l(l+1)-2k\right]\nonumber\\ &=& \sum_{k=0}^\infty c_{k+2} x^{k}(k+2)(k+1)+\sum_{k=0}^\infty c_k x^k \left[l(l+1)-k(k+1)\right]\nonumber\\ &=& \sum_{k=0}^\infty \left\{ c_{k+2} (k+2)(k+1)- c_k \left[k(k+1)-l(l+1)\right]\right\} x^{k} . \tag{24} \end{eqnarray}\] This implies \[\begin{eqnarray} c_{k+2}=\frac{k(k+1)-l(l+1)}{ (k+2)(k+1)} c_k =\frac{(k-l)(l+k+1)}{ (k+2)(k+1)} c_k , \tag{25} \end{eqnarray}\] which is the recurrence equation for the expansion coefficients.

This is a remarkable equation since it has profound consequences. Earlier in this blog, we looked at the Quantum Harmonic Oscillator and showed that for a similar series expansion to converge, we had to have the energy quantized. In this particular problem, until this point, we have no indication of \(l\) being an integer. But now, we see that it has to be an integer so that the series truncates for \(k>l\) (for every other \(k\)). That’s the first observation.

The second observation is associated with the parity symmetry of the original differential Eq. (22), which is invariant under \(x\to-x\) upto the overall sign. This shows that the solutions will also be eigenstates of the parity operator, i.e., odd and even \(k\) terms should not mix.

We have a couple of ways of terminating the series. The first one is what we have discussed above, i.e., settin \(l\) to an integer \(k^*\), which will zero out every other \(c_k\). The \(c_k\)’s, with \(k>l\), not addressed by this truncation need to be eliminated directly by their root coefficient, \(c_0\) or \(c_1\). To be more specific, take an example \(l=1\). The \(c_k\)’s with odd \(k\) quickly terminate after \(k=1\): \(c_1,0,\cdots\). The even ones will keep growing: \(c_0, \alpha c_0,\beta\alpha c_0,\cdots\). The only way to tame this series is by killing it at its root, i.e., by setting it \(a_0=0\) so that all the even terms drop out. This shows that even odd powers of \(x\) will not mix preserving respecting the parity symmetry of the original equation.

From Eq. (25), we can explicitly write the fist few Legendre polynomials: \[\begin{eqnarray} P_0(x)&=&1,\nonumber\\ P_1(x)&=&x,\nonumber\\ P_2(x)&=&\frac{1}{2}(3x^2-1),\nonumber\\ P_3(x)&=&\frac{1}{2}(5x^3-3x),\nonumber\\ P_4(x)&=&\frac{1}{8}(35x^4-30x^2+3). \tag{26} \end{eqnarray}\]

Having shown that Legendre polynomials solve the differential equation (with \(m=0\)), \[\begin{eqnarray} (1-x^2)\frac{d^2 P_l}{dx^2}-2x\frac{dP_l}{dx} +l(l+1)P_l =0. \tag{27} \end{eqnarray}\] We now need to address the full equation with \(m\neq 0\). The idea would be to differentiate \(m\) times to create the \(m^2\) term. For this, we will need the Leibniz’s formula: \[\begin{eqnarray} \frac{d^m}{dx^m}\left[f(x)g(x)\right]=\sum_{k=0}^n {n \choose k} \frac{d^k f}{dx^k}\frac{d^{n-k}g}{dx^{n-k}} \tag{28} \end{eqnarray}\] Let’s dive into the differentiation: \[\begin{eqnarray} 0&=&\frac{d^m}{dx^m}\left[(1-x^2)\frac{d^2 P_l}{dx^2}-2x\frac{dP_l}{dx} +l(l+1) P_l\right]\nonumber\\ &=&(1-x^2)\frac{d^m}{dx^m}\frac{d^2 P_l}{dx^2}+m\frac{d}{dx}(1-x^2)\frac{d^{m-1}}{dx^{m-1 }}\frac{d^2 P_l}{dx^2} +\frac{m(m-1)}{2}\frac{d^2}{dx^2}(1-x^2)\frac{d^{m-2}}{dx^{m-2 }}\frac{d^2 P_l}{dx^2} \nonumber\\ &&-2x\frac{d^m}{dx^m}\frac{dP_l}{dx}-2m\frac{d}{dx}( x) \frac{d^{m-1} }{dx^{m-1} }\frac{dP_l}{dx} +l(l+1) \frac{d^m}{dx^m} P_l\nonumber\\ &=&(1-x^2)u''-2mx u' -m(m-1)u-2xu' -2mu +l(l+1)u\nonumber\\ &=&(1-x^2)u''-2(m+1)x u' - m(m+1)u +l(l+1)u\nonumber\\ &=&(1-x^2)u''-2(m+1)x u' - (l-m)(l+m+1)u \tag{29}, \end{eqnarray}\] where \(u\equiv\frac{d^m P_l}{dx^m}\). We still need to modify the equation further so that it matches Eq.(22). First of all, we note that the equation we want to get at was self-adjoint, and we kind of destroyed it as we acted with \(\frac{d^m}{dx^m}\). Let’s restore it and see where it takes us.

Sturm–Liouville theory

We are going to use some machinery from Sturm–Liouville theory on second order differential equations. Consider the second order differential operator \(\mathcal{L}\)[4]:

\[\begin{eqnarray} \mathcal{L} u=\left(p_0(x)\frac{d^2 }{dx^2}+p_1(x) \frac{d }{dx} +p_2(x)\right)u \tag{30}. \end{eqnarray}\] We are going to define the inner product in the function space as an integral in a range \([a,b]\). \[\begin{eqnarray} \langle u| \mathcal{L}| u\rangle\equiv \int_a^b dx u(x) \mathcal{L} u(x)= \int_a^b dx u(x)\left(p_0\frac{d^2 }{dx^2}+p_1 \frac{d }{dx} +p_2\right)u= \int_a^b dx \left( p_0 u u'' + p_1 u u' +u^2 p_2\right) \tag{31}. \end{eqnarray}\] We can integrate Eq.(31) by parts. Let’s look at each term one by one:

\[\begin{eqnarray} u p_0 u'' &=& u \frac{d^2 }{dx^2}(p_0 u)-\frac{d}{dx}(up'_0 u)+2 u p'_0u' \nonumber\\ u p_1 u' &=& - u \frac{d }{dx}(p_1 u)+\frac{d}{dx}(u p_1 u) \tag{32}. \end{eqnarray}\] Putting this back in Eq.(31) gives:

\[\begin{eqnarray} \langle u| \mathcal{L}| u\rangle &=& \left[u(p_1-p_0')u \right]_a^b+ \int_a^b dx u\left(\frac{d^2 }{dx^2}(p_0 u) -\frac{d }{dx}(p_1 u) +p_2 u\right)\nonumber\\ &=&\left[u(p_1-p_0')u \right]_a^b+ \int_a^b dx u \bar{\mathcal{L}}u, \tag{33} \end{eqnarray}\] where the adjoint operator \(\bar{\mathcal{L}}\) is defined as: \[\begin{eqnarray} \bar{\mathcal{L}}u=\frac{d^2 }{dx^2}(p_0 u) -\frac{d }{dx}(p_1 u) +p_2 u. \tag{34} \end{eqnarray}\] Although \(\bar{\mathcal{L}}u\) looks pretty different from \({\mathcal{L}}u\) in Eq. (30), they can actually be the same if \(p_1=p'_0\). Such \(\mathcal{L}\) operators are self-adjoint. Furthermore, note that the boundary term also drops out for self-adjoint operators.

The good news is that if an equation is not self adjoint, it can be converted into that form if it gets multiplied by the following factor: \[\begin{eqnarray} \frac{1}{p_0(x)}\text{exp}\left\{\int^x dt \frac{p_1(t)}{p_0(t)}\right\}. \tag{35} \end{eqnarray}\]

Let’s revisit Eq. (29) to find the factor that will make the equation self-adjoint: \[\begin{eqnarray} \frac{1}{1-x^2}\text{exp}\left\{-\int^x dt \frac{2(m+1)t}{1-t^2}\right\}= \frac{1}{1-x^2}\text{exp}\left\{(m+1)\int^x \frac{d(1-t^2)}{1-t^2}\right\} =(1-x^2)^{m}. \tag{36} \end{eqnarray}\] We will take this factor, multiply Eq. (29) with it to get:

\[\begin{eqnarray} \frac{d}{dx}\left((1-x^2)^{m+1}u'\right) - (l-m)(l+m+1)u=0 \tag{37}. \end{eqnarray}\] Finally, we will want to absord half power of that coefficient into \(u\) by defining \[\begin{eqnarray} v(x)=(1-x^2)^{\frac{m}{2}} u(x) \iff u(x)=(1-x^2)^{-\frac{m}{2}} v(x) \tag{38}, \end{eqnarray}\] to get \[\begin{eqnarray} u'&=&\left[ v' +\frac{m v x}{1-x^2}\right](1-x^2)^{-\frac{m}{2}},\nonumber\\ u''&=&\left[ v'' +\frac{2m v' x}{1-x^2}+\frac{m x}{1-x^2}+\frac{m (m+2) x^2 v}{(1-x^2)^2} \right](1-x^2)^{-\frac{m}{2}} \tag{39}. \end{eqnarray}\] The new function \(v\) satisfies the following equation:

\[\begin{eqnarray} (1-x^2)\frac{d^2 v}{dx^2}-2x\frac{d v}{dx} +\left(l(l+1)-\frac{m^2}{1-x^2} \right) v =0 \tag{40}, \end{eqnarray}\] which is identical to Eq. (22). In conclusion, the angular part of the solution is given by the associated Legendre polynomials as below: \[\begin{eqnarray} P_l^m(x)=(1-x^2)^{\frac{m}{2}}\frac{d^m }{dx^m} P_l(x) \tag{41}. \end{eqnarray}\]

Full solution

Note that the highest power in \(P_l\) is \(l\), and for \(m>l\), we run out of \(x\)’s to differentiate. This automatically limits \(|m|\) to \(l\). Putting all pieces together, the full solution to the Laplace equation reads:

\[\begin{eqnarray} \psi(r,\theta,\phi)=R(r)\Theta(\theta)\Phi(\phi)= \sum_{l=0}^\infty \sum_{m=-l}^l \left(a_{lm} r^l +\frac{b_{lm}}{r^{l+1}}\right) P_l^m(\cos\theta)e^{im \phi}. \tag{42} \end{eqnarray}\] In the case of symmetric coils, only \(m=0\) will contribute. Furthermore, for the solution inside the coil, we need to set \(b_l=0\) so that the solution stays finite: \[\begin{eqnarray} \psi(r,\theta,\phi)=\sum_{l=0}^\infty a_{l} r^l P_l(\cos\theta). \tag{43} \end{eqnarray}\]

The magneto-static potential becomes a sum of Legendre polynomials[5]: \[\begin{equation} \psi(r,\theta)=a_{0}+a_{1}rP_{1}(\cos \theta)+a_{2}r^{2}P_{2}(\cos \theta)+a_{3}r^{3}P_{3}(\cos \theta)+ \ldots . \tag{44} \end{equation}\] We can evaluate Eq. (44) on the \(z-\)axis, i.e., \(\cos \theta =1\) and \(r=z\), to get: \[\begin{equation} \psi(z)=a_{0}+a_{1}z+a_{2}z^{2}+a_{3}z^{3}+ \cdots , \tag{45} \end{equation}\] which should be a familiar form since it is also the Taylor expansion. The trick here is to calculate the coefficients by matching the corresponding magnetic field with the expression we calculated earlier in Eq. (9). We first need to calculate the vector potential from Eq. (9):

\[\begin{eqnarray} \psi(z)&=&-\int^z d\tilde z B( \tilde z) =-\int^z d\tilde z\frac{ N \mu_0 I }{ 2 H }\left( \frac{H/2 -\tilde z}{\sqrt{ R^2 +(H/2-\tilde z)^2}}+\frac{H/2 +\tilde z}{\sqrt{ R^2 +(H/2+\tilde z)^2} } \right)\nonumber\\ &=&-\frac{ N \mu_0 I }{ 2 H }\left(\sqrt{ R^2 +(H/2+\tilde z)^2} -\sqrt{ R^2 +(H/2-\tilde z)^2}\right) =-\frac{ N \mu_0 I }{ \sqrt{ 4 R^2 +H^2} } z +\frac{8 N R^2 \mu_0 I }{ ( 4 R^2 +H^2)^{5/2} } z^3+\cdots \tag{46}. \end{eqnarray}\] Comparing the expansion above with Eq. (45) we can read \(a_0=0\), \(a_1=-\frac{ N \mu_0 I }{ \sqrt{ 4 R^2 +H^2} }\), \(a_2=0\), and \(a_3=\frac{8 N R^2 \mu_0 I }{ ( 4 R^2 +H^2)^{5/2} }\). Putting these back into Eq. (44) along with the Legendre polynomials from Eq. (26), we get ²:

\[\begin{equation} \psi(r,\theta)\simeq -\frac{ N \mu_0 I }{ \sqrt{ 4 R^2 +H^2} } \cos \theta\, r+\frac{4 N R^2 \mu_0 I }{ ( 4 R^2 +H^2)^{5/2} }(5\cos^3\theta-3\cos\theta)\,r^3. \tag{47} \end{equation}\]

Now we can use Eq.(12) to compute the components of the magnetic field:

\[\begin{eqnarray} B_{r} &=& -\frac{\partial\psi(r,\theta)}{\partial r} \simeq \frac{ N \mu_0 I }{ \sqrt{ 4 R^2 +H^2} } \cos \theta-\frac{12 N R^2 \mu_0 I }{ ( 4 R^2 +H^2)^{5/2} }(5\cos^3\theta-3\cos\theta)\,r^2,\nonumber\\ B_{\theta} &=& -\frac{1}{r}\frac{\partial\psi(r,\theta)}{\partial\theta} \simeq -\frac{ N \mu_0 I }{ \sqrt{ 4 R^2 +H^2} } \sin \theta-\frac{12 N R^2 \mu_0 I }{ ( 4 R^2 +H^2)^{5/2} }(\sin\theta-5\cos^2\theta\sin\theta)\,r^2. \tag{48} \end{eqnarray}\]

If we wanted to switch to cylindrical coordinates we can do that by rotating the vector: \[\begin{eqnarray} \begin{pmatrix} B_z \\ B_\rho \end{pmatrix} = \begin{bmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{bmatrix} \begin{pmatrix} B_{r} \\ B_{\theta} \end{pmatrix}. \tag{49} \end{eqnarray}\]

Single loop

Figure 6: A single loop of current. We are now interested in the field at a generic point \(\mathbf{r}\).

For a single loop sitting at \(z=0\) with radius \(r=R\), it is convenient to work in spherical coordinates. \[\begin{eqnarray} \mathbf{J}_s(\mathbf{r}') &=& \lambda \delta(r'-R)\delta(\cos \theta')\hat{\mathbf{\varphi}'} =\lambda \delta(r'-R)\delta(\cos \theta')\left(\cos\varphi' \hat{\mathbf{j}}-\sin\varphi'\hat{\mathbf{i}} \right) \tag{52}, \end{eqnarray}\] where \(\lambda\) is the current density. In order to calculate \(\lambda\) for a loop of wire carrying a current \(I\), let’s intercept the loop with an area perpendicular to it. We can do that by selecting an area on, say, positive \(x\) axis ( i.e., \(\varphi'=0\)), pointing along the \(y\) axis, i.e., \(d\mathbf{S}'= dS' \hat{\mathbf{j}}=r'dr'd\theta' \hat{\mathbf{j}}\). Integrating the current density on this area we should get the total current: \[\begin{eqnarray} \int_S d\mathbf{S}' \cdot \mathbf{J}(\mathbf{r}) &=& \int_S r'dr'd\theta' \lambda \delta(r'-R)\delta(\cos \theta')\hat{\mathbf{j}}\cdot\left(\cos\varphi' \hat{\mathbf{j}}-\sin\varphi'\hat{\mathbf{i}} \right)\bigg\rvert_{\varphi'=0}=R\lambda =I \implies \lambda=I/R \tag{53}. \end{eqnarray}\] Therefore, the properly normalized current is \[\begin{eqnarray} \mathbf{J}_s(\mathbf{r}') &=& \frac{I}{R} \delta(r'-R)\delta(z'-z_s)\hat{\mathbf{\varphi}'} =\frac{I}{R}\delta(r'-R)\delta(\cos \theta')\left(\cos\varphi' \hat{\mathbf{j}}-\sin\varphi'\hat{\mathbf{i}} \right) \tag{54}. \end{eqnarray}\] Before we continue the calculation for a generic observation point \(\mathbf{r}\), let me throw you a curve-ball. Consider the vector potential on the \(z\) axis, i.e., \(\mathbf{r}=\mathbf{z}\). In this case, \(| \mathbf {r}- \mathbf {r}'|\) across the integration as can be seen from Fig. 6. Therefore the integration in Eq. (51) becomes: \(\propto \int_0^{2\pi}d\varphi \hat{\mathbf{\varphi}'}=0\). So the vector potential vanishes on the \(z\) axis. Since \(\mathbf{B}=\mathbf{\nabla}\times \mathbf{A}\), the magnetic field should vanish too! But, we calculated the magnetic field on the \(z\) axis in Eq. (8) and it is clearly non-zero. What went wrong? Here is the problem: if you are cooking a function to take the derivative later, that function has to be defined in the vicinity of the points you will take derivatives. More specifically, you can’t compute \(\mathbf{A}\) on only the \(z\) axis and expect to be able to take its curl since curl will by definition take its derivative stepping off of the \(z\) axis.

The integral we have to deal with for a single loop is this: \[\begin{eqnarray} \mathbf{A}_s(\mathbf{r}) &=& \frac{ \mu_0 }{ 4 \pi } \int d^3 \mathbf{r'} \frac{\mathbf{J}(\mathbf{r'})}{ | \mathbf {r}- \mathbf {r}'|}=\frac{ \mu_0 I }{ 4 \pi R } \int d^3 \mathbf{r'}\frac{1 }{ | \mathbf {r}- \mathbf {r}'|} \delta(r'-R)\delta(\cos\theta')\left(\cos\varphi' \hat{\mathbf{j}}-\sin\varphi'\hat{\mathbf{i}} \right) \tag{55}, \end{eqnarray}\] where we put the subscript \(s\) to remind us that this is for a single loop. We will parameterize the points on the loop centered at \(z=0\) as \(\mathbf {r}'= r' (\cos\varphi'\hat{\mathbf{i}}+\sin\varphi'\hat{\mathbf{j}})\), and the observation point as \(\mathbf {r}= r\cos\theta \hat{\mathbf {z}}+ r\sin\theta (\cos\varphi\hat{\mathbf{i}}+\sin\varphi\hat{\mathbf{j}})\)

\[\begin{eqnarray} | \mathbf {r}- \mathbf {r}'|&=& \sqrt{r^2\cos^2\theta +(r\sin\theta\cos\varphi -r'\cos\varphi')^2+(r\sin\theta\sin\varphi -r'\sin\varphi')^2}\nonumber\\ &=&\sqrt{r^2+r'^2 -2 r r'\sin\theta\cos(\varphi'-\varphi)} \tag{56}. \end{eqnarray}\] Note that the problem has rotational symmetry. We can rotate our coordinate system such that the observation point sits on \(y=0\), i.e., \(\varphi=0\). Once we are done with the computations, we can rotate the vectors back to general \(\mathbf{r}\) point. So let’s set \(\varphi=0\) in Eq. (56) and rewrite Eq. (55) :

\[\begin{eqnarray} \mathbf{A}_s(\mathbf{r}) &=& \frac{ \mu_0 I }{ 4 \pi R } \int \sin\theta' r'^2 dr' d\varphi' \frac{1 }{ \sqrt{r^2+r'^2 -2 r r'\sin\theta\cos(\varphi'-\varphi)}} \delta(r'-R)\delta(\cos\theta')\left(\cos\varphi' \hat{\mathbf{j}}-\sin\varphi'\hat{\mathbf{i}} \right)\nonumber\\ &=& \frac{ \mu_0 I R}{ 4 \pi } \int_{0}^{2\pi} d\varphi' \frac{\cos\varphi' \hat{\mathbf{j}} }{ \sqrt{r^2+R^2 -2 r R\sin\theta\cos(\varphi'-\varphi)}} -\cancelto{0}{\frac{ \mu_0 I R}{ 4 \pi } \int_{0}^{2\pi} d\varphi' \frac{\sin\varphi'\hat{\mathbf{i}} }{ \sqrt{r^2+R^2 -2 r R\sin\theta\cos(\varphi'-\varphi)}} } \tag{57}, \end{eqnarray}\] where the second term vanishes since the integrand is odd and the integral is evaluated over the full range. Note that we evaluated the integral at \(\varphi=0\), and the resulting potential points in \(\hat{\mathbf{j}}\) direction. For generic \(\varphi\) we can simply rotate the coordinate system about the \(z\) axis by \(\varphi\). In this rotated coordinate system \(\hat{\mathbf{j}}\rightarrow \hat{\mathbf{\varphi}}\). Therefore the vector potential reads: \[\begin{eqnarray} \mathbf{A}_s(\mathbf{r}) &=& \hat{\mathbf{\varphi}}\frac{ \mu_0 I R}{ 4 \pi } \int_{0}^{2\pi} d\varphi' \frac{\cos\varphi' }{ \sqrt{r^2+R^2 -2 r R\sin\theta\cos\varphi'}} \tag{58}. \end{eqnarray}\] Let’s define \(\phi'=\pi-\varphi'\) to get \(\cos\varphi'=-\cos\phi'\) and rewrite Eq. (58) as: \[\begin{eqnarray} \mathbf{A}_s(\mathbf{r}) &=& -\hat{\mathbf{\varphi}}\frac{ \mu_0 I R}{ 4 \pi } \int_{-\pi}^{\pi} d\phi' \frac{\cos\phi' }{ \sqrt{r^2+R^2 +2 r R\sin\theta\cos\phi'}} \tag{59}. \end{eqnarray}\] Let’s also use the half angle formula: \(\cos\phi'=1-2\sin^2\frac{\phi'}{2}\) and reorganize the integral: \[\begin{eqnarray} \mathbf{A}_s(\mathbf{r}) &=& -\hat{\mathbf{\varphi}}\frac{ \mu_0 I R}{ 4 \pi } \frac{1}{\sqrt{r^2+R^2 +2 r R\sin\theta}}\int_{-\pi}^{\pi} d\phi' \frac{1-2\sin^2\frac{\phi'}{2} }{ \sqrt{1 -\frac{4 r R\sin\theta}{r^2+R^2 +2 r R\sin\theta}\sin^2\frac{\phi'}{2}}}\nonumber\\ &\equiv& -\hat{\mathbf{\varphi}}\frac{ \mu_0 I R}{ 4 \pi } \frac{1}{\sqrt{r^2+R^2 +2 r R\sin\theta}}\int_{-\pi}^{\pi} d\phi' \frac{1-2\sin^2\frac{\phi'}{2} }{ \sqrt{1 -k^2\sin^2\frac{\phi'}{2}}} \nonumber\\ &=& -\hat{\mathbf{\varphi}}\frac{ \mu_0 I R}{ 4 \pi } \frac{1}{\sqrt{r^2+R^2 +2 r R\sin\theta}}\int_{-\pi}^{\pi} d\phi' \left[\frac{1 }{ \sqrt{1 -k^2\sin^2\frac{\phi'}{2}}}-2 \frac{\sin^2\frac{\phi'}{2} }{ \sqrt{1 -k^2\sin^2\frac{\phi'}{2}}}\right] \nonumber\\ &=& -\hat{\mathbf{\varphi}}\frac{ \mu_0 I R}{ 4 \pi } \frac{1}{\sqrt{r^2+R^2 +2 r R\sin\theta}}\int_{-\pi}^{\pi} d\phi' \left[\frac{1 }{ \sqrt{1 -k^2\sin^2\frac{\phi'}{2}}}-\frac{2}{k^2}\left( \frac{1 }{ \sqrt{1 -k^2\sin^2\frac{\phi'}{2}}}- \sqrt{1 -k^2\sin^2\frac{\phi'}{2}}\right) \right] \nonumber\\ &=& -\hat{\mathbf{\varphi}}\frac{ \mu_0 I R}{ 4 \pi } \frac{1}{k^2 \sqrt{r^2+R^2 +2 r R\sin\theta}}\int_{-\pi}^{\pi} d\phi' \left[\frac{k^2-2 }{ \sqrt{1 -k^2\sin^2\frac{\phi'}{2}}}+2 \sqrt{1 -k^2\sin^2\frac{\phi'}{2}} \right] \tag{60}, \end{eqnarray}\] where \(k^2=\frac{4 r R\sin\theta}{r^2+R^2 +2 r R\sin\theta}\). Finally, we define \(\zeta'=\phi'/2\) and split the integration into two pieces to pick an overall factor of \(4\) to get: \[\begin{eqnarray} \mathbf{A}_s(\mathbf{r})=\hat{\mathbf{\varphi}}\frac{ \mu_0 I R}{ \pi \sqrt{r^2+R^2 +2 r R\sin\theta} } \frac{(2-k^2)K(k^2)-2E(k^2)}{k^2} \tag{61}, \end{eqnarray}\] where the elliptic integral are defined as follows: \[\begin{eqnarray} K(k^2)&=&\int_0^{\frac{\pi}{2}}\frac{d\theta}{\sqrt{1-k^2\sin^2\theta}},\nonumber\\ E(k^2)&=&\int_0^{\frac{\pi}{2}}d\theta\sqrt{1-k^2\sin^2\theta} \tag{62}. \end{eqnarray}\] The calculation of the magnetic field is straightforward[6]:

\[\begin{eqnarray} B_r&=&\frac{1}{r\sin\theta}\frac{\partial}{\partial\theta}(\sin\theta A_\varphi)=\frac{ \mu_0 I R^2}{ \pi \sqrt{r^2+R^2 +2 r R\sin\theta} (r^2+R^2 -2 r R\sin\theta)}\cos\theta E(k^2) ,\nonumber\\ B_\theta&=&-\frac{1}{r}\frac{\partial}{\partial r}(r A_\varphi)=\frac{ \mu_0 I }{ 2 \pi \sqrt{r^2+R^2 +2 r R\sin\theta} (r^2+R^2 -2 r R\sin\theta)\sin\theta}\left[(r^2+R^2\cos(2\theta) )E(k^2)-(r^2+R^2 -2 r R\sin\theta)K(k^2)\right] \tag{63}. \end{eqnarray}\]

We can also express the magnetic field in the cylindrical coordinates [6]: \[\begin{eqnarray} B_\rho&=&\frac{ \mu_0 I z}{ 2 \pi \sqrt{R^2+\rho^2+z^2 +2 R \rho} (R^2+\rho^2+z^2 -2 R \rho)\rho}\left[ (R^2+\rho^2+z^2)E(k^2)-(R^2+\rho^2+z^2 -2 R \rho))K(k^2)\right] ,\nonumber\\ B_z&=&\frac{ \mu_0 I }{ 2 \pi \sqrt{R^2+\rho^2+z^2 +2 R } (R^2+\rho^2+z^2 -2 R \rho)\rho}\left[ (R^2-\rho^2-z^2)E(k^2)+(R^2+\rho^2+z^2 -2 R \rho))K(k^2)\right] \tag{64}. \end{eqnarray}\]

Solenoid

Figure 7: A loop sitting at \(z=z_s\). We will integrate over these loops along the \(z\) direction from \(-L/2\) to \(L/2\).

For the solenoid of \(H\) centered at \(z=0\), the current on the wires at \(\mathbf{r}'=R\) and \(-\frac{H}{2}<z<\frac{H}{2}\), we can express it as follows: \[\begin{eqnarray} \mathbf{J}_s(\mathbf{r}') &=& \lambda \delta(r'-R)\hat{\mathbf{\varphi}'} =\lambda \delta(r'-R)\left(\cos\varphi' \hat{\mathbf{j}}-\sin\varphi'\hat{\mathbf{i}} \right), \text{ for } -\frac{H}{2}<z<\frac{H}{2} \tag{65}, \end{eqnarray}\] where \(\lambda\) is the current density. In order to calculate \(\lambda\) for a loop of wire carrying a current \(I\), let’s intercept the loop with an area perpendicular to it. The surface is a thin rectangle just around \(\mathbf{r}'=R\) and \(-\frac{H}{2}<z<\frac{H}{2}\). We can still select an area on positive \(x\) axis ( i.e., \(\varphi'=0\)), pointing along the \(y\) axis, i.e., \(d\mathbf{S}'= dS' \hat{\mathbf{j}}=dr' H \hat{\mathbf{j}}\). Integrating the current density on this area we should get the total current: \[\begin{eqnarray} \int_S d\mathbf{S}' \cdot \mathbf{J}(\mathbf{r}) &=& \int_S dr'H \lambda \delta(r'-R)\hat{\mathbf{j}}\cdot\left(\cos\varphi' \hat{\mathbf{j}}-\sin\varphi'\hat{\mathbf{i}} \right)\bigg\rvert_{\varphi'=0}=\lambda H =N I \implies \lambda=\frac{NI}{H} \tag{66}, \end{eqnarray}\] where \(N\) is the total number of turns in the coil. Therefore, the properly normalized current is \[\begin{eqnarray} \mathbf{J}_s(\mathbf{r}') &=& \frac{NI}{H} \delta(r'-R)\hat{\mathbf{\varphi}'} =nI\delta(r'-R)\left(\cos\varphi' \hat{\mathbf{j}}-\sin\varphi'\hat{\mathbf{i}} \right) \tag{67}, \end{eqnarray}\] where \(n=\frac{N}{H}\) is the number of turns per unit length. The integral we have to deal with for the coils reads: \[\begin{eqnarray} \mathbf{A}(\mathbf{r}) &=& \frac{ \mu_0 }{ 4 \pi } \int d^3 \mathbf{r'} \frac{\mathbf{J}(\mathbf{r'})}{ | \mathbf {r}- \mathbf {r}'|}=\frac{ \mu_0 n I }{ 4 \pi } \int_{-\frac{H}{2}}^{\frac{H}{2}} dz' r'dr' d\varphi' \frac{1 }{ | \mathbf {r}- \mathbf {r}'|} \delta(r'-R)\left(\cos\varphi' \hat{\mathbf{j}}-\sin\varphi'\hat{\mathbf{i}} \right) \tag{68}, \end{eqnarray}\] where we put the subscript \(s\) to remind us that this is for a single loop. We will parameterize the points on the loop centered at \(z=z'\) as \(\mathbf {r}'= z'\,\mathbf{k}+ r' (\cos\varphi'\hat{\mathbf{i}}+\sin\varphi'\hat{\mathbf{j}})\), and the observation point as \(\mathbf {r}= \mathbf {z}+ r (\cos\varphi\hat{\mathbf{i}}+\sin\varphi\hat{\mathbf{j}})\)

\[\begin{eqnarray} | \mathbf {r}- \mathbf {r}'|=|\mathbf {z}+ r (\cos\varphi\hat{\mathbf{i}}+\sin\varphi\hat{\mathbf{j}})-z'\,\mathbf{k}- r' (\cos\varphi'\hat{\mathbf{i}}+\sin\varphi'\hat{\mathbf{j}})|= \sqrt{(z-z')^2 +r^2+r'^2-2rr'\cos(\varphi-\varphi')} \tag{69}. \end{eqnarray}\] Note that the problem has rotational symmetry. We can rotate our coordinate system such that the observation point sits on \(y=0\), i.e., \(\varphi=0\). Once we are done with the computations, we can rotate the vectors back to general \(\mathbf{r}\) point. So let’s set \(\varphi=0\) in Eq. (56) and rewrite Eq. (55) :

\[\begin{eqnarray} \mathbf{A}(\mathbf{r}) &=& \frac{ \mu_0 n }{ 4 \pi } \int r'dr' dz' d\varphi' \frac{1 }{ \sqrt{(z-z')^2 +r^2+r'^2-2rr'\cos\varphi'}} \delta(r'-R)\left(\cos\varphi' \hat{\mathbf{j}}-\sin\varphi'\hat{\mathbf{i}} \right)\nonumber\\ &=& \frac{ \mu_0 n R}{ 4 \pi } \int_{-\frac{H}{2}}^{\frac{H}{2}} dz'\int_{0}^{2\pi} d\varphi' \frac{\cos\varphi' \hat{\mathbf{j}} }{ \sqrt{(z-z')^2 +r^2+R^2-2rR\cos\varphi'}} -\frac{ \mu_0 n R}{ 4 \pi }\int_{-\frac{H}{2}}^{\frac{H}{2}} dz'\cancelto{0}{ \int_{0}^{2\pi} d\varphi' \frac{\sin\varphi'\hat{\mathbf{i}} }{ \sqrt{(z-z')^2 +r^2+R^2-2rR\cos\varphi'}} } \tag{70}, \end{eqnarray}\] where the second term vanishes since the integrand is odd and the integral is evaluated over the full range. Note that we evaluated the integral at \(\varphi=0\), and the resulting potential points in \(\hat{\mathbf{j}}\) direction. For generic \(\varphi\) we can simply rotate the coordinate system about the \(z\) axis by \(\varphi\). In this rotated coordinate system \(\hat{\mathbf{j}}\rightarrow \hat{\mathbf{\varphi}}\). Therefore the vector potential reads: \[\begin{eqnarray} \mathbf{A}(\mathbf{r}) &=& \hat{\mathbf{\varphi}}\frac{ \mu_0 n R}{ 2 \pi } \int_{0}^{\pi} d\varphi' \int_{-\frac{H}{2}}^{\frac{H}{2}} dz' \frac{\cos\varphi' }{ \sqrt{(z'-z)^2 +r^2+R^2-2rR\cos\varphi'}}, \tag{71}. \end{eqnarray}\] where the integrand periodic and the range can be set to \([0,\pi]\) by inserting an overall factor of \(2\). We will want to do the the \(z'\) integral first by changing the variable as \(z'-z\equiv\sqrt{r^2+R^2-2rR\cos\varphi'} \tan\alpha\):

\[\begin{eqnarray} \mathcal{I}&=& \int dz' \frac{1}{ \sqrt{(z'-z)^2 +r^2+R^2-2rR\cos\varphi'}}= \int d\alpha \sqrt{r^2+R^2-2rR\cos\varphi'} \sec^2\alpha \frac{1}{ \sqrt{r^2+R^2-2rR\cos\varphi'} \sqrt{\tan^2\theta +1}}\nonumber\\ &=&\int d\alpha \sec\alpha =\int d\alpha \sec\alpha \frac{\sec\alpha+\tan \alpha}{\sec\alpha+\tan \alpha} =\int d\alpha \frac{\sec^2\alpha+\tan \alpha\sec \alpha}{\sec\alpha+\tan \alpha} =\int \frac{d(\sec\alpha+\sec \alpha)}{\sec\alpha+\tan \alpha}=\ln|\sec\alpha+\tan \alpha|\nonumber\\ &=&\ln\left| \frac{z'-z+\sqrt{(z'-z)^2+r^2+R^2-2rR\cos\varphi'}}{\sqrt{r^2+R^2-2rR\cos\varphi'}} \right|\nonumber\\ &=&\ln\left| z'-z+\sqrt{(z'-z)^2+r^2+R^2-2rR\cos\varphi'} \right|-\ln\left|\sqrt{r^2+R^2-2rR\cos\varphi'} \right| \tag{72} \end{eqnarray}\] where we can drop the last term since it has no \(z\) dependence and its values at the upper and lower boundary is the same. Let’s define \(z'-z=\xi\), \(\xi_\pm=z\pm\frac{H}{2}\) and go back to the \(\varphi'\) integral:

\[\begin{eqnarray} \mathbf{A}(\mathbf{r}) &=& \hat{\mathbf{\varphi}}\frac{ \mu_0 n R}{ 2 \pi } \int_{0}^{\pi} d\varphi' \cos\varphi' \ln\left| \xi+\sqrt{\xi^2+r^2+R^2-2rR\cos\varphi'} \right|_{\xi_-}^{\xi_+}\nonumber\\ &=& \hat{\mathbf{\varphi}}\frac{ \mu_0 n R}{ 2 \pi } \int_{0}^{\pi} d\left\{\sin\varphi' \ln\left| \xi+\sqrt{\xi^2+r^2+R^2-2rR\cos\varphi'} \right|_{\xi_-}^{\xi_+}\right\}\nonumber\\ &&- \hat{\mathbf{\varphi}}\frac{ \mu_0 n R}{ 2 \pi } \int_{0}^{\pi} d\varphi' \sin\varphi' \frac{d}{d\varphi'}\left\{\ln\left| \xi+\sqrt{\xi^2+r^2+R^2-2rR\cos\varphi'} \right|_{\xi_-}^{\xi_+}\right\}\nonumber\\ &=&\hat{\mathbf{\varphi}}\frac{ \mu_0 n R}{ 2 \pi } \cancelto{0}{\left.\left\{\sin\varphi' \ln\left| \xi+\sqrt{\xi^2+r^2+R^2-2rR\cos\varphi'} \right|_{\xi_-}^{\xi_+}\right\}\right\rvert_0^\pi}\nonumber\\ &&- \hat{\mathbf{\varphi}}\frac{ \mu_0 n R}{ 2 \pi } \int_{0}^{\pi} d\varphi' \left.\frac{rR\sin^2\varphi'}{ \left(\xi+\sqrt{\xi^2+r^2+R^2-2rR\cos\varphi'}\right)\sqrt{\xi^2+r^2+R^2-2rR\cos\varphi'} }\right\rvert_{\xi_-}^{\xi_+} \nonumber\\ &=&-\hat{\mathbf{\varphi}}\frac{ \mu_0 n R}{ 2 \pi } \int_{0}^{\pi} d\varphi' \left.\frac{rR\sin^2\varphi'}{ \left(\xi+\sqrt{\xi^2+r^2+R^2-2rR\cos\varphi'}\right)\sqrt{\xi^2+r^2+R^2-2rR\cos\varphi'} }\right\rvert_{\xi_-}^{\xi_+} \tag{73}. \end{eqnarray}\]

Multiply and divide by \(\sqrt{\xi^2+r^2+R^2-2rR\cos\varphi'}-\xi\) [7] to get: \[\begin{eqnarray} \mathbf{A}(\mathbf{r}) &=&-\hat{\mathbf{\varphi}}\frac{ \mu_0 n I R}{ 2 \pi } \int_{0}^{\pi} d\varphi' \left.\frac{(\sqrt{\xi^2+r^2+R^2-2rR\cos\varphi'}-\xi)rR\sin^2\varphi'}{ \left(\xi^2+r^2+R^2-2rR\cos\varphi'-\xi^2\right)\sqrt{\xi^2+r^2+R^2-2rR\cos\varphi'} }\right\rvert_{\xi_-}^{\xi_+} \nonumber\\ &=&-\hat{\mathbf{\varphi}}\frac{ \mu_0 n I R}{ 2 \pi } \int_{0}^{\pi} d\varphi' \left.\frac{(\sqrt{\xi^2+r^2+R^2-2rR\cos\varphi'}-\xi)rR\sin^2\varphi'}{ \left(r^2+R^2-2rR\cos\varphi'\right)\sqrt{\xi^2+r^2+R^2-2rR\cos\varphi'} }\right\rvert_{\xi_-}^{\xi_+} \nonumber\\ &=&-\hat{\mathbf{\varphi}}\frac{ \mu_0 n I R}{ 2 \pi } \int_{0}^{\pi} d\varphi' \left[\cancelto{0}{\left.\frac{rR\sin^2\varphi'}{ \left(r^2+R^2-2rR\cos\varphi'\right) }\right\rvert_{\xi_-}^{\xi_+}} -\left.\frac{\xi rR\sin^2\varphi'}{ \left(r^2+R^2-2rR\cos\varphi'\right)\sqrt{\xi^2+r^2+R^2-2rR\cos\varphi'} }\right\rvert_{\xi_-}^{\xi_+}\right] \nonumber\\ &=&\hat{\mathbf{\varphi}}\frac{ \mu_0 n I R^2 r}{ 2 \pi } \int_{0}^{\pi} d\varphi' \left.\frac{\xi \sin^2\varphi'}{ \left(r^2+R^2-2rR\cos\varphi'\right)\sqrt{\xi^2+r^2+R^2-2rR\cos\varphi'} }\right\rvert_{\xi_-}^{\xi_+} \tag{74}, \end{eqnarray}\] where the first term vanishes since it has no \(\xi\) dependence. As we did in the Single loop section, Let’s define \(\phi'=\pi-\varphi'\) to get \(\cos\varphi'=-\cos\phi'\). This will change the integral limits from \([0,\pi]\) to \([\pi,0]\) and the measure from \(d\varphi\) to \(-d\phi\). We can flip the integral limits with the negative sign of the measure to get:

\[\begin{eqnarray} \mathbf{A}(\mathbf{r}) &=&\hat{\mathbf{\varphi}}\frac{ \mu_0 n I R^2 r}{ 2 \pi } \int_{0}^{\pi} d\phi' \left.\frac{\xi \sin^2\phi'}{ \left(r^2+R^2+2rR\cos\phi'\right)\sqrt{\xi^2+r^2+R^2+2rR\cos\phi'} }\right\rvert_{\xi_-}^{\xi_+} \tag{75}. \end{eqnarray}\]

we define \(\zeta'=\frac{\phi'}{2}\). Noting \[\begin{eqnarray} \sin^2\phi'=\sin^2(2\zeta')=[2\sin\zeta'\cos\zeta']^2=4\sin^2\zeta'\cos^2\zeta'=4\sin^2\zeta'-4\sin^4\zeta' \tag{76}, \end{eqnarray}\] we get: \[\begin{eqnarray} \mathbf{A}(\mathbf{r}) &=&\hat{\mathbf{\varphi}}\frac{ \mu_0 n I R^2 r}{ \pi } \int_{0}^{\pi/2} d\zeta' \left.\frac{4\sin^2\zeta'-4\sin^4\zeta'}{ \left(r^2+R^2+2rR\cos(2\zeta')\right)\sqrt{\xi^2+r^2+R^2+2rR\cos(2\zeta')} }\xi\right\rvert_{\xi_-}^{\xi_+}\nonumber\\ &=&\hat{\mathbf{\varphi}}\frac{ 4 \mu_0 n I R^2 r}{ \pi } \int_{0}^{\pi/2} d\zeta' \left.\frac{\sin^2\zeta'-\sin^4\zeta'}{ \left(r^2+R^2+2rR\cos(2\zeta')\right)\sqrt{\xi^2+r^2+R^2+2rR\cos(2\zeta')} }\xi\right\rvert_{\xi_-}^{\xi_+} \nonumber\\ &=&\hat{\mathbf{\varphi}}\frac{ 4 \mu_0 n I R^2 r}{ \pi } \int_{0}^{\pi/2} d\zeta' \left.\frac{\sin^2\zeta'-\sin^4\zeta'}{ \left(r^2+R^2+2rR\cos(2\zeta')\right)\sqrt{\xi^2+r^2+R^2+2rR\cos(2\zeta')} }\xi\right\rvert_{\xi_-}^{\xi_+} \nonumber\\ &=&\hat{\mathbf{\varphi}}\frac{4 \mu_0 n I R^2 r}{ \pi } \int_{0}^{\pi/2} d\zeta' \left.\frac{\sin^2\zeta'-\sin^4\zeta'}{ \left((r+R)^2-4rR\sin^2\zeta'\right)\sqrt{\xi^2+(r+R)^2-4rR\sin^2\zeta'} }\xi\right\rvert_{\xi_-}^{\xi_+} \nonumber\\ &=&\hat{\mathbf{\varphi}}\frac{4 \mu_0 n I R^2 r}{ \pi (r+R)^2 \sqrt{\xi^2+(r+R)^2} } \int_{0}^{\pi/2} d\zeta' \left.\frac{\sin^2\zeta'-\sin^4\zeta'}{ \left(1-\frac{4rR}{(r+R)^2}\sin^2\zeta'\right)\sqrt{1-\frac{4rR}{\xi^2+(r+R)^2}\sin^2\zeta'} }\xi\right\rvert_{\xi_-}^{\xi_+} \nonumber\\ &=&\hat{\mathbf{\varphi}}\frac{4 \mu_0 n I R^2 r}{ \pi (r+R)^2 \sqrt{\xi^2+(r+R)^2} } \int_{0}^{\pi/2} d\zeta' \left.\frac{\sin^2\zeta'-\sin^4\zeta'}{ \left(1-h^2\sin^2\zeta'\right)\sqrt{1-k^2\sin^2\zeta'} }\xi\right\rvert_{\xi_-}^{\xi_+} \tag{77}. \end{eqnarray}\] where \(h^2\equiv\frac{4rR}{(r+R)^2}\) and \(k^2\equiv\frac{4rR}{\xi^2+(r+R)^2}\).

We can do partial fractions:

\[\begin{eqnarray} \mathcal{I} &=& \int_{0}^{\pi/2} d\zeta' \frac{\sin^2\zeta'-\sin^4\zeta'}{ \left(1-h^2\sin^2\zeta'\right)\sqrt{1-k^2\sin^2\zeta'}}\nonumber\\ &=& \int_{0}^{\pi/2} d\zeta' \left\{ \frac{1}{h^2}\left[\frac{1}{(1-h^2\sin^2\zeta')\sqrt{1-k^2\sin^2\zeta'}}- \frac{1}{\sqrt{1-k^2\sin^2\zeta'}} \right]+\frac{1}{h^4}\left[ \frac{1+h^2\sin^2\zeta'}{\sqrt{1-k^2\sin^2\zeta'}}- \frac{1}{(1-h^2\sin^2\zeta')\sqrt{1-k^2\sin^2\zeta'}} \right]\right\}\nonumber\\ &=& \frac{k^2+h^2-h^2k^2}{h^2k^2}K(k^2) -\frac{1}{k^2}E(k^2)- \frac{h^2-1}{h^2}\Pi(h^2,k^2) \tag{78}, \end{eqnarray}\] where \[\begin{eqnarray} K(k^2)&=&\int_0^{\frac{\pi}{2}}\frac{d\theta}{\sqrt{1-k^2\sin^2\theta}},\nonumber\\ E(k^2)&=&\int_0^{\frac{\pi}{2}}d\theta\sqrt{1-k^2\sin^2\theta},\nonumber\\ \Pi(h^2,k^2)&=&\int_0^{\frac{\pi}{2}}\frac{d\theta}{(1-h^2\sin^2\theta)\sqrt{1-k^2\sin^2\theta}} \tag{79}. \end{eqnarray}\] Putting it all together, we get:

\[\begin{eqnarray} \mathbf{A}(\mathbf{r}) &=&\left.\hat{\mathbf{\varphi}}\frac{ \mu_0 n I \sqrt{R}}{ 2 \pi\sqrt{r} } \left[\frac{k^2+h^2-h^2k^2}{h^2k^2}K(k^2) -\frac{1}{k^2}E(k^2)- \frac{1-h^2}{h^2}\Pi(h^2,k^2)\right] \xi k\right\rvert_{\xi_-}^{\xi_+} \tag{80}. \end{eqnarray}\]

We can now compute the magnetic field. This will require derivatives of the elliptic integrals[8]:

\[\begin{eqnarray} \frac{dK}{dk}&=&\frac{E}{k(1-k^2)}-\frac{K}{k},\nonumber\\ \frac{dE}{dk}&=&\frac{E}{k}-\frac{K}{k},\nonumber\\ \frac{d\Pi}{dk}&=&\frac{k \Pi}{k^2-h^2}-\frac{kK}{(1-k^2)(h^2-k^2)} \tag{81}. \end{eqnarray}\] Finally we have[9]:

\[\begin{eqnarray} B_r&=& -\partial_z \mathbf{A}(\mathbf{r})_\varphi =\left.\frac{ \mu_0 n I }{ \pi } \sqrt{\frac{R}{r}} \left[\frac{k^2-2}{k}K(k^2) +\frac{2}{k}E(k^2)\right] \xi k\right\rvert_{\xi_-}^{\xi_+} \tag{82}, \end{eqnarray}\] and, \[\begin{eqnarray} B_z&=& \frac{1}{r}\partial_r\left(r \mathbf{A}(\mathbf{r})_\varphi \right) =\left.\frac{ \mu_0 n I }{ 4 r \pi\sqrt{Rr} } \left[K(k^2) +\frac{R-r}{R+r} \Pi(h^2,k^2)\right] \xi k\right\rvert_{\xi_-}^{\xi_+} \tag{83}. \end{eqnarray}\]