Designing optimal loop, saddle, and ellipse-based magnetic coils by spherical harmonic mapping

Adaptable, low-cost, coils designed by carefully selecting the arrangements and geometries of simple primitive units are used to generate magnetic fields for diverse applications. These extend from magnetic resonance and fundamental physics experiments to active shielding of quantum devices including magnetometers, interferometers, clocks, and computers. However, finding optimal arrangements and geometries of multiple primitive structures is time-intensive and it is challenging to account for additional constraints, e.g. optical access, during the design process. Here, we demonstrate a general method to find these optimal arrangements. We encode specific symmetries into sets of loops, saddles, and cylindrical ellipses and then solve exactly for the magnetic field harmonics generated by each set. By combining these analytic solutions using computer algebra, we can use numerical techniques to efficiently map the landscape of parameters and geometries which the coils must satisfy. Sets of solutions may be found which generate desired target fields accurately while accounting for complexity and size restrictions. We demonstrate this approach by employing simple configurations of loops, saddles, and cylindrical ellipses to design target linear field gradients and compare their performance with designs obtained using conventional methods. A case study is presented where three optimized arrangements of loops, designed to generate a uniform axial field, a linear axial field gradient, and a quadratic axial field gradient, respectively, are hand-wound around a low-cost, 3D-printed coil former. These coils are used to null the background in a typical laboratory environment, reducing the magnitude of the axial field along the central half of the former's axis from $\left(7.8\pm0.3\right)$ $\mu$T (mean $\pm$ st. dev.) to $\left(0.11\pm0.04\right)$ $\mu$T.


I. INTRODUCTION
Evermore sophisticated methods of magnetic field design are required to better null unwanted variations and generate targeted biases in precision measurement devices. Techniques to control magnetic fields have been applied for various precision measurements including in magnetic resonance scanners [1]- [3], electric dipole moment experiments [4]- [6], and Kibble balances for mass determination [7]- [9]. As well as this, precision magnetic field control is a necessity for new generations of quantum-enabled devices, including atomic magnetometers [10]- [15], atom interferometers [16], [17], atom [18], [19] and ion [20] clocks, and quantum computers [21], [22]. Recently, the design of on-board coils for quantum-enabled devices has garnered extensive research interest [23]- [26] as these coils must balance magnetic constraints, such as powerefficiency and field uniformity, with significant non-magnetic considerations, including optical access and ease-of-assembly. Turner [27] pioneered the design of target coils from discretized surface currents, in which wires emulate a continuum of current flowing on a surface. Typically, the surface current is decomposed into a set of weighted orthogonal modes [28]- [30]. As these modes generate spatially orthogonal magnetic fields, the best combinations of weightings to generate a specified target field profile may be determined using simple mathematical programming [31], often least-squares minimization. However, the real-world performance of surface current-based coil systems is limited if the surface current is inaccurately emulated by the wire pattern. This is known as discretization error [32]. This error increases if the target field region is close to the coil or if manufacturing limitations prevent precise emulation of the surface current. In some cases discretization error may be calculated analytically [33], [34], but often it may only be ascertained by numerically simulating each wire pattern a posteriori until a well-represented design is found. This may be very computationally intensive. Unlike surface current-based coils, those directly designed from simple building blocks, including loops, saddles, and ellipses, hereby referred to as primitives, do not suffer from discretization error and have regular shapes, making them cheap and easy to manufacture. Romeo and Hoult [35] developed a method to design magnetic fields using these coils by exploiting the symmetries in a desired magnetic field, expressed using a spherical harmonic decomposition, with sets of primitives. They then maximized the target harmonic strength by tuning both the continuous geometric coil parameters, which determine the shape of the sets of primitives, and the discrete ratios of the number of turns of wire among different sets of primitives. Arrangements of multiple sets of primitives have the capability to produce a high quality magnetic field since each extra primitive provides additional parameters that one can optimize. However, the optimization of multiple primitives is challenging as each primitive generates many spherical harmonics. To design such systems, one may examine Taylor expansion coefficients [36]- [39] of analytic expressions [23], [40], but these do not relate simply to the basis of spherical harmonics in many cases. This can make it difficult to ascertain how the minimization of different undesired variations should be prioritized to maximize field quality. Alternatively, numerically-solved analytic formulations of the magnetic field [41]- [43] have been used effectively to design loops and saddles, but such approaches do not allow the relationship between the coil parameters and field quality to be determined a priori. As a result of these limitations, surface current-based coils are often preferred over primitive-based systems in settings where accurate target fields are required over large target regions relative to the coil size. Here, we facilitate the wider use of primitive-based coil systems by providing a generalized approach to their design. We mathematically encode sets of primitive building blocks with spherical harmonic-like symmetries such that all magnetic field variations in free space generated by the sets may be encoded as exact, closed-form expressions. To achieve this, we impose the sets directly into Turner's surface current solution [27] and apply a spherical harmonic decomposition [35]. We solve the resulting integral equations analytically to determine the expansion coefficients as simple derivatives with respect to the coil parameters. Then, we use Mathematica to construct and simplify the analytic expressions for the expansion coefficients, upon which numerical rootfinding routines search the coil parameter space for solutions that cancel field errors. By interpolating these solutions, we determine a meshed contour in the parameter space on which the solutions lie. We then rank the solutions according to a desired attribute, e.g. field fidelity, power-efficiency (field strength per unit current), inductance, or by practical concerns such as spatial extent, overlaps with access holes, or the distance between sets of primitives. We present three worked examples of this process -one of each using sets of loops, arcs, and ellipses primitives -to demonstrate the scope of our method, and compare each to standard equivalents. To finish, we design and build uniform axial, linear axial gradient, and quadratic axial gradient, field-generating loops on a 3Dprinted cylindrical coil former. We examine how this system can be used to null the background magnetic field in a typical laboratory environment. The analytic model developed in this work was initially presented at the Conference on Precision Electromagnetic Measurements 2022 [44].

A. Field harmonic basis
In free space, the magnetic field can be represented as the gradient of a magnetic scalar potential, B = −∇Ψ. The scalar potential and magnetic field both satisfy Laplace's equation, ∇ 2 B = ∇ 2 Ψ = 0. Here, we express the magnetic scalar potential as the complete basis of real spherical harmonics and so the magnetic field may be expressed in terms of the complete basis of the vector derivatives of each real spherical harmonic [45]. We refer to this as the basis of field harmonics, where each field harmonic, with magnitude f n,m , is denoted by B n,m (r) = ∇h n,m (r) and h n,m (r) is the spherical harmonic of the same order n ∈ Z + and degree m ∈ Z ∈ [−n : n]. Each spherical harmonic may be expressed as h n,m (r, θ, φ) = η n,m r n P n,|m| (cos θ) cos (mφ) sin (|m|φ) , where the upper and lower terms in the right-hand-side bracket denote variations of degree m ≥ 0 and m < 0, respectively. Each spherical harmonic has a zenith dependence described by Ferrer's associated Legendre polynomials, P n,|m| (cos θ).
The harmonics are defined with the standard normalization, on a unitary sphere, r 0 = 1, with and ζ m,m = 2 − δ m,m , where δ m,m is the Kronecker delta function.
Here, we shall only consider cos (mφ)-like variations, for m ≥ 0. Any sin (|m|φ)-like variations, for m < 0, may be generated by rotating the equivalent cos (mφ)-like variation by π/(2|m|). The field harmonic components in Cartesian coordinates, B n,m (r) = X n,m (r)x+Y n,m (r)ŷ+Z n,m (r)ẑ, are presented in appendix A. The axial component is Z n,m (r, θ, φ) = η n,m (n + m) r n−1 × P n−1,m (cos θ) cos (mφ) . (5) The symmetry of the axial component along the z-axis is determined by the parity of (n+m−1) due to its dependence on P n−1,m (cos θ) [46]. Thus, the axial component is symmetric along the z-axis if (n + m − 1) is even and is antisymmetric along the z-axis if (n + m − 1) is odd. The low order field harmonics encode simple variations when expressed in Cartesian coordinates. For example, for n = [1,2] and m = [0, 1] they are as shown in Fig. 1.

B. Harmonic matching
The magnetic field, (1), generated by the surface currents with cos (mφ)-like symmetry flowing on a cylinder of radius ρ c may be separated into contributions where the axial magnetic field has total symmetry (+) or total antisymmetry (−) along the z-axis about the origin, The axially symmetric field harmonics are of order n = 2ν + m + 1 for ν ∈ Z 0+ and degree m ∈ Z 0+ . The axially antisymmetric field harmonics are of order n = 2ν + m for ν ∈ Z + and degree m ∈ Z 0+ . As the magnetic field must relate uniquely to the scalar potential, the spatial forms of the field harmonics are preserved when examining the field generated by the surface current. Thus, the field harmonics must also be present within the wellknown Green's function integral expression for the axial field in free space generated inside a cylindrical azimuthal surface current [27]. We express this in the axially symmetric and antisymmetric cases as where I m (z) and K m (z) represent the modified Bessel functions of the first and second kinds, respectively, of order m. The Fourier transforms of the axially symmetric and antisymmetric azimuthal current flows, J + φ (r ) and J − φ (r ), respectively, are As detailed in appendix B, we can transform the cylindrical variations in equations (12) and (13) into spherical variations like those in the axial field harmonic component, (5). We find Substituting equations (15) and (16) into equations (12) and (13) and grouping terms to match the field harmonic basis, (10) and (11), we find The axial field harmonic component is zero for harmonics of equal order and degree, Z n,n (r) = 0. To derive these magnitudes, we follow the same procedure as above using the transverse field, B x (r), from the well-known solution [27] and the transverse field harmonic component X n,n (r). The resulting integrals which determine the field harmonic magnitudes are the same as equations (17) and (18). Fig. 2: Variation of primitive current patterns (red arrows indicate current flow direction) in the φz-plane on a cylinder of radius ρc to generate the field harmonics optimized in the main text. (a) Axially antisymmetric loops separated axially by 2dc generate field harmonics of order N = 2ν, for ν ∈ Z + , and degree M = 0. (b) An even number of pairs of axially symmetric sets of arcs with one-fold symmetry along φ, separated axially by 2dc1 and 2dc2 and which extend azimuthally by 2φc, may be connected to make saddles which generate field harmonics of order N = 2ν for ν ∈ Z + and degree M = 2µ + 1 for µ ∈ Z + . (c) Axially symmetric sets of ellipses with one-fold symmetry along φ, separated axially by 2dc and which extend axially by a maximum of ec, generate the same field harmonics as (b).

III. PRIMITIVE DESIGN
Now, we encode the azimuthal surface current on example sets of simple primitive building blocks: loops, saddles, and ellipses. We design the sets of primitives to have specific azimuthal and axial symmetries so that only certain field harmonics are present in the magnetic field. We then solve the field harmonic magnitude integrals, (17) and (18), analytically to exactly determine the magnitudes of the remaining field harmonics. These solutions are used to choose the coil parameters of combinations of sets of primitives to generate target field harmonics. In example III-A, we design axially antisymmetric loops to generate the B 2,0 field harmonic and, in example III-B, we design separate axially symmetric sets of saddles and ellipses to generate the B 2,1 field harmonic. The optimization methodology we apply is presented in section IV. The Mathematica programs we use to design each coil are publicly-available for non-commercial use and are stored in the repository listed in reference [47]. The default arguments in each of the examples in the repository return the coil systems optimized in the main text.

A. Linear axial gradient with respect to axial position
The B 2,0 field harmonic is contained within the axially antisymmetric field, (11), and has m = 0 azimuthal symmetry (no φ-dependence). As such, it is generated by axially antisymmetric current flows with m = 0 azimuthal symmetry. These symmetries are matched by pairs of simple loops which carry equal and opposite current, I, and are separated about the origin by an axial distance 2d c (Fig. 2a). We may represent such a pair of loops using the following azimuthal current density: The Fourier transform, (14), of equation (19) is We determine the field harmonic magnitudes by substituting equation (20) into equation (18), where the normalized separation is χ c = d c /ρ c . In appendix C, we analytically solve the class of integrals, including equation (21), which encode the field harmonic magnitudes generated by loops. We find that We design B 2,0 using equation (22) by choosing turn ratios and separations to maximize the f 2,0 field harmonic while minimizing all others. We then proceed by exactly nulling as many leading-order field harmonics as possible, which we define as the lowest order field harmonics generated by the selected primitive current pattern that are not the desired field I: Properties of the standard coils in Figs. 3a, 5a and 7a, designed using the methodology presented in Romeo and Hoult [35], and the optimized equivalents in Figs. 3b, 5b and 7b designed in this work for a coil radius of ρc = 10 cm. The following coil properties are listed: target field harmonic, Bn,m of order n and degree m; the maximum axial dimension, max (dimc), where dimc = dc in the loops and saddles cases, dimc = dc + ec in the ellipses case, dc is the separation of the coil units from the origin, and ec is the axial extent of the ellipses from their centre; the total length of wire in the coil, l, including repeated units with multiple turn ratios; the magnetic field gradient strength, B0, i.e. dBz/dz for B2,0 and dBx/dz for B2,1 along the z-axis at the centre of the coil, per unit current, I; the inductance, L; and the size of the central volume where the target field gradient is generated with less than 1% deviation from target relative to that generated by the standard coil, V 1% . The inductance is approximated numerically from the total magnetic energy [48] using COMSOL Multiphysics ® , assuming that the wire tracks are filamentary. harmonic. This is generally more effective than minimizing as many field harmonics as possible since the leading-order errors have the lowest spatial frequencies.
In this specific case, as expected, the leading-order field harmonic, f 4,0 , is nulled exactly for χ c = √ 3/2 (Fig. 3a). This coil geometry may be referred to as an anti-Helmholtz pair/gradient-field Maxwell coil.
Here, we use Mathematica (as detailed in section IV) to optimize three axially antisymmetric loop pairs, as presented in Fig. 3b. We null [f 4,0 , f 6,0 ] and simultaneously maximize the weighted sum of f 2,0 /f 8,0 , i.e. the ratio of the target field harmonic magnitude to the magnitude of the next-leadingorder-error weighted by the respective turn ratios of each set of primitives. We also constrain all optimized separations to be less than or equal to that of an anti-Helmholtz pair, χ c ≤ √ 3/2, to make the generated design equally as applicable as an anti-Helmholtz pair in settings where the maximum axial extent of the system is limited, e.g. due to experimental equipment.  Table I, in the xz-plane generated by the coils in Fig. 3a and Fig. 3b, corresponding to (a) and (b), respectively, where the coils are of radius ρc. White contours enclose the regions where dBz/dz deviates from perfect linearity by less than 1% (dot-dashed curves).
In Fig. 4, the axial magnetic field generated by an anti-Helmholtz pair is compared to that generated by the optimized coil. The properties of the coils are summarized in Table I. The size of the central volume where B 2,0 is generated with less than 1% error (less than 1% deviation from target; bounded within the central dot-dashed curves in Fig. 4) is a factor of 2.85 greater than that generated by an anti-Helmholtz pair.
However, due to the increased number of turns and alternating turn ratio polarities, compared to the anti-Helmholtz pair, the optimized coil is 1.17 times less power-efficient and has 5.24 times greater inductance. In a scenario where fast current switching is required alongside high gradient linearity, the inductance could be calculated analytically [49] and imposed as an additional minimization condition.

B. Linear transverse gradient with respect to axial position
The B 2,1 field harmonic is contained within the axially symmetric field, (10), and has m = 1 azimuthal symmetry (one line of symmetry in the ρφ-plane). To match this, we optimize sets of symmetric arcs of azimuthal extent 2φ c , which are axially separated about the origin by 2d c and have one-fold azimuthal periodicity (the arc pairs repeat every φ = π with alternating polarity; Fig. 2b). The azimuthal current density that describes one such set of arcs may be represented as where H (x) is the Heaviside step function. The Fourier transform, (14), of equation (23) is Substituting equation (24) into equation (17), we find As with the loops case, we analytically solve the class of integrals, including equation (25), which encode the field harmonic magnitudes generated by arcs in appendix C. We find where the binomial coefficient is n k = n!/ ((n − k)!k!). Now, we use equation (26) to select turn ratios, separations, and extents of four sets of arcs to maximize f 2,1 while minimizing leading-order errors. To obey current continuity, the number of sets of arcs must be even and the current in each set of arcs must correspond to an equal and opposite current in another set of arcs. As a result, we join each set of arcs to another to form double saddles (Fig. 2b). Firstly, we note that equation (26) prohibits all even degrees and shows us that all field harmonics of odd degrees are proportional to sin (mφ c ). We find optimal extents of multiple sets of arcs by substituting sin (mφ c ) for the relevant Chebyshev polynomial [46], and then finding the root of the functions where leading-order degrees are nulled. Here, to null all field harmonics with m = [3,5], the required azimuthal extents of two sets of arcs with unitary turn ratios are φ c = [7π/30, 13π/30]. We then use Mathematica to optimize four sets of arcs (each of which has two nested azimuthal extensions) to null the leading-order errors, [f 4,1 , f 6,1 , f 8,1 ], and maximize the weighted sum of f 2,1 /f 10,1 generated by the four sets of arcs. As with example III-A, we also constrain all axial separations, here to χ c ≤ 2.56, so that the maximum extent is less than or equal to that of a standard design created in reference [35] (Fig. 5a).
(a) (b) Fig. 6: Magnitude of the normalized transverse magnetic field (color scales right), Bx = Bx/B0, where B0 is the magnetic field gradient strength specified in Table I, in the xz-plane generated by the coils in Fig. 5a and Fig. 5b, corresponding to (a) and (b), respectively, where the coils are of radius ρc. White contours enclose the regions where dBx/dz deviates from perfect linearity by less than 1% (dot-dashed curves).
The optimized coil design is presented in Fig. 5b. In Fig. 6 we show that the optimized coil generates B 2,1 much more effectively than the standard coil. The power-efficiency of the optimized coil is 1.47 times greater and it generates B 2,1 with less than 1% error over a central region that is 4.55 times greater in volume (Table I). Compared with the standard coil, however, as with example III-A, the increased complexity of the optimized coil arrangement increases the comparative inductance by a factor of 4.25.
Given the acute changes in flow direction in saddles between current-carrying arcs and axial connecting wires, inaccuracy in saddle construction may generate unwanted high-order field harmonics [35]. To overcome this, we can use smoothlyvarying elliptical wire tracks on the surface of a cylinder instead of saddles to generate the same field harmonics.
Here, we demonstrate this by generating B 2,1 using sets of axially symmetric cylindrical ellipses, which each extend over an axial distance 2e c and are separated about the origin by a central axial distance 2d c > 2e c . The sets of ellipses do not cross, preserving the axial symmetry, and have one-fold azimuthal periodicity (the ellipse pairs repeat every φ = π with alternating polarity; Fig. 2c). The azimuthal current density which traces the path of one such set of ellipses may be represented as . (27) The Fourier transform, (14), of equation (27) is where J m (z) represents the Bessel function of the first kind of order m. Substituting equation (28) into equation (17), we find where the normalized extent is ψ c = e c /ρ c . We solve the class of integrals which include equation (29) analytically in appendix D. We find wherẽ and Q m−1/2 (z) represents a Legendre function of the second kind of half-integer order, (m − 1/2). Now, we use equation (30) to select turn ratios, separations, and extents of two sets of ellipses to maximize f 2,1 while minimizing leading-order errors. Unlike the saddles case, the optimization of orders and degrees cannot be separated. Thus, to maximize a field harmonic using two sets of ellipses, we completely null the leading-order magnitudes, [f 4,1 , f 4,3 , f 6,1 ], and maximize the weighted ratio of f 2,1 /f 8,1 generated by both sets of ellipses. Additionally, we impose (χ c + ψ c ) ≤ 1.32 to limit the maximum system dimension to be less than or equal to half that of the optimized saddle coils. In Figs. 7 and 8 we present schematics and performances of a standard coil design from reference [35] and the optimized design, respectively. As expected, the optimized elliptical coils are effective at generating B 2,1 even with the constraint on system extent. Compared to the standard double saddles (Fig. 5a) and ellipses (Fig. 7a), the optimized elliptical coil has a maximum axial extent 2.22 and 1.96 times smaller, but generates B 2,1 with less than 1% error over central volumes 1.36 and 6.26 times greater in size. However, due to opposite directions of current polarity between the two sets of ellipses in the optimized design, the standard coils are 5.31 and 5.65 times more power-efficient than the optimized elliptical coil.  Table I, in the xz-plane generated by the coils in Fig. 7a and Fig. 7b, corresponding to (a) and (b), respectively, where the coils are of radius ρc. White contours enclose the regions where dBx/dz deviates from perfect linearity by less than 1% (dot-dashed curves).

IV. OPTIMIZATION ROUTINE
Our goal is to find an optimal set of coil parameters which maximize a desired field harmonic while minimizing unwanted contributions. We shall only consider optimization of continuous geometric parameters, not discrete turn ratios which we specify alongside other properties like the minimum separation between primitives. To determine the best turn ratios, we rerun the optimization for simple combinations of turn ratios and then rank solutions across all runs. Each coil parameter affords a degree of freedom with which to null a field harmonic. In theory, it is possible to null as many field harmonics as there are coil parameters, but in non-trivial scenarios these solutions are difficult to find and may not exist within the constrained parameter space. Instead, we consider an underdetermined problem system, i.e. we null fewer field harmonics than there are coil parameters. In the space of coil parameters, the solutions lie on a contour of dimension equal to the number of coil parameters less the number of nulled harmonics. Finding a solution on this contour using numerical techniques is faster and more reliable than finding a single-point solution. Increasing the dimensionality of the contour makes it easier to find solutions at the expense of nulling fewer harmonics. Once solutions are found on a contour, we rank them according to a desired property, such as the ratio of the desired-to-leading-order-error harmonic magnitudes, and then select the best-ranked solution. We implement the optimization in Mathematica because of its fast and easy-to-use symbolic expression simplification (via the Simplify function), numerical root-finding (via the FindRoot function), and process parallelization. Let us consider the design presented in subsection III-A to generate B 2,0 using N loops = 3 pairs of loops. Our optimization algorithm is summarized in Fig. 9.
(a) (b) Firstly, we construct and simplify symbolic expressions for the total magnitudes of the low-order field harmonics generated by the primitive set, weighted by the predetermined turn ratios. The weighted magnitude of each field harmonic is where I i are the turn ratios of each loops pair and the individual field harmonic magnitudes, f n,0 , may be obtained using equation (22). We use equation (32) as the objective function to search for solutions of F 4,0 = F 6,0 = 0 in a coarse mesh of the N loops -dimensional parameter space. In this design, three filtering conditions are applied to the mesh to bound the search space and impose a minimum distance between adjacent coil pairs: χ ci ≥ 0.1, χ ci ≤ √ 3/2, and χ c(i+1) − χ ci ≥ 0.01. This reduces the number of mesh points to 1140 (Fig. 10a). At each mesh point, we apply FindRoot, which uses a Newton-Raphson method with step control [50] to search locally around positions in the solution space for locations where the undesired harmonics are nulled. The search over 1140 points takes 0.65 seconds to complete 1 and 21 unique points are found on the one-dimensional solution contour. We then linearly interpolate between the known points on the solution contour to estimate further solutions. The interpolated points are randomized slightly to expand the evaluation scope, and are then used to seed FindRoot once more. This process occurs at 411 interpolated points and takes 0.18 seconds. This is more efficient than the previous step because the seeds are close to the solution contour and therefore converge rapidly. In this example, all the interpolated seeds converge onto the solution contour, which is shown in Fig. 10b, but in cases where the contour is discontinuous or varies greatly, not all seeds may converge. Finally, we rank the solutions on the contour according to the absolute value of the ratio of the magnitude of the desired field harmonic, F 2,0 , to the leadingorder error field harmonic, F 8,0 . We select the solution that maximizes this ratio (denoted with an arrow in Fig. 10b), which takes 0.01 seconds for the B 2,0 example.

V. CASE STUDY: AXIAL NULLING
Now, we present a case study demonstrating the implementation of optimized magnetic field coils to null residual axial variations in free space, e.g. for residual field compensation for atomic magnetometers [13] or atom interferometers [17]. Fig. 11: Hand-wound 3D-printed coil former containing independent uniform axial, Bz, linear axial gradient, dBz/dz, and quadratic axial gradient, d 2 Bz/dz 2 , field-generating coils.
The system comprises the loop coil in Fig. 3b which generates the B 2,0 field harmonic, along with other loop coils (presented in appendix E) designed using our open-access Mathematica program to generate the B 1,0 and B 3,0 field harmonics. We choose to generate these field harmonics as they are the lowest order m = 0 field harmonics [35], and so are most likely to be present in the background field. The coils are wound by hand using 0.63 mm diameter (22 AWG) enamelled copper wire on a coil former of radius of ρ c = 50 mm (Fig. 11). During optimization, the normalized axial separations are bounded between 0.35 < χ c < 1. This allows the coil former to be only 100 mm long and leaves a large central axial region of length 35 mm where there are no wires. This area contains four optical access holes that are evenly spaced and have a height of 30 mm and an azimuthal extent of π/3 radians. The coil former is manufactured using Fused Deposition Modeling (FDM) with a commercially available UltiMaker S5 using hard PLA body material and PVA support material. The PVA support is dissolved during post-processing by soaking the print in lukewarm water for 48 hours. The FDM printer is capable of rendering grooves in the coil former with a resolution of 0.06 mm and uses cost-effective materials (< £20 in this scenario).  Fig. 12: (a) Magnitude of the axial magnetic field, Bz, per unit current, I, along the z-axis, generated by the uniform (light blue curve), linear gradient (orange curve), and quadratic gradient (light green curve) axial field-generating coils presented in Fig. 13a, Fig. 3b, and Fig. 13b, alongside experimental measurements (black scatter). The magnetic null region is shown with dashed lines (grey); (b) Measured static axial magnetic field magnitude pre-(blue curve) and post-(red curve) active nulling, alongside labels (blue and red) displaying the mean (also shown with dashed black curves) and standard deviation in the axial field magnitude within the null region.
The axial magnetic field generated by each coil is measured using a Stefan Mayer Fluxmaster magnetometer connected to a NI-USB 6212 data acquisition system. During each measurement, sinusoidal currents of amplitude 80 mA and oscillating at a frequency of 5 Hz are passed through each coil in turn for 7 s. The magnetic field produced by each coil is obtained by calculating the Fast Fourier Transform (FFT) of the magnetometer output, eliminating the need for separate measurement offsets. The measured and expected field profiles show excellent agreement and are presented in Fig. 12a. Next, the static axial magnetic field along the central 50 mm of the coil's axis is decomposed into uniform, linear gradient, and quadratic gradient contributions using the method of leastsquares, following the method in Ref. [51]. Coil currents of I = ([178, 59,9] ±1) mA are applied to the B 1,0 , B 2,0 , and B 3,0 field-generating coils, respectively, to null the measured static background. The magnitude of the axial field pre-and post-null is shown in Fig. 12b. Over the central 50 mm of the coil's axis, the mean axial field magnitude reduces by over a factor of 70 from 7.8 µT to 0.11 µT and the standard deviation in the axial field magnitude reduces by over a factor of 7 from 300 nT to 40 nT. To further diminish the axial field, several approaches could be considered, such as adding more field-generating coils to null higher order field harmonics, increasing the number of leading-orders nulled in existing coil designs by increasing the number of separations optimized, or incorporating dynamic feedback from a reference magnetometer to update the applied coil currents [52].

VI. CONCLUSION AND OUTLOOK
Here, we have presented a spherical harmonic decomposition of the magnetic field in free space generated inside sets of primitive structures based on loops, saddles, and ellipses. In each case, we solved for the magnetic field harmonic magnitudes as simple derivatives with respect to the coil parameters. These derivatives can be computed rapidly using computer algebra software, and in this work we used Mathematica. We then demonstrated a generalized approach to design simple primitive coil structures by mapping the landscape of parameters which control the geometry of each set. Using this approach, we designed three field-generating structures -nested sets of loops, arcs, and ellipses -which generate target fields accurately by completely nulling a set number of leading-order undesired contributions. To allow fair comparison with standard systems, in each worked example we chose the solutions closest to the global optimum (which null the leading-order deviations while maximizing the targetfield-to-next-leading-order-error ratio) within the constrained space of coil parameters. However, the optimization procedure could be readily adjusted to prioritize other constraints such as power-efficiency, size, and inductance. Furthermore, more diverse objective functions could also be added to ensure that designs can be manufactured easily, e.g. sensitivity to wire misplacement or number of wire overlaps. Although we have focused on specific worked examples within the main text, the Mathematica program we provide [47] may be used to design any low-order field harmonic using the integral solutions in appendices C and D. We use the loops coil optimized in the main text, alongside other loops coils generated using our program (see appendix E), to construct a 3D-printed, hand-wound axial field compensation system. By calculating the appropriate currents applied to the coil, we were able to significantly reduce both the magnitude and variation in the measured axial field along the central half of the coil's axis in a typical laboratory environment. Specifically, we were able to reduce the axial field magnitude from 7.8 µT to 0.11 µT and its standard deviation from 300 nT to 40 nT. One can design any magnetic field in free space using our program by decomposing the target field into a weighted sum of field harmonics and combining sets of axially symmetric and antisymmetric primitives. Optimized arrangements of primitive coils may be useful in diverse contexts such as the design of accurate and power-efficient linear gradient fields for magneto-optical traps [53] or the cancellation of external interference for atomic magnetometers [51] by using modified governing equations accounting for the response of external passive magnetic shielding [32]. The field harmonics generated by primitives could also be posed and solved for other surface geometries. For example, cuboidal primitives may be useful for generating magnetically-shielded enclosures [54] and toroidal primitives may be applied for plasma confinement in tokamaks [55]. Discrete building block and discretized surface current-based coils could also be implemented together to maximize their performance, e.g. broadband reduction in the noise floor using power-efficient simple building block coils combined with accurate shimming of residual field harmonics using complex surface current-based coil patterns. Alternatively, one may modify and re-solve the governing equations to impose more complex primitives, such as combinations of axially symmetric and antisymmetric units with predetermined turn ratios. This may be useful for applications where a single coil is required to generate a target magnetic field profile composed of multiple field harmonics with different symmetries. Our optimization approach may also be applied in other settings where there are diverse contributions that need to be minimized and maximized simultaneously, including aerofoil design [56] and perturbation analysis in optimization problems [57].
APPENDIX B MATCHING CYLINDRICAL AND SPHERICAL FIELD VARIATIONS Firstly, let us consider the standard series expansions [46], To match the variations in equations (36) and (37) with the arguments of equation (12), x = kρ/ρ c and y = kz/ρ c , we group the terms with common exponents of k, and then convert to spherical coordinates, ρ = r sin θ and z = r cos θ. We find We now introduce the expansion of the associated Legendre polynomials, as derived in reference [58], We substitute x = cos θ into equation (39) and re-index to n = 2ν + m and l = l + m, finding For a fixed value of ν, one may factorize out the variations proportional to sin m θ in equations (38) and (40) by examining the l = l = 0 terms. We can therefore match these variations together, with an unknown scaling, and invert to find the scaling. We find Similarly, following the same method using the standard series expansion of sin(y) [46] and re-indexing equation (39) to n = 2ν + m − 1, we find

WEIGHTING FUNCTION
Here, we obtain the complete class of integrals which we shall refer to as the symmetric and antisymmetric cases, respectively, for ν ∈ Z 0+ and m ∈ Z 0+ . We shall use the standard result [59] ∞ 0 dk k m cos(kχ c )K m (k) = π(2m)! 2 m+1 m!
First, we integrate the left-hand-side of equation (45) by parts, which we can rearrange to show We can substitute equation (45) into this to find Now, we note that the derivative of the second term on the right-hand-side of equation (48) with respect to χ c relates to equation (45) and so To solve equation (50) we perform the substitution χ c = tan θ c , and so cos θ c = 1/ 1 + χ 2 c , finding Changing variables to u = sin θ c , equation (51) can be expanded as where we have applied the binomial theorem, Therefore, equation (50) can be expressed as Then, we note that we can perform derivatives of equation (47) with respect to χ c an odd or even number of times to match equations (43) and (44), Thus, the symmetric and antisymmetric integral cases, (43) and (44), have the same result, where n = 2ν + 1 in the symmetric case and n = 2ν in the antisymmetric case. In the case where m = 0, equation (57) simplifies to
Equation (64) has the known solution [59] S m (χ c , ψ c ) = where Q m−1/2 (z) is a Legendre polynomial of the second kind of half-integer order, m − 1/2, for m ∈ Z + . Now, let us integrate the right-hand-side of equation (64) Table II, in the xz-plane generated by the coils in Fig. 13a and Fig. 13b, corresponding to (a) and (b), respectively, where the coils are of radius ρc. White contours enclose the regions where (a) Bz and (b) d 2 Bz/dz 2 deviate from perfect uniformity and curvature, respectively, by less than 1% and 10%, respectively (dot-dashed curves).  Fig. 13, here given for a coil radius of ρc = 10 cm and labelled as Table I. The magnetic field strength, B0, is the value of Bz for B1,0 and d 2 Bz/dz 2 for B3,0 and is evaluated at the centre of the coil per unit current, I. Sets of symmetric loops are presented in Fig. 13 which are optimized to generate the B 1,0 and B 3,0 field harmonics, displayed in equations (6) and These designs are optimized using our Mathematica program stored in the repository listed in reference [47]. The magnetic fields generated by the coils are displayed in Fig. 14 and their properties are summarized in Table II. All supporting data may be made available on request. Our Mathematica program used to generate the coil designs presented in this work is publicly-available for non-commercial use [47].
This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible.