Planar Coil Optimization in a Magnetically Shielded Cylinder
M. Packer, P.J. Hobson, N. Holmes, J. Leggett, P. Glover, M.J. Brookes, R. Bowtell, T.M. Fromhold
PPlanar Coil Optimization in a Magnetically Shielded Cylinder
M. Packer , P. J. Hobson , N. Holmes , J. Leggett ,P. Glover , M. J. Brookes , R. Bowtell , and T. M. Fromhold School of Physics and Astronomy, University of Nottingham, Nottingham, NG7 2RD, UK. Sir Peter Mansfield Imaging Centre, School of Physics and Astronomy,University of Nottingham, Nottingham, NG7 2RD, UK. (Dated: January 19, 2021)Hybrid magnetic shields with both active field generating components and high-permeabilitymagnetic shielding are increasingly needed for various technologies and experiments that requireprecision-controlled magnetic field environments. However, the fields generated by the active compo-nents interact with the passive magnetic shield, distorting the desired field profiles. Consequently,optimization of the active components needed to generate user-specified target fields must includecoupling to the high-permeability passive components. Here, we consider the optimization of planaractive systems, on which an arbitrary static current flows, coupled to a closed high-permeability cylin-drical shield. We modify the Green’s function for the magnetic vector potential to match boundaryconditions on the shield’s interior surface, enabling us to construct an inverse optimization problemto design planar coils that generate user-specified magnetic fields inside high-permeability shields. Wevalidate our methodology by designing two bi-planar hybrid active–passive systems, which generatea constant transverse field, B = ˆx , and a linear field gradient, B = ( − x ˆx − y ˆy + 2 z ˆz ) , respectively.For both systems, the inverse-optimized magnetic field profiles agree well with forward numericalsimulations. Our design methodology is accurate and flexible, facilitating the miniaturization ofhigh-performance hybrid magnetic field generating technologies with strict design constraints andspatial limitations. I. INTRODUCTION
Tailored, high-precision, low magnetic field environments are required for many applications, devices, and experiments.Examples include magnetic field control in quantum sensing of gravity for underground surveying and mapping [1–4],magnetic field cancellation for atomic magnetometry [5] with applications in medical imaging such as magnetoen-cephalography [6–9] and neonatal/fetal magnetocardiography [10–12], and noise suppression in fundamental physicsexperiments [13–15]. Usually, these systems are enclosed by a shield formed from high-permeability material, whichexcludes stray external magnetic fields that can limit the accuracy and sensitivity of the measurements. In particular,cylindrical shield geometries are often used because the dimensions and spacings of multiple cylindrical shield layerscan be optimized to generate a large interior shielded region [16, 17]. The magnetic field interior to the shield can thenbe adjusted by active field generating components to either cancel background fields further or to define a specific fieldenvironment. However, the surrounding passive shielding material deforms the magnetic field profiles generated by theactive components, making it hard to design wire patterns which accurately generate specified target magnetic fieldprofiles [18].Boundary element methods (BEMs) can be used to optimize magnetic fields generated by surface currents on atriangular mesh [19–22] to generate arbitrary target magnetic fields. BEMs are extremely powerful and flexible sincethey can be used to define active systems with complex geometries inside passive shields. They are, however, limitedby computational power with results depending on mesh size and on the distance of the active components from theshielding material. Alternatively, analytical methods for optimizing hybrid active–passive systems are advantageoussince they provide intuition and understanding of how magnetic fields are distorted by the presence of high-permeabilitymaterial. Analytical formulations also provide a fast and efficient route to determine the best system for generating abespoke user-specified magnetic field profile.Currently, however, analytical models for hybrid active–passive systems are restricted to a limited number of scenarios.Simple discrete coil geometries have been formulated in cylindrical high-permeability magnetic shields, where themagnetic field is decomposed into azimuthal Fourier modes [23–25] and matched at the shield boundary. Planar high-permeability materials have been incorporated into optimization procedures using the method of mirror images [26, 27]to determine the total magnetic field generated by a static current source inside a magnetically shielded room [28–31].Similar work in magnetic resonance imaging (MRI) has investigated the interaction of switched magnetic field gradientswith high-conductivity materials used for passive shielding [32]. More recently, solutions for the interaction between astatic current source and a finite closed cylindrical high-permeability shield have been formulated for the specific case ofcylindrical co-axial surface current geometries [33]. These formulations incorporate the effect of the shielding materialby explicitly solving Maxwell’s equations by matching the required boundary conditions on the surface interface tofind the total magnetic field generated by the system. a r X i v : . [ phy s i c s . a pp - ph ] J a n Building on this previous work, here we determine the effect of a closed finite length high-permeability cylinder on themagnetic field generated by an arbitrary static current distribution on an interior circular plane, oriented perpendicularto the axis of the cylinder. We calculate the effect of the high-permeability cylinder on the magnetic field produced bycurrent flow on the plane, which enables us to determine optimal current paths to generate user-specified target fieldsinside the cylinder. We first derive the vector potential generated by an arbitrary planar current source. Then, wecalculate a pseudo-current density induced on the cylindrical surface of the high-permeability material in responseto the planar current source, and, hence, derive a Green’s function for our system. Finally, we implement a Fourierdecomposition of the current paths to calculate the total magnetic field in terms of a set of weighted Fourier coefficients.This formulation allows the incorporation of a quadratic optimization procedure to determine globally optimal designsthat generate user-specified target magnetic field profiles. Through the use of this optimization procedure, we designtwo example bi-planar coil systems optimized for operation inside a finite-length closed high-permeability cylindricalmagnetic shield. In both cases, we confirm that our analytical model agrees well with the result of numerical finiteelement simulations. Our work extends the range of coil geometries which can be efficiently optimized to generatestatic user-specified target magnetic fields in the presence of high-permeability materials, expanding design flexibilityfor systems which require precision-controlled magnetic field environments.
II. THEORY
High-permeability magnetic shields, usually constructed from materials such as mumetal, are used to attenuateexternal background magnetic fields in order to shield magnetically sensitive equipment from spurious signals. Thesehigh-permeability materials generate an induced magnetization, M , in response to a incident field at their surface, S . This ensures continuity of the parallel and tangential components of the magnetic field, B , and magnetic fieldstrength, H , respectively, at the interface between air and a material. Considering this interface while working in themagnetostatic regime with no surface currents, the boundary conditions take the form ( B mat . − B air ) · ˆn = 0 on S, (1)and (cid:18) µ r B mat . − B air (cid:19) ∧ ˆn = 0 on S, (2)where ˆ n is the unit vector normal to the boundary and µ r is the relative permeability of the material. To designmagnetic fields effectively in shielded environments, we need to determine the total field generated by an active coilstructure and the high-permeability material. The total magnetic field in free space is related to the magnetic fieldstrength and the magnetization by B = µ ( H + M ) . (3)The induced magnetization can be formulated in terms of a pseudo-current density ˜j , confined to the material’s surface ∇ ∧ M = (cid:101) j , (4)with the magnetic field strength related to a current density j c , on an active structure, using Ampère’s law, ∇ ∧ H = j c . (5)By using (3)-(5), and the relation between the vector potential and the magnetic field, B = ∇ ∧ A , the total vectorpotential in free space generated by the system can be cast as the Poisson equation, ∇ A = − µ ( j c + (cid:101) j ) , (6)resulting in the integral solution in terms of an arbitrary current density A ( r ) = µ (cid:90) r (cid:48) d r (cid:48) G ( r , r (cid:48) ) j ( r (cid:48) ) , (7)where G ( r , r (cid:48) ) is the associated Green’s function for the system [26].Here, we consider the specific example of a closed cylindrical magnetic shield surrounded by free space, as shown inFig. 1. This cylinder has relative permeability, µ r (cid:29) , radius ρ s , and length L s , with planar end caps located at z = ± L s / . Inside this cylinder, an arbitrary current flows on a circular planar surface of radius ρ c < ρ s , centered onthe z -axis and lying in the z = z (cid:48) plane where | z (cid:48) | < L s / . Previous work has shown that the magnetic field generated xyz ρ c z = z (cid:48) ρ s z = L s / z = − L s / FIG. 1: A cylindrical hollow high-permeability magnetic shield of length L s , radius ρ s , and with planar end capslocated at z = ± L s / encloses an interior circular current source (bounded by the upper dashed curve) of radius ρ c < ρ s , which lies in the z = z (cid:48) plane.by a high-permeability material in response to an applied static field deviates from that of a perfect magnetic conductoron the scale of O (cid:0) µ − r (cid:1) [34, 35]. More specifically, in the case of a closed finite length cylinder with permeability µ r = 20000 , equal to industrial standard mumetal, a material thickness of d = 1 mm is sufficient to provide a responsesimilar to that of bulk material, differing from the perfect magnetic conducting response only by approximately . [33]. Therefore, we can approximate sufficiently thick materials with high-permeability µ r > , such as mumetal,to that of a perfect magnetic conductor without compromising the accuracy of a theoretical model. Assuming thesurface of the shield shown in Fig. 1 is a perfect magnetic conductor, the boundary condition at its surface can bewritten as B ρ (cid:12)(cid:12)(cid:12)(cid:12) z = ± L s / = 0 , B φ (cid:12)(cid:12)(cid:12)(cid:12) z = ± L s / ,ρ = ρ s = 0 , B z (cid:12)(cid:12)(cid:12)(cid:12) ρ = ρ s = 0 , (8)where the magnetic field B = B ρ ˆ ρ + B φ ˆ φ + B z ˆz , is expressed in cylindrical polar coordinates.As previously formulated for cylindrical current flow [33], the response generated by a finite closed cylinder withhigh-permeability, µ r (cid:29) , can be determined by matching the orthogonal modes in the magnetic field and generatinga relation between the initial current flow and the resulting magnetization induced on the surface of the cylinder. Herewe apply the same methodology to generate an analytical expression for the response of a finite closed high-permeabilitycylinder to a planar current source using a similar decomposition. The vector potential generated by a current sourcein cylindrical coordinates can be expressed as, A ρ ( r ) = µ (cid:90) r (cid:48) d r (cid:48) G ( r , r (cid:48) ) [ j ρ ( r (cid:48) ) cos ( φ − φ (cid:48) ) + j φ ( r (cid:48) ) sin ( φ − φ (cid:48) )] , (9) A φ ( r ) = − µ (cid:90) r (cid:48) d r (cid:48) G ( r , r (cid:48) ) [ j ρ ( r (cid:48) ) sin ( φ − φ (cid:48) ) − j φ ( r (cid:48) ) cos ( φ − φ (cid:48) )] , (10) A z ( r ) = µ (cid:90) r (cid:48) d r (cid:48) G ( r , r (cid:48) ) j z ( r (cid:48) ) . (11)Since there is no current flow in the z -direction for a planar current source perpendicular to the axis of the cylindricalshield, the continuity equation can be used to express the planar current flow in terms of a single scalar streamfunctiondefined on the planar surface j ρ ( r (cid:48) ) = 1 ρ (cid:48) ∂ϕ ( r (cid:48) ) ∂φ (cid:48) , j φ ( r (cid:48) ) = − ∂ϕ ( r (cid:48) ) ∂ρ (cid:48) , (12)where ϕ ( r (cid:48) ) = ϕ ( ρ (cid:48) , φ (cid:48) ) . To exploit the radial symmetries of the system, we choose to decompose the Green’s functionin cylindrical coordinates in terms of Bessel functions of the first kind, G ( r , r (cid:48) ) = 14 π ∞ (cid:88) m = −∞ e im ( φ − φ (cid:48) ) (cid:90) ∞ d k e − k | z − z (cid:48) | J m ( kρ ) J m ( kρ (cid:48) ) , (13)allowing the vector potential to be expressed in terms of cylindrical harmonics defined on a circular plane. Using(9)-(11), (12), and (13) we cast the vector potential in a simplified form A ρ ( r ) = iµ ρ ∞ (cid:88) m = −∞ e imφ (cid:90) ∞ d k me − k | z − z (cid:48) | J m ( kρ ) ϕ m ( k ) , (14) A φ ( r ) = − µ ∞ (cid:88) m = −∞ e imφ (cid:90) ∞ d k ke − k | z − z (cid:48) | J (cid:48) m ( kρ ) ϕ m ( k ) , (15) A z ( r ) = 0 , (16)where J (cid:48) m ( z ) is the derivative of J m ( z ) with respect to z , and ϕ m ( k ) is defined as the m th order Hankel transform ϕ m ( k ) = 12 π (cid:90) ∞ dρ (cid:48) (cid:90) π dφ (cid:48) e − imφ (cid:48) ρ (cid:48) J m ( kρ (cid:48) ) ϕ ( ρ (cid:48) , φ ) . (17)We now consider the vector potential generated by the magnetic shield. To do this, we first introduce a pseudo-currentdensity induced on an infinite cylinder, ˜ j = ˜ j φ ( φ (cid:48) , z (cid:48) ) ˆ φ + ˜ j z ( φ (cid:48) , z (cid:48) ) ˆz . Next, we seek a Fourier representation of thispseudo-current density, which satisfies the boundary condition over the entire domain of the shield. In particular, wewish to equate the shared azimuthal Fourier modes at the radial boundary of the shield cylinder. This is achievedusing a combination of methods that must be applied sequentially, since each method satisfies the condition at anorthogonal boundary. The radial condition is satisfied by equating the magnetic field generated by the cylindricalpseudo-current density and planar current flow, generating a relation between the response of an infinite cylindricalshield and the initial current source. Then, the boundary condition at the end caps can simultaneously be satisfiedby applying the method of mirror images [26]. These methods can be combined because the infinite pseudo-currentdensity and planar current flow are spatially orthogonal to the end caps, meaning that any reflections generated by theapplication of the method of mirror images continue to satisfy the radial condition. The components of the vectorpotential generated by a pseudo-current density induced on an infinite cylinder [32] in the region ρ < ρ s are given by A ρ ( ρ, φ, z ) = − iµ ρ s π ∞ (cid:88) m = −∞ (cid:90) ∞−∞ d k e imφ e ikz (cid:104) I m − ( | k | ρ ) K m − ( | k | ρ s ) − I m +1 ( | k | ρ ) K m +1 ( | k | ρ s ) (cid:105) ˜ j φm ( k ) , (18) A φ ( ρ, φ, z ) = µ ρ s π ∞ (cid:88) m = −∞ (cid:90) ∞−∞ d k e imφ e ikz (cid:104) I m − ( | k | ρ ) K m − ( | k | ρ s )+ I m +1 ( | k | ρ ) K m +1 ( | k | ρ s ) (cid:105) ˜ j φm ( k ) , (19) A z ( ρ, φ, z ) = µ ρ s π ∞ (cid:88) m = −∞ (cid:90) ∞−∞ d k e imφ e ikz I m ( | k | ρ ) K m ( | k | ρ s ) ˜ j zm ( k ) , (20)where the Fourier transforms of the pseudo-currents are defined by j mφ ( k ) = 12 π (cid:90) π d φ (cid:48) e − imφ (cid:48) (cid:90) ∞−∞ d z (cid:48) e − ikz (cid:48) ˜ j φ ( φ (cid:48) , z (cid:48) ) , (21) j mz ( k ) = 12 π (cid:90) π d φ (cid:48) e − imφ (cid:48) (cid:90) ∞−∞ d z (cid:48) e − ikz (cid:48) ˜ j z ( φ (cid:48) , z (cid:48) ) . (22)The corresponding inverse transforms are given by ˜ j φ ( φ (cid:48) , z (cid:48) ) = 12 π ∞ (cid:88) m = −∞ (cid:90) ∞−∞ d k e imφ (cid:48) e ikz (cid:48) j mφ ( k ) , (23) ˜ j z ( φ (cid:48) , z (cid:48) ) = 12 π ∞ (cid:88) m = −∞ (cid:90) ∞−∞ d k e imφ (cid:48) e ikz (cid:48) j mz ( k ) . (24)Therefore, adding the contributions from the planar current flow, (14)-(15), and the infinite pseudo-current density,(18)-(20), while using (21)-(24), and applying the method of mirror images for two infinitely large parallel planes, wecan write the total magnetic field generated by the system in the region ρ < ρ s as B ρ ( ρ, φ, z ) = − µ π ∞ (cid:88) p = −∞ ∞ (cid:88) m = −∞ e imφ (cid:34) (cid:90) ∞ d k k ( z − ( − p z (cid:48) + pL s ) | z − ( − p z (cid:48) + pL s | e − k | z − ( − p z (cid:48) + pL s | J (cid:48) m ( kρ ) ϕ m ( k ) − iρ s π (cid:90) ∞−∞ d k ke ikz I (cid:48) m ( | k | ρ ) K (cid:48) m ( | k | ρ s ) j mpφ ( k ) (cid:35) , (25) B φ ( ρ, φ, z ) = − iµ ρ ∞ (cid:88) p = −∞ ∞ (cid:88) m = −∞ me imφ (cid:34) (cid:90) ∞ d k k ( z − ( − p z (cid:48) + pL s ) | z − ( − p z (cid:48) + pL s | e − k | z − ( − p z (cid:48) + pL s | J m ( kρ ) ϕ m ( k ) − iρ s π (cid:90) ∞−∞ d k | k | k e ikz I m ( | k | ρ ) K (cid:48) m ( | k | ρ s ) j mpφ ( k ) (cid:35) , (26) B z ( ρ, φ, z ) = µ ∞ (cid:88) p = −∞ ∞ (cid:88) m = −∞ e imφ (cid:34) (cid:90) ∞ d k k e − k | z − ( − p z (cid:48) + pL s | J m ( kρ ) ϕ m ( k ) − ρ s π (cid:90) ∞−∞ d k | k | e ikz I m ( | k | ρ ) K (cid:48) m ( | k | ρ s ) j mpφ ( k ) (cid:35) , (27)where j mpφ ( k ) is the p th reflected Fourier transformed azimuthal pseudo-current density induced on the cylindricalsurface of the magnetic shield with I (cid:48) m ( z ) and K (cid:48) m ( z ) defined as the derivatives with respect to z of the modified Besselfunctions of the first and second kind I m ( z ) and K m ( z ) , respectively. By applying the boundary condition at the radialsurface, (8), we can match the shared m th azimuthal Fourier mode generated by each p th reflected pseudo-current andstreamfunction, resulting in the relation (cid:90) ∞−∞ d k | k | e ikz I m ( | k | ρ s ) K (cid:48) m ( | k | ρ s ) j mpφ ( k ) = πρ s (cid:90) ∞ d k k e − k | z − ( − p z (cid:48) + pL s | J m ( kρ s ) ϕ m ( k ) . (28)Physically, due to the formulation of the response in terms of a pseudo-current density, there must be a uniquesolution that is independent of the axial position that satisfies the boundary condition over the infinite domain ofthe cylindrical shield. Therefore, we perform an inverse Fourier transform with respect to z to generate an integralrepresentation of the p th reflected Fourier pseudo-current density, j mpφ ( k ) , in terms of the m th order Hankel transformof the streamfunction defined on the planar surface, ϕ m ( k ) , j mpφ ( k ) = e − ik ( ( − p z (cid:48) + pL s ) ρ s | k | I m ( | k | ρ s ) K (cid:48) m ( | k | ρ s ) (cid:90) ∞ d˜ k ˜ k ˜ k + k J m (˜ kρ s ) ϕ m (˜ k ) . (29)This expression for j mpφ ( k ) can now be substituted into (25)-(27) to determine the total magnetic field in terms of ϕ m ( k ) . Next, we must choose an appropriate expansion of the streamfunction, ϕ ( ρ (cid:48) , φ (cid:48) ) . Although the choice oforthogonal basis for the expansion of the streamfunction is somewhat arbitrary, a choice of basis which considers thesymmetries between the Hankel transform, coordinate system, and the integral representation of the pseudo-currentyields a simpler solution. Here, we choose to decompose the radial component of the planar current flow into aFourier–Bessel series while using a Fourier series representation of the azimuthal dependence, ϕ ( ρ (cid:48) , φ (cid:48) ) = ( H ( ρ (cid:48) ) − H ( ρ − ρ c )) ρ c N (cid:88) n =1 M (cid:88) m (cid:48) =0 J m (cid:48) (cid:18) ρ nm (cid:48) ρ (cid:48) ρ c (cid:19) ( W nm (cid:48) cos( m (cid:48) φ (cid:48) ) + Q nm (cid:48) sin( m (cid:48) φ (cid:48) )) , (30)where ( W nm (cid:48) , Q nm (cid:48) ) are Fourier coefficients and ρ nm (cid:48) is the n th zero of the m (cid:48) th Bessel function of the first kind, J m (cid:48) ( ρ nm (cid:48) ) = 0 . Therefore, using (25)-(27), (29), and (30), the total magnetic field generated by an arbitrary currentflow on the planar surface inside the closed finite high-permeability cylinder can be written as B ρ ( ρ, φ, z ) = µ ρ c N (cid:88) n =1 M (cid:88) m =0 ρ nm J (cid:48) m ( ρ nm ) ( W nm cos ( mφ ) + Q nm sin ( mφ )) B nmρ ( ρ, z ) , (31) B φ ( ρ, φ, z ) = µ ρ c ρ N (cid:88) n =1 M (cid:88) m =0 mρ nm J (cid:48) m ( ρ nm ) ( W nm sin ( mφ ) − Q nm cos ( mφ )) B nmφ ( ρ, z ) , (32) B z ( ρ, φ, z ) = µ ρ c N (cid:88) n =1 M (cid:88) m =0 ρ nm J (cid:48) m ( ρ nm ) ( W nm cos ( mφ ) + Q nm sin ( mφ )) B nmz ( ρ, z ) , (33)where B nmρ ( ρ, z ) = − (cid:90) ∞ dk k σ ( k ; z, z (cid:48) , L s ) J (cid:48) m ( kρ ) J m ( kρ c ) k ρ c − ρ nm − ∞ (cid:88) p =1 ˜ p | ˜ p | λ p ( z, z (cid:48) , L s ) I (cid:48) m ( | ˜ p | ρ ) I m ( | ˜ p | ρ c ) K m ( | ˜ p | ρ s ) I m ( | ˜ p | ρ s )( | ˜ p | ρ c + ρ nm ) , (34) B nmφ ( ρ, z ) = (cid:90) ∞ d k kσ ( k ; z, z (cid:48) , L s ) J m ( kρ ) J m ( kρ c ) k ρ c − ρ nm + ∞ (cid:88) p =1 ˜ p λ p ( z, z (cid:48) , L s ) I m ( | ˜ p | ρ ) I m ( | ˜ p | ρ c ) K m ( | ˜ p | ρ s ) I m ( | ˜ p | ρ s )( | ˜ p | ρ c + ρ nm ) , (35) B nmz ( ρ, z ) = (cid:90) ∞ d k k γ ( k ; z, z (cid:48) , L s ) J m ( kρ ) J m ( kρ c ) k ρ c − ρ nm − ∞ (cid:88) p =1 ˜ p τ p ( z, z (cid:48) , L s ) I m ( | ˜ p | ρ ) I m ( | ˜ p | ρ c ) K m ( | ˜ p | ρ s ) I m ( | ˜ p | ρ s )( | ˜ p | ρ c + ρ nm ) , (36)and γ ( k ; z, z (cid:48) , L s ) = e − k | z − z (cid:48) | + 2 e kL s − (cid:20) e kL s cosh ( k ( z + z (cid:48) )) + cosh ( k ( z − z (cid:48) )) (cid:21) , (37) σ ( k ; z, z (cid:48) , L s ) = ( z − z (cid:48) ) | z − z (cid:48) | e − k | z − z (cid:48) | − e kL s − (cid:20) e kL s sinh ( k ( z + z (cid:48) )) + sinh ( k ( z − z (cid:48) )) (cid:21) , (38) λ p ( z, z (cid:48) , L s ) = 2 L s (( − p sin (˜ p ( z + z (cid:48) )) + sin (˜ p ( z − z (cid:48) ))) , (39) τ p ( z, z (cid:48) , L s ) = 2 L s (( − p cos (˜ p ( z + z (cid:48) )) + cos (˜ p ( z − z (cid:48) ))) , (40)with ˜ p = pπ/L s . A full derivation of these expressions is given in Appendix A. Solving for the unknown Fouriercoefficients, ( W nm , Q nm ) , to generate a desired magnetic field using the system of governing equations (31)-(33) is anill-conditioned problem due to the formulation of the vector potential through the integral representation in (9)-(11).As in previous work by Carlson et al. [36], this may be solved by a least squares minimization with the addition ofa penalty term that acts as a regularization parameter. The choice of regularization parameter is arbitrary since itexists only to facilitate the inversion. Well-regularized solutions yield more simplistic coil designs at a cost to the fieldfidelity [37]. Here, we choose the regularization parameter to be the power, P , dissipated by a circular planar currentsource of thickness, t , and resistivity, (cid:37) , P = (cid:37)t (cid:90) ρ c d ρ (cid:48) ρ (cid:48) (cid:90) π d φ (cid:48) | J ρ ( ρ (cid:48) , φ (cid:48) ) | + | J φ ( ρ (cid:48) , φ (cid:48) ) | . (41)Substituting (30) into (12) and then (41), and integrating over the planar surface we find that, for m = 0 , P = (cid:37)t πρ c N (cid:88) n =1 W n ρ n J ( ρ n ) , (42)and, for m ∈ Z + , P = (cid:37)t πρ c N (cid:88) n =1 M (cid:88) m =1 (cid:0) W nm + Q nm (cid:1) (cid:16) ρ nm (cid:17) m m !( m − (cid:34) ˜ F (cid:18) m, m + 12 ; m + 1 , m + 1 , m + 1; − ρ nm (cid:19) − ρ nm m + 1) ˜ F (cid:18) m + 12 , m + 1 , m + 1; m, m + 2 , m + 2 , m + 1; − ρ nm (cid:19) (cid:35) , (43)where i ˜ F j is the regularized hypergeometric function: see Appendix B for a full derivation. We can now formulate acost function for the least squares minimization, Φ = K (cid:88) k (cid:2) B desired ( r k ) − B ( r k ) (cid:3) + βP, (44)where β is a weighting parameter chosen to adjust the physical constraints of the system. The cost function isminimized using a least squares fitting to calculate the optimal Fourier coefficients to generate the desired magneticfield at K target points. The minimization is achieved by taking the differential of the cost function with respect tothe Fourier coefficients, ∂ Φ ∂W ij = 0 , ∂ Φ ∂Q ij = 0 , i ≥ , j ≥ , (45)which enables the optimal Fourier coefficients to be found for any given physical target magnetic field profile throughmatrix inversion. The inversion process yields the optimal continuous streamfunction defined on the planar surfacefor a finite number of Fourier coefficients. In the ideal case, the number of Fourier coefficients would be infinite.However, a finite number of terms can provide accurate solutions in well-regularized problems. This number can beapproximated by ensuring that the distance between the planar current-carrying surface and the closest target fieldpoint is much larger than the smallest spatial frequency. This is due to the decreased contribution of higher-orderterms at distances much larger than their spatial frequency. The final objective is to design a coil that generates thedesired magnetic field to a specified accuracy. To do this, a discrete approximation of the field profile may be found bycontouring the streamfunction at N ϕ discrete levels. The contours of the streamfunction generate streamlines wherewires should be laid to replicate the desired target magnetic field. This is achieved by discretizing the streamfunctioninto N ϕ contours, where ϕ j = min ϕ + ( j − / ϕ, j = 1 , ..., N ϕ , separated by ∆ ϕ = max ϕ − min ϕN ϕ , and the totalcurrent through each wire is I = ∆ ϕ . The number of contours should be maximized, limited only by manufacturingsince the accuracy of the theoretical model depends on the quality of the discrete approximation of the streamfunction.The approximation to the streamfunction can be determined by using the elemental Biot–Savart law to calculate theerror as N ϕ is increased. In the case that multiple current-carrying planes are designed, the contours should be definedevenly between the global maximum and minimum of the streamfunction across all the planes so that the currentthrough each wire is equal. III. RESULTS
We now verify our analytical model by designing hybrid active–passive systems composed of bi-planar coils inside aclosed high-permeability magnetic shield and compare the resulting magnetic field profiles with forward numericalsimulations of each optimized system. We consider two distinct systems, each containing current confined to twodisks of radius ρ c = 0 . m and symmetrically placed at z (cid:48) = ± . m. Both systems are interior to a perfect closedcylindrical magnetic shield of radius ρ s = 0 . m and length L s = 1 m centered on the origin, as shown in Fig. 2. Thefirst system is designed to generate a constant transverse field, B = B ˆx , and the second system creates a linear fieldgradient, B = G ( − x ˆx − y ˆy + 2 z ˆz ) , within an optimization region defined by − z (cid:48) / ≤ z ≤ z (cid:48) / and ≤ ρ ≤ ρ c / .The magnetic field profiles chosen, i.e. the uniform transverse and linear field gradient fields, are examples of tesseral( m (cid:54) = 0 ) and zonal ( m = 0 ) harmonics, respectively, which exhibit m -fold and complete azimuthal symmetry [38],which facilitates analysis of the shield’s particular response to tesseral and zonal harmonic fields generated by planarcurrent sources. The two systems that we consider here are chosen to illustrate targeted magnetic field compensation insituations where compact systems are required, but space inside the central cylindrical cavity of the system is limitedby the presence of experimental equipment (e.g. magnetic sensors). Coil designs generating other field profiles, forexample using combinations of planar coils with varying sizes, can be found in the Appendix C. All designs can bereplicated using our open-access Python code[42].In Fig. 3a and Fig. 5a, we show the optimized contoured streamfunctions for the constant transverse and linear fieldgradient systems, respectively. Due to the symmetric placement of the bi-planar coil systems and optimized fieldregions within the shield (see Fig. 2), the magnitudes of the streamfunctions defined on both coils are identical. Thecurrent flow directions, however, are opposite, due to the form of the desired fields. Figures. 3b–d show the transversemagnetic field, B = 1 µ T, along the x -axis and z -axis, respectively, generated by the bi-planar coil design shown inFig. 3a calculated in three different ways: analytically using (31)-(33) (solid red curves); numerically using COMSOLMultiphysics ® with the shield treated as a perfect magnetic conductor (blue dotted curves); numerically in free space,i.e. excluding the high-permeability material and evaluating the magnetic fields through the Biot–Savart law fordiscretized bi-planar coils with N ϕ = 250 (dashed green curve in b). Furthermore, color maps of the transverse and axialmagnetic field components in the xz -plane are presented in Fig. 4a and Fig. 4b, respectively, calculated numericallyusing COMSOL Multiphysics ® with the shield treated as a perfect magnetic conductor. Similarly, Fig. 5b–c show thelinear axial field gradient, G = 1 µ T/m, along the x -axis and z -axis, respectively, generated by the bi-planar coil systemshown in Fig. 5a in the same three cases with similar color maps of the magnetic field components in the xz -planeshown in Fig. 6a and Fig. 6b. For both designs, the analytical field profiles agree well with numerical simulations. Themaximum errors between the analytical model and numerical simulation are . % and . %, for the constanttransverse field and the gradient of the linear gradient field, respectively, along the z -axis of the optimization region.To quantify the performance of our optimization procedure, we can analyze the deviations between the magnetic fieldsgenerated by our theoretical model and the desired target fields. Examining the fidelity of the fields generated byboth systems along x -axis of the optimization region, the maximum absolute deviations from the constant transverseand linear axial gradient fields are . and . , respectively. Along the z -axis of the optimization region, themaximum deviations are . and . , respectively. Clearly, the axial field gradient field is generated moreaccurately than the uniform transverse field. This can be seen from Fig. 3d, which reveals small oscillations in thetransverse field over the z -axis of the optimization region. To understand the difference between the fidelity of the fieldprofiles generated by the two systems, we must analyze how the passive magnetic shield affects the fields from theplanar current distributions, decomposing its response into zonal ( m = 0 ) and tesseral ( m (cid:54) = 0 ) harmonic components.These harmonic responses relate to the variations in the induced pseudo-current density, corresponding to the planarstreamfunction (29), required to satisfy the boundary condition at the shield wall. Schematic approximations ofthe surface currents for the m = 0 zonal and m = 1 tesseral harmonic responses can be seen in Fig. 7. For theboundary condition on the cylindrical surface to be satisfied, the induced azimuthal pseudo-currents must mirrorazimuthal current paths on the planar coil surfaces, so that the associated fields cancel in the region ρ s > ρ > ρ c .Since zonal responses are composed of simple circular loops, which can be formed in either a cylindrical or planarbasis, the response of the magnetic shield enhances the magnetic field in the optimization region. Consequently, theshield amplifies the axial magnetic field by a factor of . at the shield’s center, as shown in Fig. 5b–c. The uniformfield gradient therefore exhibits superior fidelity because the continuum response of the passive shield approximatesa distributed cylindrical coil. The resulting system, composed of both the coil and shield, completely encloses theinterior region, producing a high-fidelity magnetic field gradient. The tesseral responses are more complicated, withthe cylindrical surface of the magnetic shield acting to oppose the magnetic field generated by the planar system in theregion ρ < ρ s due to the formation of saddle-type currents [38]. These currents result from the required continuity ofthe induced azimuthal pseudo-current density in a cylindrical basis. Due to the restricted current flow on the planarsurface, the only way to mitigate the opposing field generated by the cylindrical surface of the shield is to minimizethe magnetic field at the boundary, resulting in field coil designs that are oscillatory, as seen in Fig. 3a. The responseof the passive shield, for this configuration, reduces the magnetic field by a factor of . at the shield’s center, asshown in Fig. 3b–c.The shield’s response not only makes the optimization of tesseral harmonic fields more difficult in regions close to thecylindrical surface of the magnetic shield, but also has an effect on the level of fidelity that can be achieved over anyregion when the coils are in close proximity to the magnetic shield. To demonstrate this, we investigate the fidelityof the constant transverse field generated by optimized bi-planar coils, similar to those in Fig. 3, within a magneticshield whose radius varies over the range ρ s = [ ρ c , ρ c ] . In Fig. 8a we plot field profiles at selected shield radii toshow how the uniformity of the transverse field along the x -axis improves as the radius of the shield increases, asexpected, due to the reduced inhomogeneities introduced by the cylindrical surface of the magnetic shield at distancesfurther from the field coils. To evaluate the deviation in the field uniformity, we can calculate the mean RMS error,which represents the averaged deviation of the transverse field from perfect uniformity. In Fig 8b, we evaluate themean RMS error in the transverse field along the z -axis of the optimization region as the radius of the magnetic shieldincreases. The uniformity of the optimized fields in Fig. 8a and Fig. 8b initially decreases rapidly as the radius of themagnetic shield increases, but approaches the wide shield limit when ρ s ∼ ρ c , at which point the field generated bythe cylindrical surface of the magnetic shield becomes negligible. In the wide shield limit the predominant contributionfrom the magnetic shield can then be approximated through the use of the method of mirror images from two infiniteparallel perfect magnetic conductors.Although the pseudo-currents do not exist physically, analysis of the magnetic fields generated by the shield in responseto the planar coils gives insight into how such coils should be designed in order to best realize a specified target field. Forexample, generating magnetic fields that require coil geometries with m -fold azimuthal symmetry causes the cylindricalsurface of the magnetic shield to oppose the magnetic field generated by the planar coil. Consequently, if a tesseral fieldis required over a large radial region, the distance from this region to the planar coils and from the cylindrical shieldsurface should be minimized and maximized, respectively. The opposite is true for the fields generated by the planarend caps of the shield. When the bi-planar coils are located near the end caps, where − L s / < − z (cid:48) < z < z (cid:48) < L s / ,but radially distant from the cylindrical surface of the magnetic shield, the reflected pseudo-currents generated bythe method of mirror images provide a field similar to that of the planar coils, resulting in a field that is magnified.Consequently, bi-planar coils are desirable in magnetic shields with large radii and small aspect ratios, L s /ρ s < .In comparison, the magnetic field generated by the shield in response to coils contained on a cylindrical surface [33]shows the opposite effect. Tesseral fields generated by cylindrical coils are enhanced by the cylindrical surface and arereduced by the planar end caps, favoring long magnetic shields with L s /ρ s > . Due to the conflicting conditions onthe geometries (i.e. planar or cylindrical) of the coils and the magnetic shield, a system composed of coupled planarand cylindrical coils may have advantages for generating desired tesseral field profiles with the greatest accuracy.Intrinsic inaccuracies in calculating the optimum current density and approximating the current continuum by discretewires are not the only sources of error that must be considered when designing fields using our method. Particularly forthe tesseral harmonic field generating coils, wire patterns may be highly meandering, making them hard to manufactureand, due to their high resistance, power consuming. When manufacturing these systems, the achievable field fidelity islimited by the discretization of the continuum current flow pattern and by the creation of magnetic dipoles via theinteraction between separate current streamlines and the wires that follow them. Fortunately, new ways to realize thecurrent continuum are emerging due to advances in foldable PCBs [40] and 3D-printing technologies [41], which enable0more complex coils to be made accurately. It should also be noted that approximating a perfect magnetic shield by amaterial of finite permeability µ r = 20000 and thickness d = 1 mm only introduces small deviations in the model of . % [33], which is much less than the error introduced in the desired field by the coil designs. As a result, given anaccurate representation of the current continuum, our methodology can be used reliably to generate target magneticfields in high-permeability environments. The python code used to design arbitrary the planar coils in a magneticallyshielded cylinder is openly available from GitHub and can be cited at http://doi.org/10.5281/zenodo.4442661 [42].Verification using COMSOL Multiphysics ® requires a valid license. z = L s / z = − L s / z = z (cid:48) z = − z (cid:48) z = z (cid:48) / z = − z (cid:48) / ρ s ρ c ρ c FIG. 2: Schematic diagram showing the areas occupied by the bi-planar coils (light gray) inside a cylindrical closedhigh-permeability magnetic shield (black outline). The high-permeability shield is of length L s = 1 m and radius ρ s = 0 . m, with planar end caps located at z = ± L s / . The shield encloses interior conducting planes of radius ρ c = 0 . m which are located at z (cid:48) = ± . m. The optimization region (red, enclosed by hatched lines) is boundedalong the z -axis between the coil planes with top and bottom positions z = ± z (cid:48) / , respectively, and extends radiallybetween ρ = [0 , ρ c / .1FIG. 3: Wire layouts (a) and performance (b–d) of a hybrid active–passive system optimized to generate a constanttransverse magnetic field, B x . Current flows on the surface of two disks of radius ρ c = 0 . m, which are separatedsymmetrically from the origin and lie at z (cid:48) = ± . m. The wire layouts are optimized to generate a constanttransverse field, B = 1 µ T, across the cylinder and normal to its axis of symmetry. The current-carrying planes areplaced symmetrically inside a perfect closed magnetic shield of radius ρ s = 0 . m and length L s = 1 m and themagnetic field is optimized between ρ = [0 , ρ c / and z = ± z (cid:48) / , as shown in Fig. 2. The least squares optimizationwas performed with parameters N = 50 , M = 1 , β = 1 . × − T /W, t = 0 . mm, and ρ = 1 . × − Ω m. (a)Color map of the optimal current streamfunction calculated for the upper current-carrying plane in Fig. 2. Blue andred shaded regions correspond to current counter flows and their intensity shows the streamfunction magnitude fromlow (white) to high (intense color). Solid and dashed black curves represent discrete wires with opposite senses ofcurrent flow, approximating the current continuum with N ϕ = 12 global contour levels. Streamfunction on the lowerlight gray plane in Fig. 2 is geometrically identical but the current direction is reversed. (b) B x , calculated versustransverse position, x , for y = z = 0 , from the current continuum in (a) in three ways: analytically using (31)-(33)(solid red curve); numerically using COMSOL Multiphysics ® Version 5.3a, modeling the high-permeability cylinder asa perfect magnetic conductor (blue dotted curve); numerically without the high-permeability cylinder and using theBiot–Savart law with N ϕ = 100 contour levels (dashed green curve). Black lines enclose the regions where thecalculated field deviates from the target field by (dashed) and (dot-dashed). (c) B x , calculated versus axialposition, z , for x = y = 0 , from the current continuum in the same three ways as for (b). (d) Enlarged section of (c)showing good agreement between the numerical and analytical results throughout the optimization region.2FIG. 4: Color maps showing the magnitude of the magnetic field, in the xz -plane inside a closed finite length perfectmagnetic conductor generated by the active–passive system depicted in Fig. 3: (a) transverse component, B x , and (b)axial component, B z . The field profiles were calculated numerically using COMSOL Multiphysics ® Version 5.3a. Themagnetic field is optimized between ρ = [0 , ρ c / and z = ± z (cid:48) / ; dashed red lines in (a–b). Contours, in (a) only, showwhere the field deviates from the target field by % (dashed curves) and % (dot-dashed curves).3FIG. 5: Wire layouts (a) and performance (b–c) of a hybrid active–passive system optimized to generate a linearvariation of B z with z -position. Current flows on the surface of two disks of radius ρ c = 0 . m, which are separatedsymmetrically from the origin and lie at z (cid:48) = ± . m. The wire layouts are optimized to generate a linear axial fieldgradient, d B z / d z = 2 µ T/m, along the z -axis of the cylinder through the harmonic B = G ( − x ˆx − y ˆy + 2 z ˆz ) . Thecurrent-carrying planes are placed symmetrically inside a perfect closed magnetic shield of radius ρ s = 0 . m andlength L s = 1 m and the magnetic field is optimized between ρ = [0 , ρ c / and z = ± z (cid:48) / , as shown in Fig. 2. Theleast squares optimization was performed with parameters N = 50 , M = 0 , β = 1 . × − T /W, t = 0 . mm, and ρ = 1 . × − Ω m. (a) Color map of the optimal current streamfunction calculated for the upper current-carryingplane in Fig. 2. Intensity of red shaded regions shows the streamfunction magnitude from low (white) to high (intensecolor). Solid and dashed black curves represent discrete wires with opposite senses of current flow, approximating thecurrent continuum with N ϕ = 12 global contour levels. Streamfunction on the lower light gray plane in Fig. 2 isgeometrically identical but the current direction is reversed. (b) Transverse magnetic field, B x , calculated versustransverse position, x , for y = z = 0 , from the current continuum in (a) in three ways: analytically using (31)-(33)(solid red curve); numerically using COMSOL Multiphysics ® Version 5.3a, modelling the high-permeability cylinder asa perfect magnetic conductor (blue dotted curve); numerically without the high-permeability cylinder and using theBiot–Savart law with N ϕ = 100 contour levels (dashed green curve). Black lines enclose the regions where thecalculated field deviates from the target field by (dashed) and (dot-dashed). (c) Axial magnetic field, B z ,calculated versus axial position, z , for x = y = 0 , from the current continuum in (a) in the same three ways as for (b).4FIG. 6: Color maps showing the magnitude of the magnetic field, in the xz -plane inside a closed finite length perfectmagnetic conductor generated by the active–passive system depicted in Fig. 5: (a) transverse component, B x , and (b)axial component, B z . The field profiles were calculated numerically using COMSOL Multiphysics ® Version 5.3a. Themagnetic field is optimized between ρ = [0 , ρ c / and z = ± z (cid:48) / ; dashed red lines in (a–b). Contours show where thefield deviates from the target field by % (dashed curves) and % (dot-dashed curves).5...... ( a ) z = − L s / z = L s / stuf f z = ( − p z (cid:48) + pL s z = ( − p z (cid:48) + pL s z = z (cid:48) z = z (cid:48) p = 0 p = − p = 1 ...... ( b ) stuf f p = 0 p = − p = 1 FIG. 7: Schematic diagram showing how the approximate zonal (a) and tesseral (b) harmonic responses are generatedby simple current loops with cylindrical and m -fold azimuthal symmetries, respectively. The discrete current loops(black) are located at z = z (cid:48) and z = z (cid:48) with their p th reflected pseudo-image-currents (red) generated by the surfaceof the cylindrical shield at ρ = ρ s and by the planar end caps at z = ± L s / , in accordance with the method of imagesdescribed by (29). Images of the two planar (black) coils resulting from the end caps are located at z = ( − p z (cid:48) + pL s and z = ( − p z (cid:48) + pL s , respectively, where p ∈ Z (two such image coils are shown blue). For all odd reflections theaxial current direction is reversed.6FIG. 8: Improvement in the performance of a hybrid active–passive system, optimized to generate a constanttransverse field B = 1 µ T, upon increasing the radius of the passive shield. Current flows on the surface of two disksof radius ρ c = 0 . m, which are separated symmetrically from the origin and lie at z (cid:48) = ± . m. The wire layoutsare optimized to generate constant B = 1 µ T between ρ = [0 , ρ c / and z = ± z (cid:48) / , as shown in Fig. 2. Thecurrent-carrying planes lie within the bore of a perfect closed magnetic shield of radius ρ s = [ ρ c , ρ c ] and length L s = 1 m. Each least squares optimization was performed with parameters N = 50 , M = 1 , β = 1 . × − T /W, t = 0 . mm, and ρ = 1 . × − Ω m. (a) Transverse magnetic field, B x , calculated analytically versus transverseposition, x , using (31)-(33) (red curves), where: ρ s /ρ c = 1 . (dotted curve), ρ s /ρ c = 1 . (dash-dotted curve), ρ s /ρ c = 1 . (dashed curve), and ρ s /ρ c = 2 . (solid curve). (b) Root mean square (RMS) deviation, ∆ B RMS x , betweenthe calculated and target fields evaluated along the axis of the optimized field region and plotted as a function of ρ s /ρ c . Light blue crosses show ∆ B RMS x values calculated analytically using (31)-(33). Horizontal dashed red lineshows the analytical value of ∆ B RMS x = 0 . obtained in the wide shield limit ( ρ s /ρ c >> ).7 IV. CONCLUSION
Here we have formulated an analytical model to calculate the magnetic field generated by an arbitrary static currentdistribution confined to a plane whose normal is parallel to the axis of a finite closed high-permeability cylinder. Ourformalism is based on a Green’s function expansion of the vector potential generated by a planar coil system. Tosatisfy the boundary conditions at the surface of the cylinder, we assume that it is a perfect magnetic conductor andexpress the induced magnetization in terms of a pseudo-current density. Due to the shared azimuthal symmetries ofthe planar current flows and the induced pseudo-current on the cylindrical surface of the magnetic shield, the responseof the magnetic shield can be expressed in terms of the planar current distribution. We used this formalism to enablea priori optimization of the planar current distribution required to generate specified target magnetic field profiles bycombining a Fourier–Bessel decomposition of the current flow with a quadratic minimization procedure.We demonstrated this optimization methodology by using it to design bi-planar coils that accurately generate aconstant transverse field, B x , across the cylinder, and a linear gradient field, d B z / d z , along the axis of the cylinder.Predictions from our analytical model agree well with with subsequent forward finite element simulations, validatingour methodology. We quantified the interaction of planar systems with a closed finite high-permeability cylinder andfound that all coil designs without complete azimuthal symmetry induced magnetic fields in the cylindrical surface ofthe shield that opposed the field generated by the planar current source. Conversely, planar end caps have the oppositeeffect and amplify the field from the planar current source.Our analytical model and optimization procedure extends the range of current geometries that can be tailored withinfinite closed cylindrical magnetic shields in order to accurately generate specified target field profiles. This enables thedevelopment of systems that require the best magnetic field control that can be achieved subject to size, weight, powerand cost constraints. It may also enable hybrid active/passive planar coils to be retrofitted to existing cylindricalmagnetic shields in order to improve their performance for a low cost and without disrupting existing experimentalsetups. In the future, combining planar and cylindrical coils could enable ultra-precise magnetic fields to be generatedthrough even larger interior fractions of magnetic shields.The PYTHON code used to design arbitrary planar coils in a magnetically shielded cylinder is openly available fromGitHub [42]. Verification using COMSOL Multiphysics requires a valid license. ACKNOWLEDGMENTS
We acknowledge support from the UK Quantum Technology Hub Sensors and Timing, funded by the Engineering andPhysical Sciences Research Council (EP/M013294/1).The authors M.P., P.J.H., T.M.F., M.J.B., and R.B. declare that they have a patent pending to the UK GovernmentIntellectual Property Office (Application No.1913549.0) regarding the magnetic field optimization techniques describedin this work. N.H, M.J.B, and R.B also declare financial interests in a University of Nottingham spin-out company,Cerca. The authors have no other competing financial interests.
APPENDIX A: MAGNETIC FIELD DERIVATION
In order to calculate the magnetic field, we must first determine both the Fourier transform of the streamfunction andthe induced Fourier pseudo-current density. Substituting (30) into (17) and integrating over the azimuthal angle yields ϕ m ( k ) = ρ c N (cid:88) n =1 M (cid:88) m (cid:48) =0 ( W nm (cid:48) ( δ mm (cid:48) + δ m, − m (cid:48) ) + iQ nm (cid:48) ( δ mm (cid:48) − δ m, − m (cid:48) )) × (cid:90) ρ c d ρ (cid:48) ρ (cid:48) J m ( kρ (cid:48) ) J m (cid:48) (cid:18) ρ nm (cid:48) ρ (cid:48) ρ c (cid:19) . (46)Upon substituting (46) into the expressions for the magnetic field, (25)-(27), the δ − functions collapse the infinitesummation, resulting in every Bessel function of order m being replaced by one of order m (cid:48) . Hence, without loss of8generality, applying a Bessel function product identity (5.54.1 in [39]), we may now write ϕ m ( k ) = ρ c N (cid:88) n =1 M (cid:88) m (cid:48) =0 ( W nm (cid:48) ( δ mm (cid:48) + δ m, − m (cid:48) ) + iQ nm (cid:48) ( δ mm (cid:48) − δ m, − m (cid:48) )) × ρ nm (cid:48) k ρ c − ρ nm (cid:48) J m (cid:48) ( kρ c ) J (cid:48) m (cid:48) ( ρ nm (cid:48) ) . (47)Similarly, we may also calculate the induced Fourier pseudo-current density. Substituting (30) into (29) and integratingover the azimuthal coordinate yields J mpφ ( k ) = ρ c e − ik ( ( − p z (cid:48) + pL s )2 ρ s | k | I m ( | k | ρ s ) K (cid:48) m ( | k | ρ s ) N (cid:88) n =1 M (cid:88) m (cid:48) =0 ( W nm (cid:48) ( δ mm (cid:48) + δ m, − m (cid:48) ) + iQ nm (cid:48) ( δ mm (cid:48) − δ m, − m (cid:48) )) × (cid:90) ∞ d˜ k ˜ k ˜ k + k J m (˜ kρ s ) (cid:90) ρ c d ρ (cid:48) ρ (cid:48) J m (cid:48) (˜ kρ (cid:48) ) J m (cid:48) (cid:18) ρ nm (cid:48) ρ (cid:48) ρ c (cid:19) . (48)To simplify these expressions, we evaluate the integral over ˜ k first and then separate the integrand using partialfractions to yield (cid:90) ∞ d˜ k ˜ k ˜ k + k J m (˜ kρ s ) J m (˜ kρ (cid:48) ) = (cid:90) ∞ d˜ k ˜ kJ m (˜ kρ s ) J m (˜ kρ (cid:48) ) − k (cid:90) ∞ d˜ k ˜ k ˜ k + k J m (˜ kρ s ) J m (˜ kρ (cid:48) ) . (49)Inspecting (49) we see that the first component is simply the orthogonality relation between Bessel functions, and thesecond has a known solution (6.541.1 in [39]). Hence, we can state that (cid:90) ∞ d˜ k ˜ k ˜ k + k J m (˜ kρ s ) J m (˜ kρ (cid:48) ) = δ ( ρ (cid:48) − ρ s ) ρ (cid:48) − k I m (cid:48) ( | k | ρ (cid:48) ) K m (cid:48) ( | k | ρ s ) . (50)However, since ρ (cid:48) < ρ s , the first term in (50) can be neglected. Inserting (50) into (48) and integrating over ρ (cid:48) resultsin the expression J mpφ ( k ) = ρ c k e − ik ( ( − p z (cid:48) + pL s )2 ρ s | k | I m ( | k | ρ s ) K (cid:48) m ( | k | ρ s ) N (cid:88) n =1 M (cid:88) m (cid:48) =0 ( W nm (cid:48) ( δ mm (cid:48) + δ m, − m (cid:48) ) + iQ nm (cid:48) ( δ mm (cid:48) − δ m, − m (cid:48) )) × ρ nm (cid:48) k ρ c + ρ nm (cid:48) J (cid:48) m (cid:48) ( ρ nm (cid:48) ) I m (cid:48) ( | k | ρ c ) K m (cid:48) ( | k | ρ s ) . (51)Substituting (47) and (51) into (25)-(27), we can express the summation over the exponentials associated with theinfinite reflections of the planar streamfuction as ∞ (cid:88) p = −∞ e − k | z − ( − p z (cid:48) + pL s | = e − k | z − z (cid:48) | + 2 e kL − (cid:0) e kL cosh ( k ( z + z (cid:48) )) + cosh ( k ( z − z (cid:48) )) (cid:1) , (52)and ∞ (cid:88) p = −∞ ( z − ( − p z (cid:48) + pL s ) | z − ( − p z (cid:48) + pL s | e − k | z − ( − p z (cid:48) + pL s | = ( z − z (cid:48) ) | z − z (cid:48) | e − k | z − z (cid:48) | − e kL − (cid:0) e kL sinh ( k ( z + z (cid:48) )) + sinh ( k ( z − z (cid:48) )) (cid:1) . (53)Similarly, the summation over the infinite pseudo-current reflections can be written as a Fourier series expansion, ∞ (cid:88) p = −∞ e ik ( z − ( − p z (cid:48) + pL s ) = πL ∞ (cid:88) p = −∞ δ (cid:16) k − pπL (cid:17) (cid:16) e ik ( z + z (cid:48) − L ) + e ik ( z − z (cid:48) ) (cid:17) . (54)9Using these expansions we may finally write B ρ ( ρ, φ, z ) = µ ρ c N (cid:88) n =1 M (cid:88) m =0 ρ nm J (cid:48) m ( ρ nm ) ( W nm cos ( mφ ) + Q nm sin ( mφ )) B nmρ ( ρ, z ) , (55) B φ ( ρ, φ, z ) = µ ρ c ρ N (cid:88) n =1 M (cid:88) m =0 mρ nm J (cid:48) m ( ρ nm ) ( W nm sin ( mφ ) − Q nm cos ( mφ )) B nmφ ( ρ, z ) , (56) B z ( ρ, φ, z ) = µ ρ c N (cid:88) n =1 M (cid:88) m =0 ρ nm J (cid:48) m ( ρ nm ) ( W nm cos ( mφ ) + Q nm sin ( mφ )) B nmz ( ρ, z ) , (57)where B nmρ ( ρ, z ) = − (cid:90) ∞ dk k σ ( k ; z, z (cid:48) , L s ) J (cid:48) m ( kρ ) J m ( kρ c ) k ρ c − ρ nm − ∞ (cid:88) p =1 ˜ p | ˜ p | λ p ( z, z (cid:48) , L s ) I (cid:48) m ( | ˜ p | ρ ) I m ( | ˜ p | ρ c ) K m ( | ˜ p | ρ s ) I m ( | ˜ p | ρ s )( | ˜ p | ρ c + ρ nm ) , (58) B nmφ ( ρ, z ) = (cid:90) ∞ d k kσ ( k ; z, z (cid:48) , L s ) J m ( kρ ) J m ( kρ c ) k ρ c − ρ nm + ∞ (cid:88) p =1 ˜ p λ p ( z, z (cid:48) , L s ) I m ( | ˜ p | ρ ) I m ( | ˜ p | ρ c ) K m ( | ˜ p | ρ s ) I m ( | ˜ p | ρ s )( | ˜ p | ρ c + ρ nm ) , (59) B nmz ( ρ, z ) = (cid:90) ∞ d k k γ ( k ; z, z (cid:48) , L s ) J m ( kρ ) J m ( kρ c ) k ρ c − ρ nm − ∞ (cid:88) p =1 ˜ p τ p ( z, z (cid:48) , L s ) I m ( | ˜ p | ρ ) I m ( | ˜ p | ρ c ) K m ( | ˜ p | ρ s ) I m ( | ˜ p | ρ s )( | ˜ p | ρ c + ρ nm ) , (60)where γ ( k ; z, z (cid:48) , L s ) = e − k | z − z (cid:48) | + 2 e kL s − (cid:20) e kL s cosh ( k ( z + z (cid:48) )) + cosh ( k ( z − z (cid:48) )) (cid:21) , (61) σ ( k ; z, z (cid:48) , L s ) = ( z − z (cid:48) ) | z − z (cid:48) | e − k | z − z (cid:48) | − e kL s − (cid:20) e kL s sinh ( k ( z + z (cid:48) )) + sinh ( k ( z − z (cid:48) )) (cid:21) , (62) λ p ( z, z (cid:48) , L s ) = 2 L s (( − p sin (˜ p ( z + z (cid:48) )) + sin (˜ p ( z − z (cid:48) ))) , (63) τ p ( z, z (cid:48) , L s ) = 2 L s (( − p cos (˜ p ( z + z (cid:48) )) + cos (˜ p ( z − z (cid:48) ))) , (64)and ˜ p = pπ/L s .0 APPENDIX B: POWER DISSIPATION
The power dissipated in the surface of a circular planar current source of thickness t and resistivity (cid:37) is given by P = (cid:37)t (cid:90) ρ c d ρ (cid:48) ρ (cid:48) (cid:90) π d φ (cid:48) | J ρ ( ρ (cid:48) , φ (cid:48) ) | + | J φ ( ρ (cid:48) , φ (cid:48) ) | . (65)We can substitute the streamfunction (30) into the continuity relations (12) to obtain J ρ ( ρ (cid:48) , φ (cid:48) ) = ( H ( ρ (cid:48) − ρ c ) − H ( ρ (cid:48) )) ρ c ρ (cid:48) N (cid:88) n =1 M (cid:88) m =0 mJ m (cid:18) ρ nm ρ (cid:48) ρ c (cid:19) ( W nm sin( mφ (cid:48) ) − Q nm cos( mφ (cid:48) )) , (66) J φ ( ρ (cid:48) , φ (cid:48) ) = ( H ( ρ (cid:48) − ρ c ) − H ( ρ (cid:48) )) N (cid:88) n =1 M (cid:88) m =0 ρ nm J (cid:48) m (cid:18) ρ nm ρ (cid:48) ρ c (cid:19) ( W nm cos( mφ (cid:48) ) + Q nm sin( mφ (cid:48) )) . (67)Inserting (66)-(67) into (65) and integrating over the azimuthal component results in P = (cid:37)t πρ c N (cid:88) n =1 M (cid:88) m =0 (cid:2) (1 + δ m ) W nm + (1 − δ m ) Q nm (cid:3) × (cid:90) d˜ ρ (cid:34) m ˜ ρ J m ( ρ nm ˜ ρ ) − ρ nm ˜ ρ J m − ( ρ nm ˜ ρ ) J m +1 ( ρ nm ˜ ρ ) (cid:35) , (68)which, when integrated over the radial component, can be expressed as P = (cid:37)t πρ c N (cid:88) n =1 W n ρ n J ( ρ n ) , (69)for m = 0 , and P = (cid:37)t πρ c N (cid:88) n =1 M (cid:88) m =1 (cid:0) W nm + Q nm (cid:1) (cid:16) ρ nm (cid:17) m m !( m − (cid:34) ˜ F (cid:18) m, m + 12 ; m + 1 , m + 1 , m + 1; − ρ nm (cid:19) − ρ nm m + 1) ˜ F (cid:18) m + 12 , m + 1 , m + 1; m, m + 2 , m + 2 , m + 1; − ρ nm (cid:19) (cid:35) , (70)for m ∈ Z + , where i ˜ F j is the regularized hypergeometric function. Some specific evaluations of (70) are: m = 1 : P = (cid:37) t πρ c N (cid:88) n =1 (cid:0) W n + Q n (cid:1) ρ n J ( ρ n ) , (71) m = 2 : P = (cid:37) t πρ c N (cid:88) n =1 (cid:0) W n + Q n (cid:1) (cid:34) (cid:0) ρ n − (cid:1) J ( ρ n ) − (cid:0) ρ n − (cid:1) ρ n J ( ρ n ) J ( ρ n )+ (cid:0) ρ n − (cid:1) ρ n J ( ρ n ) (cid:35) , (72)1 m = 3 : P = (cid:37) t πρ c N (cid:88) n =1 (cid:0) W n + Q n (cid:1) (cid:34) (cid:0) ρ n − (cid:1) ρ n J ( ρ n ) − (cid:0) ρ n − (cid:1) ρ n J ( ρ n ) J ( ρ n )+ (cid:0) ρ n − ρ n + 96 ρ n − (cid:1) ρ n J ( ρ n ) (cid:35) . (73) APPENDIX C: EXAMPLE COIL DESIGNS
FIG. 9: Wire layouts (a) and performance (b) of a hybrid active–passive system optimized to generate to generate alinear variation of B x with z -position. Current flows on the surface of two disks of radius ρ c = 0 . m, which areseparated symmetrically from the origin and lie at z (cid:48) = ± . m. The wire layouts are optimized to generate a lineartransverse field gradient, d B x / d z = 1 µ T/m, along the z -axis of the cylinder through the harmonic field profile B = ( z ˆx + x ˆz ) . The current-carrying planes are placed symmetrically inside a perfect closed magnetic shield ofradius ρ s = 0 . m and length L s = 1 m and the magnetic field is optimized along the z -axis between z = ± z (cid:48) / . Theleast squares optimization was performed with parameters N = 100 , M = 1 , β = 1 . × − T /W, t = 0 . mm, and ρ = 1 . × − Ω m. (a) Color map of the optimal current streamfunction calculated for the upper current-carryingplane in Fig. 2 [blue and red shaded regions correspond to current counter flows and their intensity shows thestreamfunction magnitude from low (white) to high (intense color)]. Solid and dashed black curves represent discretewires with opposite senses of current flow, approximating the current continuum with N ϕ = 12 global contour levels.Streamfunction on the lower light brown plane in Fig. 2 is geometrically identical but the current direction is reversed.(b) Transverse magnetic field, B x , calculated versus axial position, z , from the current continuum in (a) in three ways:analytically using (31)-(33) (solid red curve); numerically using COMSOL Multiphysics ® Version 5.5a, modeling thehigh-permeability cylinder as a perfect magnetic conductor (blue dotted curve); numerically without thehigh-permeability cylinder and using the Biot–Savart law with N ϕ = 100 contour levels (dashed green curve). Blacklines enclose the regions where the axial field gradient deviates by % (dashed) and % (dot-dashed).To further validate our design process, in Figs. 9–10 we show additional bi-planar hybrid active–passive systems andanalyze their behavior. These systems are designed using our open-access Python code [42]. The coordinate axes andmagnetic field plots are labeled in the same way as the systems presented in the main text. For both systems, the coilslie in the z (cid:48) = ± . m planes, as in the main text and shown in Fig. 2, and the fields are optimized along the z -axisbetween z = ± z (cid:48) / . The design in Fig. 9 generates linear transverse field gradient, d B x / d z , whose spatial variationmatches the harmonic B = ( z ˆ x + x ˆ z ) between two identical plates of radius ρ c = 0 . m, inside a magnetic shieldof radius ρ s = 0 . m and length L s = 1 m. The field gradient is amplified by a factor of . at the shield’s center.If the shield is lengthened, the field gradient is diminished to . at the shield’s center in the long shield limit, asexpected from our analysis in the main body of the text. The design in Fig. 10 generates a uniform field, B z , in the2FIG. 10: Wire layouts (a–b) and performance (c) of a hybrid active–passive linear system optimized to generate aconstant axial field, B z . Current flows on the surface of two disks of radii ρ c = 0 . m (on the upper current-carryingplane in Fig. 2) and ρ c = 0 . m (on the lower current-carrying plane in Fig. 2), respectively, which are separatedsymmetrically from the origin and lie at z (cid:48) = ± . m. The wire layouts are optimized to generate a constant axialfield, B = 1 µ T, along the cylinder and parallel to its axis of symmetry. The current-carrying planes are placedsymmetrically inside a perfect closed magnetic shield of radius ρ s = 1 m and length L s = 1 m and the magnetic field isoptimized along the z -axis between z = ± z (cid:48) / . The least squares optimization was performed with parameters N = 100 , M = 0 , β = 1 . × − T /W, t = 0 . mm, and ρ = 1 . × − Ω m. (a–b) Color maps of the optimalcurrent streamfunction on the upper and lower planar disks, respectively [blue and red shaded regions correspond tocurrent counter-flows and their intensity shows the streamfunction magnitude from low (white) to high (intense color)].Solid and dashed black curves represent discrete wires with opposite senses of current flow, approximating the currentcontinuum with N ϕ = 6 global contour levels. (c) Axial magnetic field, B z , versus axial position, z , calculated fromthe current continuum in (a–b) in three ways: analytically using (31)-(33) (solid red curve); numerically usingCOMSOL Multiphysics ® Version 5.5a, modeling the high-permeability cylinder as a perfect magnetic conductor (bluedotted curve); numerically without the high-permeability cylinder and using the Biot–Savart law with N ϕ = 100 contour levels (dashed green curve). Black lines enclose the regions where the axial field gradient deviates by %(dashed) and % (dot-dashed). z -direction and extending between two planes of different sizes, with upper and lower coil sizes ρ c = 0 . m and3 ρ c = 0 . m, respectively, inside a magnetic shield of radius ρ s = 1 m and length L s = 1 m. Since the shield andcoils have a relatively wide form factor, the field is highly uniform, with the region of axial field deviation below %extending over more than % of the distance between the planes. REFERENCES [1] X. Wu, Z. Pagel, B. S. Malek, T. H. Nguyen, F. Zi, D. S. Scheirer, and H. Müller, Gravity surveys using a mobile atominterferometer, Sci. Adv. , 9 (2019).[2] V. Menoret, P. Vermeulen, N. L. Moigne, S. Bonvalot, P. Bouyer, A. Landragin, and B. Desruelle, Gravity measurementsbelow − g with a transportable absolute quantum gravimeter, Sci. Rep. , 1 (2018).[3] M. J. Snadden, J. M. McGuirk, P. Bouyer, K. G. Haritos, and M. A. Kasevich, Measurement of the Earth’s gravity gradientwith an atom interferometer-based gravity gradiometer, Phys. Rev. Lett. , 971 (1998).[4] K. Bongs, M. Holynski, J. Vovrosh, P. Bouyer, G. Condon, E. Rasel, C. Schubert, W. P. Schleich, and A. Roura, Takingatom interferometric quantum sensors from the laboratory to real-world applications, Nat. Rev. Phys. , 731 (2019).[5] V.K Shah and R.T. Wakai, A compact, high performance atomic magnetometer for biomedical applications, Phys. Med.Biol , 8153–8161 (2013).[6] E. Boto et al. , Moving magnetoencephalography towards real-world applications with a wearable system, Nature ,657–661 (2018).[7] R.M. Hill et al. , Multi-channel whole-head OPM-MEG: Helmet design and a comparison with a conventional system,NeuroImage , 116995 (2020).[8] T.M. Tierney et al. , Optically pumped magnetometers: From quantum origins to multi-channel magnetoencephalography,NeuroImage , 598–608 (2019).[9] A. Borna et al. , Non-Invasive Functional-Brain-Imaging with an OPM-based Magnetoencephalography System, PLOSONE , 1–24 (2020).[10] H. Eswaran, D Escalona-Vargas, E.H. Bolin, J.D. Wilson, and C.L. Lowery, Fetal magnetocardiography using opticallypumped magnetometers: a more adaptable and less expensive alternative?, Prenat Diagn. , 193–196 (2017).[11] S. Lew, M.S. Hämäläinen, and Y. Okada, Toward noninvasive monitoring of ongoing electrical activity of human uterus andfetal heart and brain, Clinical Neurophysiology , 2470–2481 (2017).[12] S. Strand, W. Lutter, J.F. Strasburger, V.k. Shah, O. Baffa, and R.T. Wakai, Low-Cost Fetal Magnetocardiography: AComparison of Superconducting Quantum Interference Device and Optically Pumped Magnetometers, Journal of theAmerican Heart Association , 16 (2019).[13] L. Liang et al. , In-orbit operation of an atomic clock based on laser-cooled Rb atoms, Nature Comms. , 16 (2018).[14] T.E. Chupp, P. Fierlinger, M.J. Ramsey-Musolf, and J.T. Singh, Electric dipole moments of atoms, molecules, nuclei, andparticles, Rev. Mod. Phys. , 015001 (2019).[15] I. Altarev et al. , A magnetically shielded room with ultra low residual field and gradient, Review of Scientific Instruments , 075106 (2014).[16] J. Li, W. Quan, B. Han, Z. Wang, and J. Fang, Design and Optimization of Multilayer Cylindrical Magnetic Shield forSERF Atomic Magnetometer Application , 1793–1800 (2020).[17] J. Li, Z. Wang, and W. Quan, Multi-objective optimization of multilayer passive magnetic shield based on genetic algorithm,AIP Advances , 125210 (2019).[18] S. Celozzi, R. Araneo, and G. Lovat, Appendix B: Magnetic Shielding, John Wiley & Sons Ltd. (2008).[19] S. Pissanetzky, Minimum energy MRI gradient coils of general geometry, Meas. Sci. Technol. , 7 (1992).[20] M. Poole and R. Bowtell, Novel gradient coils designed using a boundary element method, Concept Magn. Reson. B ,162–175 (2007).[21] R. Zettera, A.J. Mäkinen, J. Iivanainen, K.C.J. Zevenhoven, R.J. Ilmoniemi, and L. Parkkonen, Magnetic field modelingwith surface currents. Part II. Implementation and usage of bfieldtools, J. Appl. Phys , 063905 (2020).[22] A.J. Mäkinen, R. Zetter, J. Iivanainen, K.C.J. Zevenhoven, L. Parkkonen, and R.J. Ilmoniemi, Magnetic field modelingwith surface currents. Part I. Implementation and usage of bfieldtools, J. Appl. Phys , 063906 (2020).[23] C.-Y. Liu, T.Andalib, D.C.M. Ostapchuk, and C.P. Bidinosti, Analytic models of magnetically enclosed spherical andsolenoidal coils, Nucl. Instrum. Methods Phys. Res. , 162837 (2020).[24] R. Lambert and C. Uphoff, Magnetically shielded solenoid with field of high homogeneity, Rev. Sci. Instrum. , 3 (1975).[25] R.J. Hanson and F.M. Pipkin, Magnetically shielded solenoid with field of high homogeneity, Rev. Sci. Instrum. , 2(1965).[26] J.D. Jackson, Classical electrodynamics, (Wiley, New York, 1998) 3rd ed.[27] A. Kordyuk, Magnetic levitation for hard superconductors, J. Appl. Phys , 610–612 (1985).[28] Q. Cao, D. Pan, J. Li, Y. Jin, Z. Sun, S. Lin, G. Yang, and L. Li, Optimization of a coil system for generating uniformmagnetic fields inside a cubic magnetic shield, Energies, , 3 (2018).[29] T. Liu, J. Voigt, Z. Sun, A. Schnabel, K. Grueneberg, I. Fan, and L. Li, Two-step mirror model for the optimization of themagnetic field generated by coils inside magnetic shielding, Conference on Precision Electromagnetic Measurements, 1–2(2018). [30] T. Liu, A. Schnabel, Z. Sun, J. Voigt, and L. Li, Approximate expressions for the magnetic field created by circular coilsinside a closed cylindrical shield of finite thickness and permeability, J. Magn. Magn. Mater. , 166846 (2020).[31] N. Holmes, T. M. Tierney, J. Leggett, E. Boto, S. Mellor, G. Roberts, R.M. Hill, V. Shah, G.R. Barnes, M.J. Brookes,and R. Bowtell, Balanced, bi-planar magnetic field and field gradient coils for field compensation in wearable magnetoen-cephalography, Sci. Rep. , 14196 (2019).[32] R. Turner and R.M. Bowley, Passive screening of switched magnetic field gradients, J. Phys. E Sci. Instrum. , 876 (1986).[33] M. Packer, P.J. Hobson, N. Holmes, J. Leggett, P. Glover, M.J. Brookes, R. Bowtell, and T.M. Fromhold, Optimal InverseDesign of Magnetic Field Profiles in a Magnetically Shielded Cylinder, Phys. Rev. Applied , 054004 (2020).[34] P. Hammond, Electric and magnetic images, Proc. IEE Part C Monogr. , 306 (1960).[35] W.A. Roshen, Effect of finite thickness of magnetic substrate on planar inductors, IEEE Trans. Magn. , 270 (1990).[36] J.W. Carlson, K.A. Derby, K.C. Hawryszko, and M.Weideman, Design and evaluation of shielded gradient coils, Magn.Reson. Med. , 2, (1992).[37] A.N. Tikhonov, On the stability of inverse problems, Dokl. Akad. Nauk SSSR , 195–198, (1943).[38] F. Roméo and D.I. Hoult, Magnet field profiling: analysis and correcting coil design, Magn. Reson. Med. , 1 (1984).[39] I.S. Gradshteyn and I.M. Ryzhik, Table of integrals, series, and products, Academic Press, Amsterdam , 1 (2007).[40] J. R. Corea, A. M. Flynn, B. Lechene, G. Scott, G. D. Reed, P. J. Shin, M. Lustig, and A. C. Ariasb, Screen-printed flexibleMRI receive coils, Nat. Commun. , 10839 (2016).[41] T. D. Ngo, A. Kashani, G. Imbalzano, K. T. Q. Nguyen, and D. Hui, Additive manufacturing (3D printing): A review ofmaterials, methods, applications and challenges, Compos. B. Eng.143