Optimal Inverse Design of Magnetic Field Proﬁles in a Magnetically Shielded Cylinder

Magnetic shields that use both active and passive components to enable the generation of a tailored low-ﬁeld environment are required for many applications in science, engineering, and medical imaging. Until now, accurate ﬁeld nulling, or ﬁeld generation, has only been possible over a small fraction of the overall volume of the shield. This is due to the interaction between the active ﬁeld-generating components and the surrounding high-permeability passive shielding material. In this paper, we formulate the interaction between an arbitrary static current ﬂow on a cylinder and an exterior closed high-permeability cylinder. We modify the Green’s function for the magnetic vector potential and match boundary conditions on the shield’s interior surface to calculate the total magnetic ﬁeld generated by the system. We cast this formulation into an inverse optimization problem to design active–passive magnetic ﬁeld shaping systems that accurately generate any physical static magnetic ﬁeld in the interior of a closed cylindrical passive shield. We illustrate this method by designing hybrid systems that generate a range of magnetic ﬁeld proﬁles to high accuracy over large interior volumes, and simulate them in real-world shields whose passive components have ﬁnite permeability, thickness, and axial entry holes. Our optimization procedure can be adapted to design active–passive magnetic ﬁeld shaping systems that accurately generate any physical user-speciﬁed static magnetic ﬁeld in the interior of a closed cylindrical shield of any length, enabling the development and miniaturization of systems that require accurate magnetic shielding and control. DOI:


I. INTRODUCTION
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 [1][2][3][4]; atomic magnetometry [5][6][7] for applications including material characterization [8] and magnetoencephalography [9,10]; noise reduction in fundamental physics experiments [11,12] such as timing using high-precision atomic clocks [13,14], measuring the electric dipole moment of fundamental systems [15][16][17][18], and testing Lorentz invariance by observing the limits of spin precession [19][20][21]. 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 [22][23][24]. * mark.fromhold@nottingham.ac.uk 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 [25]. 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 [26][27][28] 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. Pseudoanalytical formulations using a Green's function expansion and harmonic minimization methods have been developed to generate accurate userspecified magnetic fields in free space [29][30][31] to facilitate the need for fast and accurate design procedures. However, these methods did not incorporate the interaction of the high-permeability passive shielding material, hindering the accurate generation of specified target magnetic field profiles [32] in shielded environments. Consequently, optimization of magnetic field cancelation, 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 [33,34] to evolve the coil parameters iteratively, including the effect of the passive material on the magnetic field generated actively at each iteration. Because of 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 nonlinearly 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 [35]. 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 [36][37][38]. 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 [39] has been used to formulate the response of a planar material with high magnetic permeability to a magnetic field generated by a current source [40][41][42][43]. 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 [44], 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 previous work [45], we first derive the 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 that 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) [46,47]. 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, μ r 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, finding globally optimal solutions required for practical, low size, weight, power, and cost technologies for the applications listed above.

II. THEORY
The interface conditions for a magnetic field along a boundary S between two materials are that the normal component of the magnetic field 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 and wheren is the unit vector normal to the boundary and μ 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 Physically, the magnetization of the subdomains of a material change in response to an applied magnetic field to satisfy boundary condition (2) at the material's surface.
Here, we choose to formulate this response in terms of a pseudocurrent density J confined to the surface of the material, which relates to the curl of the magnetization, The magnetic field strength resulting from a current source J c can be expressed using Ampère's law Substituting Eqs. (4) and (5) As shown in many papers and textbooks [39], the solution to the Poisson equation, for an arbitrary current density J is given by where G(r, r ) is the associated Green's function.
Let us now use the boundary condition on the magnetic field, Eq. (2), to examine the effect of a high-permeability material on the magnetic field generated by a specific current distribution. Consider a hollow high-permeability, μ r → ∞, cylinder of radius ρ s , length L s , and thickness d with high-permeability planar end caps located at z = ±L s /2 that is surrounded by free space (Fig. 1). An arbitrary static current flows over an internal open coaxial cylinder of radius ρ c < ρ s and length L 1 − L 2 = L c < L s , where −L s /2 < L 2 < L 1 < L s /2, as shown in Fig. 1. Since high-permeability materials, such as mumetal, can have μ r values greater than 100 000 times that of air, the magnetic field must abruptly change direction at the boundary between its surface and air to satisfy boundary conditions (1)- (2). When the shield is of sufficient thickness and permeability, the tangential components of the magnetic field at its boundary must be approximately zero. Boundary condition (2) on the interior surface of the exterior closed cylinder for the example depicted in Fig. 1 requires that where the magnetic field B = B ρρ + B φφ + B 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 O(μ −1 r ) [48,49]. Therefore, we may assume that, for high-permeability materials, such as mumetal, the response can be approximated to that of a perfect magnetic conductor without introducing significant errors.
Because of the symmetries of the system, it is natural to work in cylindrical coordinates. Following the formulation of Turner [45], we express vector potential (7) due to an arbitrary current distribution Since the current has been restricted to flow over an interior cylindrical shell centered radially about the origin, its radial components may be set to zero, resulting in the simplified vector potentials Substituting the Green's function solution from Eq. (7) in cylindrical coordinates [39], 054004-3 into Eqs. (11)- (13), the components of the vector potential may be expressed as where ρ is the radius of the cylinder and ρ > , ρ < is either ρ, ρ if ρ > ρ or ρ , ρ if ρ < ρ . Equations (15)- (17) contain Fourier transforms of the current densities, with their corresponding inverse transforms given by 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 [50], 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. Because of 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 pseudocurrent density induced on an infinite cylinder. As the planar end caps are spatially orthogonal to the pseudocurrent 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 pseudocurrent density induced on the high-permeability cylindrical shell and, hence, satisfy the boundary conditions over the entire domain of the closed finite highpermeability cylinder. As a result, using Eqs. (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 ρ c ≤ ρ ≤ ρ s , including the effect of the high-permeability cylinder, may be written in terms of an infinite summation over the reflected image currents, where are the Fourier transforms of the pth reflected image current and J mp φ,z (k) is the Fourier transform of the pth reflected pseudocurrent induced on the high-permeability cylinder. In Fig. 2 we depict how azimuthal, J p φ (φ , z ), and axial, J p z (φ , z ), image currents are generated by two parallel planar perfect magnetic conductors.
Writing the magnetic field in cylindrical coordinates, boundary conditions (8) on the magnetic field at ρ = ρ s are Substituting Eqs. (22)-(24) into Eqs. (28)- (29), the Fourier transformed pseudocurrent density on the cylindrical wall of the high-permeability cylinder is found to be where I m (z) and K m (z) are the derivatives with respect to z of I m (z) and K m (z), respectively. Using the continuity equation, Eqs. (20)- (24), and Eq. (30), the magnetic field interior to the conducting cylinder, for all points ρ < ρ c , is given by

054004-5
where 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 this end, we use a modified Fourier basis, defining the pth reflected azimuthal current to be in which where the functions 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, Eqs. (37)- (39), to determine the unknown Fourier coefficients (W n0 , W nm , Q nm ) for a specified target magnetic field. Following work done by Carlson et al. [46], 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 is given by which, when integrated over the surface of the cylinder, using the continuity equation and Eq. (34), gives We now construct a cost function using a least-squares optimization procedure, evaluated at K target field points, where β 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, 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 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 higherorder 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 dis- 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 ϕ is adjusted.

III. RESULTS
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 [36][37][38] and shown in Appendix B, 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 Figs We quantify this agreement by analyzing the deviation from the target fields in the optimization region, B x and dB x /dz, 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. 3 ; numerically using COMSOL Multiphysics version 5.3a and modeling the high-permeability cylinder as a perfect magnetic conductor (blue dotted curve); numerically without the high-permeability cylinder and using the Biot-Savart law with N ϕ = 100 contour levels (dashed green). (c) Enlarged section of (b) emphasizing the high level of field uniformity and the agreement between the numerical and analytical results over the optimization region.
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. 3(a)] and linear transverse field gradient [ Fig. 4(a)] systems, using COMSOL Multiphysics with the shield treated as a perfect magnetic conductor. We summarize the performance of both systems in Table I, calculating the cylindrical shield fractions-the ratio of the radial and axial extents 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%.
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  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 [51], and it is possible to manufacture structures that approximate the current continuum accurately [45]. 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 setup. 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. 3(a) changes when the perfect closed magnetic conductor is replaced by one that is imperfect. To construct an imperfect shield, we define the magnetic permeability μ r = 20 000 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 ρ 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 (rms) field deviation B rms x 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. 6(a) show that B rms x values calculated numerically from the current continuum in Fig. 3(a) 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 B rms x = 0.0232% calculated using Eqs. (37)- (39). The numerical B rms x values decrease asymptotically below this analytical limit and approach the difference O(μ −1 r ) ≈ 0.005% [48,49] that we predicted for our model in Sec. 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 that reduce these system errors, such as screen-printed foldable printed circuit boards [52] and three-dimensional printing technologies [53], it may become more relevant to develop a model that accounts ab initio for magnetic shields of finite permeability and thickness.
We see from Fig. 6(a) that the asymptotic limit is reached at approximately d = 1 mm, where B rms x = 0.019%. At this point, regarding the accuracy of our model, there is little advantage to increasing d further. Consequently, in Fig. 6(b) 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 B rms x increases as the hole radius ρ h increases from no hole, ρ h = 0, to when there is no end cap, ρ h = ρ 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 ρ h = 0.25ρ s while only increasing B rms x 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. 3(a)] and linear transverse field gradient [ Fig. 4(a)] 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: μ r = 20 000, d = 1 mm, and ρ h = 0.25ρ s . The performance of both systems in terms of the cylindrical shield fraction is summarized in Table I.
Finally, we see from Table I 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 054004-11 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 timedependent 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 Appendix C. All designs can be found in our open access PYTHON code that can be used to design systems to generate a specific physical target field using our model (see below for more details).

IV. CONCLUSION
In this paper, we develop 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 modify Green's function for the magnetic vector potential, match the radial and planar boundary conditions of the magnetic field through the introduction of a pseudocurrent density, and incorporate a harmonic minimization procedure to design optimal user-specified magnetic fields by using a modified cylindrical Fourier basis. We then verify this optimization procedure by designing coils to generate a constant transverse field B x across the cylinder, and a linear transverse field gradient dB x /dz along the length of the cylinder. Our analytical calculations of these field profiles agree well with numerical simulations. The optimization procedure generates highly accurate B x and dB x /dz field profiles, inside a high-permeability cylinder, with peak-topeak deviations from the target profiles below 0.11% and 0.24%, respectively. The analytically predicted deviations agree with the numerical simulations to within 0.002% and 0.003% for the constant and linear gradient systems, respectively.
We further investigate the validity of our model by analyzing the behavior of the constant transverse fieldgenerating system inside a magnetic shield of permeability μ r = 20 000, finite thickness, and with circular axial entry holes in the end caps. We find 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 are enclosed by a passive magnetic shield with realistic experimental parameters (μ r = 20 000, 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 are less than 0.1% and 0.5%, respectively, over more than 40% of the central radial and axial extents of this simulated real-world magnetic shield.
Our flexible optimization procedure enables the design of new active-passive magnetic field shaping systems to accurately generate 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 that 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.
The PYTHON code used to design arbitrary cylindrical coils in a magnetically shielded cylinder is openly available from GitHub [54]. Verification using COMSOL Multiphysics requires a valid license.

APPENDIX B: SOLENOIDAL COIL
As shown in previous work [36][37][38], a solenoid of the same length as the high-permeability cylinder provides the most optimal solution for generating a constant axial field. Because of the image currents, the finite solenoid effectively acts as one of infinite extension, resulting in the most homogeneously possible magnetic field in the z direction. A perfect finite solenoid of length L c carries a total current of where I is the current density, i.e., the current per unit solenoid length, resulting in the azimuthal current I t /L c . Using Eq. (25), the Fourier transform of the current through a finite solenoid is given by This summation can be simplified as the only contributing term is p = 0, which, when evaluated, results in This result might seem counterintuitive because the magnetic field is identical to a long solenoid in free space with N turns, i.e., B z (ρ, φ, z) = μ 0 IN /L c , with no field created by the cylindrical surface of the passive shield. This is, however, entirely physical as an infinite solenoid generates zero magnetic field outside of the solenoid itself.
Consequently, there is no response due to the surface of the cylindrical wall and, therefore, no magnetic field generated by the cylindrical surface of the perfect magnetic conductor, i.e., the high-permeability cylinder.

APPENDIX C: EXAMPLE COIL DESIGNS
In Figs. 8 and 9, we present further hybrid active-passive systems, which generate more complex magnetic field landscapes beyond the constant field and field gradient considered in the main body of the paper. The coordinate axes and magnetic field plots are labeled in the same way as the systems presented in the main text. The design in Fig. 8 generates uniform axial magnetic fields B z with different strengths in two different regions of a magnetic shield. The design in Fig. 9 generates a quadratic gradient field whose spatial variation matches the harmonic B = [2xzx − 2yzŷ + (x 2 − y 2 )ẑ]. Further systems can be designed using the PYTHON code in the GitHub repository [54].