Optimal Inverse Design of Magnetic Field Profiles in a Magnetically Shielded Cylinder

Regions of space that require precisely-controlled magnetic field environments are essential for a wide range of experiments, devices, and applications. These include quantum sensing of gravity and gravity gradients for geophysical survey and mapping Wueaax0800 ; doi:10.1038/s41598-018-30608-1 ; SnaddenGradiometer ; AIQuantumSensors ; atomic magnetometry PhysRevLett.113.013001 ; quspin ; 10.1155/2015/491746 for applications including material characterization ROMALIS2011258 and magnetoencephalography 10.1364/BOE.3.000981 ; nature ; noise reduction in fundamental physics experiments doi:10.1063/1.4919366 ; 10.1007/s10751-014-1109-5 such as timing using high-precision atomic clocks MORIC2014287 ; LiangClock , measuring the electric dipole moment of fundamental systems RevModPhys.62.541 ; HindsEDM ; ACMEEDM ; doi:10.1063/1.4922671 , and testing Lorentz-invariance by observing the limits of spin precession PhysRevLett.105.151604 ; PhysRevA.100.010501 ; PhysRevLett.112.110801 . Typically, these systems are enclosed within high-permeability materials, such as mumetal, to shield magnetically sensitive components from undesired magnetic fields. Cylindrical geometries, in particular, are widely used due to their comparatively high shielding factor relative to their low manufacturing cost 10.1063/1.1656455 ; GRABCHIKOV201649 ; 6217348 .

Romeo and Hoult first laid the mathematical framework for the design of magnetic fields in free space through the use of current loops and arcs RandH . These simple elements were expanded in a spherical harmonic basis and unwanted field profiles removed by adjusting the relative separations of the discrete coils. Inverse methods based on formulating a current density in terms of simple triangular boundary elements Pissanetzky_1992 ; doi:10.1002/cmr.b.20091 ; doi:10.1002/cmr.b.20040 enabled the design of coils of arbitrary complexity, giving greater flexibility in the design of systems at the cost of computational power and numerical instability. Pseudo-analytical formulations using a Green’s function expansion and harmonic minimization methods have been developed to generate accurate user-specified magnetic fields in free space Forbes_2002 ; Forbes_2001 ; niall0 to facilitate the need for fast and accurate design procedures. However, these methods did not incorporate the interaction of high-permeability passive shielding material, hindering the accurate generation of specified target magnetic field profiles doi:10.1002/9780470268483.app2 in shielded environments. Consequently, optimization of magnetic field cancellation, or other field shaping, systems in the presence of a material with high magnetic permeability is a long-standing challenge in electromagnetism.

Finite element methods (FEMs) can be used to develop models of hybridized active and passive shielding systems. One method of optimizing active–passive systems using FEMs would be to use genetic algorithms genetic ; genetic1 to evolve the coil parameters iteratively, including the effect of the passive material on the magnetic field generated actively at each iteration. Due to their computational complexity, however, FEMs can be slow, and, if coupled with forward optimization procedures provide only locally optimal designs since desired magnetic field profiles depend highly non-linearly on the system parameters. Analytical formulations of the magnetic field generated by hybrid active–passive systems have the distinct advantage that optimal solutions can be calculated at a range of target points with minimal computation CACIAGLI2018423 . Currently, analytical solutions for coils in high-permeability cylindrical magnetic shields have only been formulated in the specific cases of magnetically-shielded solenoids and discrete current loops Solenoid1 ; solenoid2 ; doi:10.1063/1.1719514 . The geometric parameters of the active structure in these systems, such as the separation distances of wire loops, are adapted to account for the interaction between the active and passive systems. Previously, the method of mirror images jackson has been used to formulate the response of a planar material with high magnetic permeability to a magnetic field generated by a current source cubefem ; niall ; 8500866 ; LIU2020166846 . In these formulations, Maxwell’s equations are solved explicitly by including the reflections of the current source generated by the high-permeability material in order to match the required boundary conditions. More generally, Green’s function solutions to boundary value problems have been calculated for an extensive range of electromagnetic systems greensey , but have not previously been found generally for the case of a finite closed cylindrical high-permeability shield in the presence of an arbitrary cylindrical coaxial static current source.

To address this, here, we incorporate ab initio the effect of a finite length high-permeability cylinder on the magnetic field generated by an arbitrary static current flow on a conducting cylinder, and use our results to determine the flow required to generate specified static target fields optimally. This enables the geometry of the active components on the surface of a cylinder to be determined entirely from a desired magnetic field profile. Guided by turner , we first derive a Green’s function for a hybrid active–passive field-generating system described in cylindrical coordinates. We then utilize a modified Fourier basis to define an adjusted current density distribution which accounts for the effect of the high-permeability material. From this, we determine globally-optimal current density maps using a quadratic optimization procedure, akin to magnetic field design methodologies used previously in Magnetic Resonance Imaging (MRI) doi:10.1002/mrm.1910260202 ; hoult . To construct practical current sources from these current density maps, we define the streamfunction of the current continuum and discretize it into a set of closed-loop current-carrying wire geometries. Using this formulation, we present two example coil designs optimized in the interior of a closed cylindrical magnetic shield of finite length and high permeability, $\mu_{r}\gg 1$, and show that our analytical model agrees well with FEM simulations performed for the optimized current flow patterns. We then show that our design methodology can be used in a real-world situation where the cylindrical magnetic shield has finite thickness and permeability as well as axial entry holes in the end caps to allow experimental access. Using this formulation, we transform the performance of systems designed to generate user-specified magnetic field profiles in the magnetostatic regime, reducing the amount of high-permeability material required, and find globally-optimal solutions required for practical, low Size, Weight, Power, and Cost (SWaP-C) technologies for the applications listed above.

The interface conditions for a magnetic field along a boundary, $S$, between two materials are that the normal component of the magnetic field, $\mathbf{B}$, and tangential component of the magnetic field strength, H, are continuous. Considering the boundary between air and a material, working in the magnetostatic regime, where no eddy currents are induced, and in the case where no surface currents are present, these interface conditions are

$$(\mathbf{B_{\mathrm{mat.}}}-\mathbf{B_{\mathrm{air}}})\cdot\mathbf{\hat{n}}=0\qquad\textnormal{on $S$},$$

(1)

and

$$\left(\frac{1}{\mu_{r}}\mathbf{B_{\mathrm{mat.}}}-\mathbf{B_{\mathrm{air}}}\right)\wedge\mathbf{\hat{n}}=0\qquad\textnormal{on $S$},$$

(2)

where $\hat{\mathbf{n}}$ is the unit vector normal to the boundary, and $\mu_{r}$ is the relative permeability of the material. In free space, the magnetic field is related to the magnetic field strength and the magnetization, M, by

$$\mathbf{B}=\mu_{0}(\mathbf{H}+\textbf{M}).$$

(3)

Physically, the magnetization of the sub-domains of a material change in response to an applied magnetic field to satisfy the boundary condition, (2), at the material’s surface. Here, we choose to formulate this response in terms of a pseudo-current density, $\widetilde{\mathbf{J}}$, confined to the surface of the material, which relates to the curl of the magnetization,

$$\mathbf{\nabla}\wedge\mathbf{M}=\widetilde{\textbf{J}}.$$

(4)

The magnetic field strength resulting from a current source, $\mathbf{J_{c}}$, can be expressed using Ampère’s law

$$\nabla\wedge\mathbf{H}=\mathbf{J_{c}}.$$

(5)

Substituting (4) and (5) into the curl of (3), noting $\mathbf{B}=\mathbf{\nabla}\wedge\mathbf{A}$ where $\mathbf{A}$ is the vector potential, results in the Poisson equation in free space,

$$\nabla^{2}\mathbf{A}=-\mu_{0}(\mathbf{J_{c}}+\widetilde{\textbf{J}}).$$

(6)

As shown in many papers and textbooks jackson , the solution to the Poisson equation, for an arbitrary current density, $\mathbf{J}$, is given by

$$\mathbf{A\left(r\right)}=\mu_{0}\int_{r^{\prime}}\mathrm{d}^{3}\mathbf{r^{\prime}}\ G(\mathbf{r},\mathbf{r^{\prime}})\mathbf{J(r^{\prime})},$$

(7)

where $G(\mathbf{r},\mathbf{r^{\prime}})$ is the associated Green’s function.

Since high-permeability materials, such as mumetal, can have $\mu_{r}$ values >100,000 times that of air, the magnetic field must abruptly change direction at the boundary between its surface and air to satisfy the boundary conditions (1)-(2). When the shield is of sufficient thickness, the tangential components of the magnetic field at its boundary must be approximately zero. The boundary condition, (2), on the interior surface of the exterior closed cylinder for the example depicted in Fig. 1, requires that

$$B_{\rho}\bigg{\rvert}_{z=\pm L_{s}/2}\hskip-30.0pt\approx 0,\quad\quad B_{\phi}\bigg{\rvert}_{z=\pm L_{s}/2,\rho=\rho_{s}}\hskip-51.0pt\approx 0,\quad\hskip 35.0ptB_{z}\bigg{\rvert}_{\rho=\rho_{s}}\approx 0,$$

(8)

where the magnetic field, $\mathbf{B}=B_{\rho}~\bm{\hat{\rho}}+B_{\phi}~\bm{\hat{\phi}}+B_{z}~\mathbf{\hat{z}}$, is expressed in cylindrical polar coordinates. Previously, it has been found that the response of a high-permeability material deviates from that of a perfect magnetic conductor on the scale of $\mathcal{O}\left(\mu_{r}^{-1}\right)$ mu ; mu1 . Therefore, we may assume for high-permeability materials, such as mumetal, the response can be approximated to that of a perfect magnetic conductor without introducing significant errors.

Due to the symmetries of the system it is natural to work in cylindrical coordinates. Following the formulation of Turner turner , we express the vector potential (7) due to an arbitrary current distribution, $\mathbf{J}=J_{\rho}(\mathbf{r^{\prime}})~\bm{\hat{\rho}}+J_{\phi}(\mathbf{r^{\prime}})~\bm{\hat{\phi}}+J_{z}(\mathbf{r^{\prime}})~\hat{\mathbf{z}}$, as

$$A_{\rho}\left(\mathbf{r}\right)=\mu_{0}\int_{r^{\prime}}\mathrm{d}^{3}\mathbf{r^{\prime}}\ G(\mathbf{r},\mathbf{r^{\prime}})\left(J_{\rho}\left(\mathbf{r^{\prime}}\right)\cos\left(\phi-\phi^{\prime}\right)+J_{\phi}(\mathbf{r^{\prime}})\sin\left(\phi-\phi^{\prime}\right)\right),$$

(9)

$$A_{\phi}\left(\mathbf{r}\right)=-\mu_{0}\int_{r^{\prime}}\mathrm{d}^{3}\mathbf{r^{\prime}}\ G(\mathbf{r},\mathbf{r^{\prime}})\left(J_{\rho}\left(\mathbf{r^{\prime}}\right)\sin\left(\phi-\phi^{\prime}\right)-J_{\phi}(\mathbf{r^{\prime}})\cos\left(\phi-\phi^{\prime}\right)\right),$$

(10)

$$A_{z}\left(\mathbf{r}\right)=\mu_{0}\int_{r^{\prime}}\mathrm{d}^{3}\mathbf{r^{\prime}}\ G(\mathbf{r},\mathbf{r^{\prime}})J_{z}(\mathbf{r^{\prime}}).$$

(11)

Since the current has been restricted to flow over an interior cylindrical shell centred radially about the origin, its radial components may be set to zero, resulting in the simplified vector potentials

$$A_{\rho}\left(\mathbf{r}\right)=\mu_{0}\int_{r^{\prime}}\mathrm{d}^{3}\mathbf{r^{\prime}}\ G(\mathbf{r},\mathbf{r^{\prime}})J_{\phi}(\mathbf{r^{\prime}})\sin\left(\phi-\phi^{\prime}\right),$$

(12)

$$A_{\phi}\left(\mathbf{r}\right)=\mu_{0}\int_{r^{\prime}}\mathrm{d}^{3}\mathbf{r^{\prime}}\ G(\mathbf{r},\mathbf{r^{\prime}})J_{\phi}(\mathbf{r^{\prime}})\cos\left(\phi-\phi^{\prime}\right).$$

(13)

Substituting the Green’s function solution from (7) in cylindrical coordinates jackson ,

$$G(\mathbf{r,r^{\prime}})=\frac{1}{4\pi^{2}}\sum_{m=-\infty}^{\infty}e^{im(\phi-\phi^{\prime})}\int_{-\infty}^{\infty}\mathrm{d}k\ e^{ik\left(z-z^{\prime}\right)}I_{m}(|k|\rho_{<})K_{m}(|k|\rho_{>}),$$

(14)

into (11)-(13), the components of the vector potential may be expressed as

	$\displaystyle A_{\rho}\left(\rho,\phi,z\right)=-\frac{i\mu_{0}\rho^{\prime}}{4\pi}\sum_{m=-\infty}^{\infty}\int_{-\infty}^{\infty}\mathrm{d}k\ e^{im\phi}e^{ikz}\Big{[}I_{m-1}(|k|\rho_{<})K_{m-1}(|k|\rho_{>})\qquad\qquad\qquad\qquad\qquad\qquad$
	$\displaystyle-I_{m+1}(|k|\rho_{<})K_{m+1}(|k|\rho_{>})\Big{]}J_{\phi}^{m}(k),$		(15)

	$\displaystyle A_{\phi}\left(\rho,\phi,z\right)=\frac{\mu_{0}\rho^{\prime}}{4\pi}\sum_{m=-\infty}^{\infty}\int_{-\infty}^{\infty}\mathrm{d}k\ e^{im\phi}e^{ikz}\Big{[}I_{m-1}(|k|\rho_{<})K_{m-1}(|k|\rho_{>})\ \ \ \ \qquad\qquad\qquad\qquad\qquad\qquad$
	$\displaystyle+I_{m+1}(|k|\rho_{<})K_{m+1}(|k|\rho_{>})\Big{]}J_{\phi}^{m}(k),$		(16)

$\displaystyle A_{z}\left(\rho,\phi,z\right)=\frac{\mu_{0}\rho^{\prime}}{2\pi}\sum_{m=-\infty}^{\infty}\int_{-\infty}^{\infty}\mathrm{d}k\ e^{im\phi}e^{ikz}I_{m}(|k|\rho_{<})K_{m}(|k|\rho_{>})J_{z}^{m}(k),$

(17)

where $\rho^{\prime}$ is the radius of the cylinder and $\rho_{>}$, $\rho_{<}$ is either $\rho$, $\rho^{\prime}$ if $\rho>\rho^{\prime}$ or $\rho^{\prime}$, $\rho$ if $\rho<\rho^{\prime}$, respectively. Equations (15)-(17) contain Fourier transforms of the current densities,

$$J_{\phi}^{m}(k)=\frac{1}{2\pi}\int_{0}^{2\pi}\mathrm{d}\phi^{\prime}\ e^{-im\phi^{\prime}}\int_{-\infty}^{\infty}\mathrm{d}z^{\prime}\ e^{-ikz^{\prime}}J_{\phi}\left(\phi^{\prime},z^{\prime}\right),$$

(18)

$$J_{z}^{m}(k)=\frac{1}{2\pi}\int_{0}^{2\pi}\mathrm{d}\phi^{\prime}\ e^{-im\phi^{\prime}}\int_{-\infty}^{\infty}\mathrm{d}z^{\prime}\ e^{-ikz^{\prime}}J_{z}\left(\phi^{\prime},z^{\prime}\right),$$

(19)

with their corresponding inverse transforms given by

$$J_{\phi}(\phi^{\prime},z^{\prime})=\frac{1}{2\pi}\sum_{m=-\infty}^{\infty}\int_{-\infty}^{\infty}\mathrm{d}k\ e^{im\phi^{\prime}}e^{ikz^{\prime}}J_{\phi}^{m}\left(k\right),$$

(20)

$$J_{z}(\phi^{\prime},z^{\prime})=\frac{1}{2\pi}\sum_{m=-\infty}^{\infty}\int_{-\infty}^{\infty}\mathrm{d}k\ e^{im\phi^{\prime}}e^{ikz^{\prime}}J_{z}^{m}\left(k\right).$$

(21)

As a result of this formulation, we can now combine methods for matching the boundary conditions at the radial interface, akin to the formulation of passive screening of magnetic field gradients for MRI ap , with the method of mirror images, accounting for the effect of the planar end caps, to determine the effect of the high-permeability material on the overall magnetic field profile. Due to the cylindrical symmetry of the system, the radial boundary condition may be satisfied by matching the azimuthal Fourier modes in the magnetic field, generated by the cylindrical current source, through the use of a pseudo-current density induced on an infinite cylinder. As the planar end caps are spatially orthogonal to the pseudo-current generated by the infinite cylinder, any image current formed by applying the method of mirror images continues to satisfy the radial boundary condition. Consequently, we can decouple the boundary conditions on the high-permeability cylinder and at the planar end cap boundaries. This must be done sequentially to generate mirror images of the pseudo-current density induced on the high-permeability cylindrical shell and, hence, satisfy the boundary conditions over the entire domain of the closed finite high-permeability cylinder. As a result, using (15)-(17) and the usual formulation of the method of mirror images with two infinite parallel planes, depicted in Fig. 2, the vector potential in the region $\rho_{c}\leq\rho\leq\rho_{s}$, including the effect of the high-permeability cylinder, may be written in terms of an infinite summation over the reflected image currents,

	$\displaystyle A_{\rho}\left(\rho,\phi,z\right)=-\frac{i\mu_{0}}{4\pi}\sum_{m=-\infty}^{\infty}\sum_{p=-\infty}^{\infty}\int_{-\infty}^{\infty}\mathrm{d}k\ e^{im\phi}e^{ikz}\Bigg{\{}\rho_{c}\Big{[}I_{m-1}(|k|\rho_{c})K_{m-1}(|k|\rho)\qquad\qquad\qquad\qquad\qquad$
	$\displaystyle-I_{m+1}(|k|\rho_{c})K_{m+1}(|k|\rho)\Big{]}J_{\phi}^{mp}(k)+\rho_{s}\Big{[}I_{m-1}(|k|\rho)K_{m-1}(|k|\rho_{s})\qquad\qquad$
	$\displaystyle-I_{m+1}(|k|\rho)K_{m+1}(|k|\rho_{s})\Big{]}\widetilde{J_{\phi}^{mp}}(k)\Bigg{\}},$		(22)

	$\displaystyle A_{\phi}\left(\rho,\phi,z\right)=\frac{\mu_{0}}{4\pi}\sum_{m=-\infty}^{\infty}\sum_{p=-\infty}^{\infty}\int_{-\infty}^{\infty}\mathrm{d}k\ e^{im\phi}e^{ikz}\Bigg{\{}\rho_{c}\Big{[}I_{m-1}(|k|\rho_{c})K_{m-1}(|k|\rho)\qquad\qquad\qquad\qquad\qquad\ \ \ \$
	$\displaystyle+I_{m+1}(|k|\rho_{c})K_{m+1}(|k|\rho)\Big{]}J_{\phi}^{mp}(k)+\rho_{s}\Big{[}I_{m-1}(|k|\rho)K_{m-1}(|k|\rho_{s})\qquad\qquad$
	$\displaystyle+I_{m+1}(|k|\rho)K_{m+1}(|k|\rho_{s})\Big{]}\widetilde{J_{\phi}^{mp}}(k)\Bigg{\}},$		(23)

	$\displaystyle A_{z}\left(\rho,\phi,z\right)=\frac{\mu_{0}}{2\pi}\sum_{m=-\infty}^{\infty}\sum_{p=-\infty}^{\infty}\int_{-\infty}^{\infty}\mathrm{d}k\ e^{im\phi}e^{ikz}\Big{[}\rho_{c}I_{m}(|k|\rho_{c})K_{m}(|k|\rho)J_{z}^{mp}(k)\qquad\qquad\qquad\qquad\qquad\ \ \ \$
	$\displaystyle+\rho_{s}I_{m}(|k|\rho)K_{m}(|k|\rho_{s})\widetilde{J_{z}^{mp}}(k)\Big{]},$		(24)

where

$$J_{\phi}^{mp}(k)=\frac{1}{2\pi}\int_{0}^{2\pi}\mathrm{d}\phi^{\prime}\ e^{-im\phi^{\prime}}\int_{-\infty}^{\infty}\mathrm{d}z^{\prime}\ e^{-ikz^{\prime}}J_{\phi}^{p}\left(\phi^{\prime},z^{\prime}\right),$$

(25)

$$J_{z}^{mp}(k)=\frac{1}{2\pi}\int_{0}^{2\pi}\mathrm{d}\phi^{\prime}\ e^{-im\phi^{\prime}}\int_{-\infty}^{\infty}\mathrm{d}z^{\prime}\ e^{-ikz^{\prime}}J_{z}^{p}\left(\phi^{\prime},z^{\prime}\right),$$

(26)

are the Fourier transforms of the $p^{th}$ reflected image current and $\widetilde{J^{mp}_{\phi,z}}(k)$ is the Fourier transform of the $p^{th}$ reflected pseudo-current induced on the high-permeability cylinder. Fig. 2 depicts how azimuthal, $J_{\phi}^{p}\left(\phi^{\prime},z^{\prime}\right)$, and axial, $J_{z}^{p}\left(\phi^{\prime},z^{\prime}\right)$, image currents are generated by two parallel planar perfect magnetic conductors.

Writing the magnetic field in cylindrical coordinates,

$$\mathbf{B}=\nabla\wedge\mathbf{A}=\left(\frac{1}{\rho}\frac{\partial A_{z}}{\partial\phi}-\frac{\partial A_{\phi}}{\partial z}\right)\bm{\hat{\rho}}+\left(\frac{\partial A_{\rho}}{\partial z}-\frac{\partial A_{z}}{\partial\rho}\right)\bm{\hat{\phi}}+\frac{1}{\rho}\left(\frac{\partial}{\partial\rho}(\rho A_{\phi})-\frac{\partial A_{\rho}}{\partial\phi}\right)\mathbf{\hat{z}},$$

(27)

the boundary conditions, (8), on the magnetic field at $\rho=\rho_{s}$ are

$\displaystyle\frac{1}{\rho}\left(\frac{\partial}{\partial\rho}(\rho A_{\phi})-\frac{\partial A_{\rho}}{\partial\phi}\right)\Bigg{\rvert}_{\rho=\rho_{s}}=0,$

(28)

$\displaystyle\left(\frac{\partial A_{\rho}}{\partial z}-\frac{\partial A_{z}}{\partial\rho}\right)\Bigg{\rvert}_{\rho=\rho_{s}}=0.$

(29)

Substituting (22)-(24) into (28)-(29), the Fourier transformed pseudo-current density on the cylindrical wall of the high-permeability cylinder is found to be

$$\widetilde{J_{\phi,z}^{mp}}(k)=-\frac{\rho_{c}I^{\prime}_{m}(|k|\rho_{c})K_{m}(|k|\rho_{s})}{\rho_{s}I_{m}(|k|\rho_{s})K^{\prime}_{m}(|k|\rho_{s})}J_{\phi,z}^{mp}(k),$$

(30)

where $I^{\prime}_{m}(z)$ and $K^{\prime}_{m}(z)$ are the derivatives with respect to $z$ of $I_{m}(z)$ and $K_{m}(z)$, respectively. Using the continuity equation, (20)-(24), and (30), the magnetic field interior to the conducting cylinder, for all points $\rho<\rho_{c}$, is given by

$\displaystyle B_{\rho}\left(\rho,\phi,z\right)=\frac{i\mu_{0}\rho_{c}}{2\pi}\sum_{m=-\infty}^{\infty}\sum_{p=-\infty}^{\infty}\int_{-\infty}^{\infty}\mathrm{d}k\ ke^{im\phi}e^{ikz}I^{\prime}_{m}(|k|\rho)R_{m}(k,\rho_{c},\rho_{s})J_{\phi}^{mp}(k),$

(31)

$\displaystyle B_{\phi}\left(\rho,\phi,z\right)=-\frac{\mu_{0}\rho_{c}}{2\pi\rho}\sum_{m=-\infty}^{\infty}\sum_{p=-\infty}^{\infty}\int_{-\infty}^{\infty}\mathrm{d}k\ m\frac{|k|}{k}e^{im\phi}e^{ikz}I_{m}(|k|\rho)R_{m}(k,\rho_{c},\rho_{s})J_{\phi}^{mp}(k),$

(32)

$\displaystyle B_{z}\left(\rho,\phi,z\right)=-\frac{\mu_{0}\rho_{c}}{2\pi}\sum_{m=-\infty}^{\infty}\sum_{p=-\infty}^{\infty}\int_{-\infty}^{\infty}\mathrm{d}k\ |k|e^{im\phi}e^{ikz}I_{m}(|k|\rho)R_{m}(k,\rho_{c},\rho_{s})J_{\phi}^{mp}(k),$

(33)

where $R_{m}(k,\rho_{c},\rho_{s})=K^{\prime}_{m}(|k|\rho_{c})-\frac{I^{\prime}_{m}(|k|\rho_{c})K_{m}(|k|\rho_{s})}{I_{m}(|k|\rho_{s})}$. In order to determine the magnetic field generated by an arbitrary cylindrical current source, we must construct an orthogonal basis defined on a finite cylinder, which accounts for the mirror images. To do this we use a modified Fourier basis, defining the $p^{\textnormal{th}}$ reflected azimuthal current to be

	$\displaystyle J_{\phi}^{p}(\phi^{\prime},z^{\prime})$	$\displaystyle=\left(T^{p}_{e}(z^{\prime};L_{1},L_{2},L_{s})+T^{p}_{o}(z^{\prime};L_{1},L_{2},L_{s})\right)\Bigg{[}\sum_{n=1}^{N}W_{n0}\sin\left(\frac{n\pi\left((-1)^{p}\left(z^{\prime}-pL_{s}\right)-L_{2}\right)}{L_{c}}\right)$
		$\displaystyle\qquad\qquad+\sum_{n=1}^{N}\sum_{m=1}^{M}\left(W_{nm}\cos(m\phi^{\prime})+Q_{nm}\sin(m\phi^{\prime})\right)\cos\left(\frac{n\pi\left((-1)^{p}\left(z^{\prime}-pL_{s}\right)-L_{2}\right)}{L_{c}}\right)\Bigg{]},$		(34)

in which

$$T^{p}_{e}(z^{\prime};L_{1},L_{2},L_{s})=\left(H\left(z^{\prime}-L_{2}-pL_{s}\right)-H\left(z^{\prime}-L_{1}-pL_{s}\right)\right)\left(\frac{1+(-1)^{p}}{2}\right),$$

(35)

$$T^{p}_{o}(z^{\prime};L_{1},L_{2},L_{s})=\left(H\left(z^{\prime}+L_{1}-pL_{s}\right)-H\left(z^{\prime}+L_{2}-pL_{s}\right)\right)\left(\frac{1-(-1)^{p}}{2}\right),$$

(36)

where $H(x)$ is the Heaviside function and $\left(W_{n0},W_{nm},Q_{nm}\right)$ are Fourier coefficients to be determined. Substituting (34) into (31)-(33), we derive a set of governing equations which relate the magnetic field to the set of weighted Fourier coefficients,

$$B_{\rho}\left(\rho,\phi,z\right)=\sum_{n=1}^{N}W_{n0}F_{n}\left(\rho,z\right)+\sum_{n=1}^{N}\sum_{m=1}^{M}\left(W_{nm}G^{w}_{nm}\left(\rho,\phi,z\right)+Q_{nm}G^{q}_{nm}\left(\rho,\phi,z\right)\right),$$

(37)

$$B_{\phi}\left(\rho,\phi,z\right)=\sum_{n=1}^{N}\sum_{m=1}^{M}\left(W_{nm}H^{w}_{nm}\left(\rho,\phi,z\right)+Q_{nm}H^{q}_{nm}\left(\rho,\phi,z\right)\right),$$

(38)

$$B_{z}\left(\rho,\phi,z\right)=\sum_{n=1}^{N}W_{n0}D_{n}\left(\rho,z\right)+\sum_{n=1}^{N}\sum_{m=1}^{M}\left(W_{nm}S^{w}_{nm}\left(\rho,\phi,z\right)+Q_{nm}S^{q}_{nm}\left(\rho,\phi,z\right)\right),$$

(39)

where the functions $F_{n}\left(\rho,z\right)$, $G^{w,q}_{nm}\left(\rho,\phi,z\right)$, $H^{w,q}_{nm}\left(\rho,\phi,z\right)$, $D_{n}\left(\rho,z\right)$, and $S^{w,q}_{nm}\left(\rho,\phi,z\right)$ are defined in Appendix A. Having determined the magnetic field produced by an arbitrary cylindrical current, we may now use an inverse method to solve the system of governing equations, (37)-(39), to determine the unknown Fourier coefficients $\left(W_{n0},W_{nm},Q_{nm}\right)$ for a specified target magnetic field. Following work done by Carlson et al. doi:10.1002/mrm.1910260202 , this may be done by a least squares minimization with the addition of a penalty term to regularize the problem. This regularization term may take many forms, with individual contributions to it representing, for example, the curvature of a given wire geometry, the power consumption, or any other physical parameter that depends quadratically on the geometry of the coil. In this work we focus on the overall power dissipated in the cylindrical current flow, but this choice is somewhat arbitrary since all of the regularization parameters act to achieve the same general goal. If the regularization term is large, the result is a well-conditioned inverse problem that yields a simple current flow, but reduced field fidelity. On the other hand, if the regularization term is small, then the result is a less well conditioned inverse problem that yields a more intricate pattern of current flow but a higher-fidelity magnetic field. The power dissipation in the conducting cylinder of thickness, $t$, and resistivity, $\varrho$, is given by

$$P=\frac{\rho_{c}\varrho}{t}\int_{L_{2}}^{L_{1}}\mathrm{d}z^{\prime}\int_{0}^{2\pi}\mathrm{d}\phi^{\prime}\ |J_{z}(\phi^{\prime},z^{\prime})|^{2}+|J_{\phi}(\phi^{\prime},z^{\prime})|^{2},$$

(40)

which, when integrated over the surface of the cylinder, using the continuity equation and (34), gives

$$P=\frac{\rho_{c}\varrho}{t}\left[\sum_{n=1}^{N}W_{n0}^{2}\pi L_{c}+\sum_{n=1}^{N}\sum_{m=1}^{M}\left(W_{nm}^{2}+Q_{nm}^{2}\right)\left(\frac{\pi L_{c}}{2}+\frac{m^{2}L_{c}^{3}}{2\pi n^{2}\rho_{c}^{2}}\right)\right].$$

(41)

We now construct a cost function using a least squares optimization procedure,

$$\Phi=\sum_{k}\left[\mathbf{B}^{\mathrm{desired}}\left(\mathbf{r}_{k}\right)-\mathbf{B}\left(\mathbf{r}_{k}\right)\right]^{2}+\beta P,$$

(42)

evaluated at $K$ target field points, where $\beta$ is a weighting parameter chosen such that the physical parameters may be adjusted to achieve specified physical constraints. The minimization is achieved by taking the derivative of the cost function with respect to the Fourier coefficients,

$$\frac{\partial\Phi}{\partial W_{i0}}=0,\qquad\frac{\partial\Phi}{\partial W_{ij}}=0,\qquad\frac{\partial\Phi}{\partial Q_{ij}}=0,\qquad i,j>1,$$

(43)

allowing the optimal Fourier coefficients to be found by matrix inversion for any physically attainable target magnetic field. The current density is then related, through the continuity equation, to the streamfunction

	$\displaystyle\varphi(\phi^{\prime},z^{\prime})=\left(H\left(z^{\prime}-L_{1}\right)-H\left(z^{\prime}-L_{2}\right)\right)\Bigg{[}\sum_{n=1}^{N}\frac{L_{c}}{n\pi}W_{n0}\cos\left(\frac{n\pi\left(z^{\prime}-L_{2}\right)}{L_{c}}\right)\qquad\qquad\qquad\qquad\qquad\qquad$
	$\displaystyle-\sum_{n=1}^{N}\sum_{m=1}^{M}\frac{L_{c}}{n\pi}\left(W_{nm}\cos(m\phi^{\prime})+Q_{nm}\sin(m\phi^{\prime})\right)\sin\left(\frac{n\pi\left(z^{\prime}-L_{2}\right)}{L_{c}}\right)\Bigg{]},$		(44)

on the inner conducting cylindrical surface.

Ideally, the Fourier series would have an infinite number of coefficients. However, truncating the series at a large finite number of terms still provides accurate solutions for sufficiently regularized problems. This is because higher-order terms contribute less at distances much larger than their spatial period. When designing fields in regions close to the conducting surface, if a sufficiently large number of coefficients is chosen, the real-world field fidelity is limited by the accuracy to which the current continuum can be approximated by discretized wire configurations, rather than by the number of Fourier coefficients. Typically, setting N=200 and equating M to the degree of the desired field harmonic is adequate. An approximate solution to the current continuum is then found by discretizing the streamfunction into $N_{\varphi}$ values, where $\varphi_{j}=\textnormal{min}\ \varphi+(j-1/2)\Delta\varphi,j=1,...,N_{\varphi}$, separated by $\Delta\varphi=\frac{\textnormal{max}\ \varphi-\textnormal{min}\ \varphi}{N_{\varphi}}$. The number of contours should be maximized according to the constraints in manufacturing because the accuracy of the model depends on the approximation of the continuum. For situations where only coarse discretization is possible, we recommend that a form of discrete optimization is used instead. When physically manufacturing these structures, the distance of the wires from the high-permeability material is important. Without the passive shield, one can determine how accurately the discrete coil represents the current continuum by using the elemental Biot–Savart law to calculate the error as $N_{\varphi}$ is adjusted.

We now analyze our theoretical model by designing and testing hybrid active–passive magnetic field-generating systems. Regarding the validation of our calculations, we first note that, as expected from previous work Solenoid1 ; solenoid2 ; doi:10.1063/1.1719514 and shown in the Supplementary Material, our calculations confirm that the optimal coil design for generating a constant axial field inside a closed cylindrical perfect magnetic shield is a perfect solenoid that runs along the full length of the cylindrical shield.

In Fig. 3 and Fig. 4 respectively, we show active–passive systems for generating a constant transverse field, $B_{x}$, normal to the axis of the cylinder, and a linear transverse field gradient, $\mathbf{B}=(z~\mathbf{\hat{x}}+x~\mathbf{\hat{z}})$, along the axis of the cylinder. Each of these systems has a cylindrical surface of length $L_{c}=0.95$ m and radius $\rho_{c}=0.245$ m, which carries the coil current distribution. This current distribution is interior to and centred about the origin of a closed perfect magnetically conducting cylinder of length $L_{s}=1$ m and radius $\rho_{s}=0.25$ m. The field is optimized over the central cylindrical region spanning half the radius and length of the coil cylinder. Fig. 3a and Fig. 4a show the respective streamfunctions on the surface of the cylinder and their discretization into wire patterns. The magnetic fields shown in Fig. 3b–c and Fig. 4b are calculated in three ways: analytically, using our theoretical model in (37)-(39); numerically, using COMSOL Multiphysics^® with the shield treated as a perfect magnetic conductor; and numerically in free space, i.e. excluding the high-permeability material and calculating the magnetic fields by applying the elemental Biot–Savart law directly to the discrete coil geometries in Fig. 3a and Fig. 4a. It is clear from Figs. 3–4 that our design methodology is capable of generating highly accurate user-specified target magnetic fields inside the optimized field region with good agreement between the theoretical model and numerical simulations.

We quantify this agreement by analyzing the deviation from the target fields in the optimization region, $\Delta{B_{x}}$ and $\Delta{\mathrm{d}B_{x}/\mathrm{d}z}$, for the constant transverse field and linear transverse field gradient systems, respectively. Specifically, along the $z$-axis of the optimized region the maximum absolute deviations from the target fields are $0.11$% and $0.24$%, respectively. Over the same region, the maximum absolute deviations between the numerically-simulated and analytically-calculated field profiles are $0.002$% and $0.003$%, respectively. We can also see the hybrid nature of our optimization by the improved performance of the active systems when inside, and coupled to, the passive shield. For example, the strength of the $B_{x}$ field is nearly doubled, and its uniformity is improved by a factor of $20$, when the high-permeability cylinder is added to the constant transverse field-generating system (see Fig. 3b).

In Fig. 5 we show the numerically-calculated color maps of the $B_{x}$ field in the $y$-$z$ plane generated by the constant transverse field (Fig. 3a) and linear transverse field gradient (Fig. 4a) systems, using COMSOL Multiphysics^® with the shield treated as a perfect magnetic conductor. We summarize the performance of both systems in Table 1, calculating the cylindrical shield fractions – the ratio of the radial and axial extent of the central region to that of the passive shield – where the maximum deviations from the target fields are less than $0.01\%$, $0.05\%$, $0.1\%$, $0.5\%$, $1\%$, and $5\%$, respectively.

When physically constructing active–passive structures for real-world experiments, additional limitations must be taken into consideration in order to generate accurate magnetic fields using our design methodology. These limitations originate either from the theoretical model itself or from experimental practicalities. The limitations in the theoretical model are primarily associated with how accurately the high-permeability cylinder approximates a perfect magnetic conductor. This depends on the value of the finite permeability, the thickness of the shielding material, and the required experimental access holes in the shielding system. The errors introduced by these parameters depend on the lengths, radii, and positions of the conducting and high-permeability cylinders relative to the location of the optimization region. The experimental limitations on the field fidelity relate to the stability of the experimental equipment and errors in manufacturing an accurate representation of the current continuum. These errors include coarse discretization of the current continuum, inexact wire placement and construction, and imprecise positioning of the active structure inside the high-permeability shield. In practice, highly stable experimental equipment is available, particularly power supplies and current drivers Ac_n_2018 , and it is possible to manufacture structures that approximate the current continuum accurately turner . Consequently, the error introduced in the formulation of our model must be calculated to determine the accuracy of any future potential designs in a real-world experimental set-up. Our active–passive systems must, therefore, be analyzed for the case of a high-permeability cylinder that is not a closed perfect magnetic conductor. To do this, we use COMSOL Multiphysics^® working in the magnetostatic regime, to determine how the uniformity of the $B_{x}$ field generated by the constant transverse field-generating system in Fig. 3a changes when the perfect closed magnetic conductor is replaced by one that is imperfect. To construct an imperfect shield, we define the magnetic permeability $\mu_{r}=20000$ to match the permeability of industrial standard mumetal regularly used as a passive shield. We then vary the wall thickness, $d$, of a closed magnetic shield and determine how small $d$ can be while maintaining high field uniformity. Finally, using this minimum thickness, we introduce a circular axial entry hole of radius $\rho_{h}$ at the center of each end cap and determine how large the hole can be to preserve field uniformity.

We use the root mean square field deviation, ${\Delta}B_{x}^{\mathrm{RMS}}$, of the attained field from the uniform target field, calculated along the $z$-axis of the optimized field region, to evaluate the performance of the active–passive system. The light blue crosses in Fig. 6a show that ${\Delta}B_{x}^{\mathrm{RMS}}$ values calculated numerically from the current continuum in Fig. 3a using COMSOL Multiphysics^®, decrease as $d$ increases in the range $(0.05-2.5)$ mm with interval size $0.05$ mm, converging to the thick material limit where the thickness is assumed to be infinite. The horizontal dashed red line shows the analytical value of ${\Delta}B_{x}^{\mathrm{RMS}}=0.0232\%$ calculated using (37)-(39). The numerical ${\Delta}B_{x}^{\mathrm{RMS}}$ values decrease asymptotically below this analytical limit and approach the difference $\mathcal{O}(\mu_{r}^{-1})\approx 0.005\%$ mu ; mu1 that we predicted for our model in Section II. This intrinsic error, resulting from the small difference between thick high-permeability materials and a perfect magnetic conductor, sets the hard limit on the accuracy of any magnetic field that can by designed using our theoretical model. In reality, however, this limit is so small that for a thick material with a high permeability, such as mumetal, the errors in manufacturing and construction will be much more significant. As technologies advance which reduce these system errors, such as screen-printed foldable PCBs PCBCoils and 3D-printing technologies 3DPrintingPaper , it may become more relevant to develop a model which accounts ab initio for magnetic shields of finite permeability and thickness.

We see from Fig. 6a that the asymptotic limit is reached at approximately $d=1$ mm, where ${\Delta}B_{x}^{\mathrm{RMS}}=0.019\%$. At this point, regarding the accuracy of our model, there is little advantage to increasing $d$ further. Consequently, in Fig. 6b we take $d=1$ mm and examine the effect of introducing circular axial entry holes in both end caps of the high-permeability cylinder. Although ${\Delta}B_{x}^{\mathrm{RMS}}$ increases as the hole radius, $\rho_{h}$, increases from no hole, $\rho_{h}=0$, to when there is no end cap, $\rho_{h}=\rho_{s}$, we see that small holes in both end caps allow experimental access without significantly reducing the fidelity of the desired field profile. In particular, for our system, the hole radius can be made as large as $\rho_{h}=0.25\rho_{s}$ while only increasing ${\Delta}B_{x}^{\mathrm{RMS}}$ by $0.0012\%$ when compared to the no hole case (horizontal dashed light blue line).

In Fig. 7 we show the numerically-calculated color maps of the $B_{x}$ field in the $y$-$z$ plane generated by the constant transverse field (Fig. 3a) and linear transverse field gradient (Fig. 4a) systems. Both of these systems are simulated with the same imperfect high-permeability cylindrical magnetic shield that has had its properties determined in the above analysis: $\mu_{r}=20000$, $d=1$ mm, and $\rho_{h}=0.25\rho_{s}$. The performance of both systems in terms of the cylindrical shield fraction is summarized in Table 1.

Finally, we see from Table 1 that this imperfect magnetic shield does not introduce significant magnetic field deviations above $0.1\%$ when compared to a perfect shield with the same geometry. In particular, the maximum difference between the perfect and imperfect cylindrical shield fractions for deviations above $0.1\%$ is only $0.028$. Large deviations in the field accuracy below $0.1\%$ can be graphically seen when comparing the color maps for the perfect (Fig. 5) and imperfect (Fig. 7) cases. The contours showing field deviations of $0.01\%$ and $0.001\%$ are strongly perturbed, as expected from the analysis in Fig. 6, demonstrating the hard intrinsic limit on our model when generating target field profiles inside real-world magnetic shields. Similar analysis should be applied when designing other active–passive systems using our theoretical model in order to quantify its accuracy for a specific experimental setup. Further analysis could also be performed to determine the low-frequency limit in which a time-dependent current source could be included. However, if the magnitudes of any induced eddy currents are much less than the magnitude of the coil current, such effects will be negligible. Further designs generating more exotic field profiles can be found in the Supplementary Material. All designs can be found in our open access Python code which can be used to design systems to generate specific physical target field using our model (see addendum for more details).

In this paper, we have developed an analytical model to calculate the total magnetic field generated by an arbitrary current flow on a cylinder that is coaxially nested within a finite closed high-permeability cylinder. We modified the Green’s function for the magnetic vector potential, matched the radial and planar boundary conditions of the magnetic field through the introduction of a pseudo-current density, and incorporated a harmonic minimization procedure to design optimal user-specified magnetic fields by using a modified cylindrical Fourier basis. We then verified this optimization procedure by designing coils to generate a constant transverse field, $B_{x}$, across the cylinder, and a linear transverse field gradient, $\mathrm{d}B_{x}/\mathrm{d}z$, along the length of the cylinder. Our analytical calculations of these field profiles agreed well with numerical simulations. The optimization procedure generated highly accurate $B_{x}$ and $\mathrm{d}B_{x}/\mathrm{d}z$ field profiles, inside a high-permeability cylinder, with peak-to-peak deviations from the target profiles below $0.11\%$ and $0.24\%$, respectively. The analytically-predicted deviations agreed with the numerical simulations to within $0.002\%$ and $0.003\%$ for the constant and linear gradient systems, respectively.

We further investigated the validity of our model by analyzing the behavior of the constant transverse field-generating system inside a magnetic shield of permeability $\mu_{r}=20000$, finite thickness, and with circular axial entry holes in the end caps. We found a range of parameters where the analytical predictions for a perfect cylindrical magnetic conductor remain close to numerical simulations for a cylindrical shield with finite permeability and thickness that includes entry holes; showing that the designs generated by our model are applicable to real-world magnetic shields. Notably, when the active field-generating systems were enclosed by a passive magnetic shield with realistic experimental parameters ($\mu_{r}=20000$, thickness $1$ mm, and entry holes of radius equal to $25$% of the shield’s radius), the deviation from the desired constant and linear gradient field profiles were less than $0.1\%$ and $0.5\%$, respectively, over more than $40\%$ of the central radial and axial extent of this simulated real-world magnetic shield.

Our flexible optimization procedure enables the design of new active–passive magnetic field shaping systems to generate accurately any static magnetic field profile in the interior of a finite closed magnetic shield, subject to satisfying Maxwell’s equations in free space and not saturating the shielding material. This facilitates the development and miniaturization of systems and technologies which require such control, including quantum sensors, fundamental physics experiments, and medical technologies. Further investigation could consider an analytical treatment of finite magnetic shield thickness and permeability and interactions with an open magnetic shield topology.

Acknowledgements  We acknowledge support from the UK Quantum Technology Hub for Sensing and Timing, funded by the Engineering and Physical Sciences Research Council (EP/M013294/1).
Competing Interests  The authors M. Packer, P. J. Hobson, T. M. Fromhold, M. J. Brookes, and R. Bowtell declare that they have a patent pending to the UK Government Intellectual Property Application office (Application Number 1913549.0) regarding the magnetic field optimization techniques described in this work. The authors have no other competing financial interests.
Correspondence  Correspondence and requests should be addressed to Mark.Fromhold@nottingham.ac.uk.
Open Access  The Python code used to design arbitrary cylindrical coils in a magnetically shielded cylinder can be found in the following GitHub repository: https://github.com/peterjhobson/fields_in_shields. Verification using COMSOL Multiphysics^® requires a valid license.