A Three-step Model for Optimizing Coil Spacings Inside Cuboid-shaped Magnetic Shields

A three-step model for calculating the magnetic field generated by coils inside cuboid-shaped shields like magnetically shielded rooms (MSRs) is presented. The shield is modelled as two parallel plates of infinite width and one tube of infinite height. We propose an improved mirror method which considers the effect of the parallel plates of finite thickness. A reaction factor is introduced to describe the influence of the vertical tube, which is obtained from finite element method (FEM) simulations. By applying the improved mirror method and then multiplying the result with the reaction factor, the magnetic flux density within the shielded volume can be determined in a fast computation. The three-step model is verified both with FEM and measurements of the field of a Helmholtz coil inside an MSR with a superconducting quantum interference device. The model allows a fast optimization of shield-coupled coil spacings compared to repetitive time-consuming FEM calculations. As an example, we optimize the distance between two parallel square coils attached to the MSR walls. Measurements of a coil prototype of 2.75~m in side length show a magnetic field change of 18~pT over the central 5~cm at the field strength of 2.7~\textmu T. This obtained relative field change of 6~ppm is a factor of 5.4 smaller than our previously used Helmholtz coil.


I. INTRODUCTION
Extremely uniform magnetic fields with a high temporal stability are essential to various precision measurements, such as electric dipole moment (EDM) measurements 1-4 or magnetic field detector calibration [5][6][7] . For the challenging EDM measurements a nonzero magnetic field gradient enhances vibration noise seen by magnetometers 8,9 and shortens the possible measurement time for a single experiment, thus deteriorating its statistical sensitivity [10][11][12] . The final unavoidable field gradient also causes systematic uncertainties like geometric phase shift 13,14 , which is the dominant error contribution to the present neutron EDM upper limit 15 . More details about the influence of the magnetic field gradient can be found in Refs. 16 and 17. The common method of creating a homogeneous magnetic field is to place a coil set into a multi-layer magnetic shield [18][19][20][21] which serves to reduce external field perturbations. A number of classic coil configurations have been proposed for generating homogeneous magnetic fields in air, such as the Helmholtz coil or the solenoidal coil 22,23 . However, the high permeability shielding material alters the field distribution and normally worsens the expected uniformity. By optimizing the coil spacings, the ferromagnetic shield can increase the strength and also the uniformity of the magnetic field produced by coils compared to the same setup placed in an unshielded environment 18,21,[24][25][26] . The result of spacing optimizations depends on the accuracy of the field calculation including the distortion caused by the shielding material. To date the field generated by coils inside a non-spherical magnetic shield can be solved analytically only for some ideal a) Electronic mail: silasliutianhao@gmail.com b) Electronic mail: 23hnhosava@163.com situations, such as assuming an infinite length or an infinite permeability of the shield 21,24,25 . The real shield is often of a cuboid shape with no known analytical solution.
For realistic setups, the finite element method (FEM) allows to calculate the magnetic field distribution, but the calculation time required for an accurate description is still quite long for present computers, due to a large length-to-thickness ratio of the shielding material. Sweeping the coil parameters in FEM models to find an optimal solution takes even longer. For this reason, it is expedient to use a reasonable and sufficiently accurate simplification of the field calculation in order to optimize coil spacings in an efficient way.
For a coil in front of an infinitely large and infinitely thick shielding plate, an accepted treatment is substituting the plate by an imaginary coil on the opposite side of the plate as a reflection of the original coil, known as the method of mirror images 27,28 . Pan et al. applied the standard images method to optimize square coils enclosed in a magnetically shielded room (MSR) 29 . For the case where the thickness of the plates is several mm, the current of the imaginary coil is, however, influenced by the geometry of the setup as well as the position of the observation point, which is not considered in the method of mirror images. The reduced image current due to the finite width of the ferromagnetic plate has been accounted for in the calculation of several specific cases via FEM 30 . The accuracy was significantly improved, but the effects of the finite thickness were still neglected.
In this paper, we propose a three-step model to calculate the magnetic field generated by coils inside a cuboid-shaped magnetic shield, which is suitable for rapid optimizations of the coil spacings. The proposed three-step model is verified with a complete FEM calculation and a measurement of a 1.6-mdiameter Helmholtz coil placed inside the Berlin Magnetically Shielded Room 2 (BMSR-2) 31 . Furthermore, the model is applied to optimize the distance between two square coils attached to the walls of BMSR-2 with respect to the field homo-geneity in the central area of the MSR. The three-step model turns out to be 10 times faster than the pure FEM optimization and delivers exactly the same optimal distance. For the newly designed coil pair prototype of 2.75 m in side length with the calculated optimal distance, the measured maximal relative magnetic field change in a 5 cm region around the local extremum near the center of the coil is 6 ppm, which is an improvement by a factor of 5.4 compared to the value of our previously used Helmholtz coil.

II. THREE-STEP MODEL DESCRIPTION
It is assumed that the windings of the coil set are restricted to parallel planes, e.g. Helmholtz coils or Braunbek coils 32,33 . The six faces of a cubic shield as in Fig. 1(a) are divided into two parts, two horizontal plates parallel to the coil plane and a vertical tube comprising all side walls, as illustrated in Fig. 1(b) and Fig. 1(c).
According to the images method, the two horizontal plates of high permeability are treated as mirrors and create multiple mirror images of the coil and also of the vertical tube, thereby extending the tube to a square tube of infinite length as shown in Fig. 1(d). In our model we assume that the mirror tubes have the same permeability as the original tube, and that the effect of the horizontal plates and of the vertical tube on the field distribution can be considered separately. Therefore, the field generated by a coil inside an MSR of finite thickness and finite permeability can be determined in three steps: Step 1: Calculate the magnetic flux density B air , generated by the coil in air, with the Biot-Savart law.
Step 2: Calculate the coordinate r i and the current ratio τ * i of the ith mirror coil, which is introduced to replace the effect of the two parallel plates.
Step 3: Calculate the reaction factor curve η k in dependence of the position with an FEM simulation once, which is introduced to describe the effect of the surrounding vertical tube in relation to the coil in air.
We choose a coordinate system with the origin in the center of the magnetic shield, and the position of the original loop is denoted as r 0 = (x 0 , y 0 , z 0 ). The final expression for the k = . . . x, y or z component of the magnetic flux density at a position r = (x, y, z) inside the enclosure given by the three-step model is where r i = r − r i . The sign of i indicates the location of the mirror coil, i.e., positive and negative mean above and below the original coil. The sum in Eq. (1) is over the original coil (i = 0) and all mirror coils. Without loss of generality, in the later analysis we focus on the z component of the magnetic field since the considered coil is positioned in the xy-plane and thus B z is the most important component regarding field uniformity. The other components can be derived accordingly.
For convenience, we also let the axis of the coil be colinear with the z axis, leading x 0 = y 0 = 0. The principle of the three-step model can be applied to a coil with any shape in one plane. Here, the two most frequently used coil shapes, a circular coil and a square coil, are used as examples to illustrate our model.

III. STEP 1: COILS IN AIR
The magnetic flux density B z created by a circular coil loop carrying a current I in air at the observation point r(x, y, z) with x 2 + y 2 < a 2 is where µ 0 is the permeability of the vacuum, a is the radius of the coil, z 0 is the z coordinate of the coil plane and J 0(1) is the zeroth (first) order Bessel function of the first kind 34 . Here we use the solution with Bessel functions instead of the more familiar equivalent expression involving elliptic integrals 35 . We need this solution for the comparison done in section IV B. A square coil located in the plane z = z 0 is composed out of four identical current segments of equal length L coil . The center of the jth wire segment is denoted as S j (x j , y j , z 0 ). For the segment with a current in x direction, the magnetic flux density is solved according to the Biot-Savart law as 36 where ρ 2 = (y−y j ) 2 +(z−z 0 ) 2 , l c = L coil /2 and I 1 = −I 3 = I. The x coordinate and y coordinate in Eq. (3) have to be interchanged for segments carrying current in y direction, leading to where ρ 2 = (x − x j ) 2 + (z − z 0 ) 2 , I 2 = −I 4 = I. The total magnetic field generated by a square coil is the sum of the contributions of all four wire segments In step 2, the impact of the two parallel plates on the field distribution is analyzed. It is necessary to obtain the z coordinate z i and current ratio τ * i of the ith mirror coil generated by the two parallel plates.
A. Classical mirror model for parallel plates The classical idea of mirror images for a wire in front of one permeable plate is shown in Fig. 2(a). The effect of the highpermeability plate on magnetic fields can be substituted by the effect of a mirror wire at the same distance to the plate surface as the original wire but on the opposite side. If the thickness ∆ of the plate is large enough, the value of the mirror current I m can be written as 27 where µ r is the relative permeability of the plate, and I 0 is the source current. The magnetic field B a on the upper side of the plate in Fig. 2(a) is the sum of B 0 , the field generated by the source current I 0 , and B m , generated by the mirror current I m .
Mimicking the transportation of light in two opposite mirrors, a wire in between two parallel plates that are treated as mirrors will be repetitively reflected, thus generating multiple mirror wires 29,30 . The ith mirror wire is labelled as I i . The absolute value |i| is the number of reflections needed to create this mirror wire. One of the two reflection roads starting from the original wire can be described as: I 0 → I −1 → I 2 → I −3 , . . . , and is illustrated in Fig. 2(b). The second road, not shown, combines the other reflections. The current I i for the ith mirror wire is and the z coordinate of the ith mirror wire is where z 0 is the z coordinate of the original wire and L denotes the distance between the two plates. Herein, the origin of the z coordinate is in the middle of the two plates. Eqs. (7) and (8) are independent of the shape of the original current path and thus are valid for circular coils, too. The total field is the sum of the contributions from the original wire and all the mirror wires. In general, the multiple mirror method provides an intuitive and easy-to-calculate solution to the considered problem. However, the accuracy of this model is only guaranteed for infinitely thick plates. For a finite-thickness plate, the value of the mirror current is reduced compared to Eq. (6).

B. Analytical calculation of the mirror current ratio for plates of finite thickness
To calculate the magnetic field of coils inside a practical MSR, the assumption of an infinite thickness of the shielding material is hard to justify. The material is only a few mm thick. Together with realistic values of µ r of 15000 -30000 (the effective µ r for BMSR-2 is 17500 31 ) the influence of the finite thickness on the field could reach several percent. Therefore, our model includes the influence of the finite thickness.
In multiple mirror models the parallel plates are treated sequentially and each time only one plate is considered as a mirror. Therefore, the analysis with a single plate is applicable to obtain the current value I i for all mirror wires.
Consider a coil loop of radius a and current I 0 located at z 0 = 0 in the xy-plane. A plate of thickness ∆ is located at z = −h. We choose a Cartesian coordinate system whose origin is placed in the center of the loop. The solution for the magnetic flux density in the region z > −h and x 2 + y 2 < a 2 is solved in

Ref. 34 as
with The first term B air z in Eq. (9) is the field generated by the loop in air (see Eq. (2)) so that the second term is the impact of the permeable plate, noted as B where B air,z m z is the field generated by an imaginary air coil of the same radius a and current I 0 as the original coil, but located at z m = −2h (the mirror position of the original coil). τ is a proportionality factor. Substituting B plate z and B air z (Eq. (9)) into Eq. (11), one obtains (12) τ is independent of the source current I 0 and it can also be interpreted as the current ratio of the mirror loop to the original loop. It can be deduced from Eq. (12) that τ < 1, which means that the mirror current value is always smaller than the source current. In the limit ∆ → ∞, T defined in Eq. (10) reduces to T = (µ r − 1)/(µ r + 1), which is independent of the integral variable ζ . Inserting this term into the definition of τ in Eq. (12), one finds that τ = T = (µ r − 1)/(µ r + 1), which is identical to the constant mirror current ratio defined in the classical images method (see Eq. (6)).
After the transformation in Eq. (11), the magnetic field can be written as B z = B air z + τB air,z m z , which is an alternative to the analytic expression in Eq. (9). The analytical expression for the magnetic field generated by a coil in front of a plate does not exist for most of coil geometries, such as the square coil loop. In this case, the mirror method can still be applied as an approximate solution via taking the current ratio τ of a circular loop as that of other coil geometry. Fig. 3 shows τ as a function of an observation point on the z-axis and the yaxis. For a plate of µ r = 20000 and ∆ = 2 mm, τ is about 5% smaller than the ratio in the classical mirror method, which is 0.9999. It was found via FEM simulations that τ changes by less than 1% if the square loop of side length L coil = 1 m is approximated by an internal tangent circular loop (a = 0.5 m) or external tangent circular loop (a = 0.707 m). For the further calculation, we always choose the average of these two radii as the radius for an equivalent circular loop, that is r equ = L coil (1 + √ 2)/4 substituting a in Eq. (12). The current ratio defined in Eq. (12) for a plate of finite thickness can be also applied to the multiple mirror method, which is used to model a coil in between two parallel plates. For the ith mirror loops, which are created by reflecting the original loop |i| times, the overall current ratios are the successive product of |i| independent current ratios as where the z coordinate of the nth mirror loop z n is given in Eq. (8) and τ i is the current ratio for a coil located at z i in front of a single plate, as defined in Eq. (12). Compared to the classical mirror current ratio in Eq. (7), τ * i takes the finitethickness of plates into consideration, thus named as finitethickness current ratio.The magnetic flux density at the observation point (x, y, z) between the two parallel plates considering 2M mirror loops is where i = 0 for the original loop and τ * 0 = 1. For a mirror loop with an increasing order i the current ratio τ * i decreases while the distance to the observation point |z − z i | increases, making its contributions to the total magnetic field in Eq. (14) decrease rapidly. For a certain accuracy, this infinite sum can be truncated and the resulting numerical error will be discussed in section VI.

V. STEP 3: COILS INSIDE A VERTICAL TUBE
The effect of the high-permeability infinite tube is quantitatively characterized by a reaction factor η , which is a function of the observer position (x, y, z) inside the tube and can be where B s and B air are the z component of the magnetic flux density generated by the coil with the tube and without the tube (in air). η > 1 means the magnetic field strength is augmented by the high-permeability tube at the observation point. The definition of η is independent from the location of the coil. Due to the difficulties in analytically determining B s , we applied FEM simulations in COMSOL v5.4 to obtain η. In the FEM model, the height of the tube is set 10 times larger than its width to approximate the infinite length. The shielding tube is modelled as a surface assigned with the magnetic shielding boundary condition with a certain thickness ∆. The usage of this boundary condition avoids the meshing in the direction of the thickness of the material which is a problem due to the large length-to-thickness ratio. By assuming a coil located in the xy-plane at the position z 0 =0, η can be calculated faster using the symmetries (η(x, y, z 0 ) = η(−x, y, z 0 ) = η(x, −y, z 0 ) = η(x, y, −z 0 )) with a one-eighth model in COMSOL as described in Ref. 37. For a coil located at any other z 0 position, η can be obtained with a coordinate transformation. The FEM result for a square coil in the xy-plane at z 0 = 0 with 2.75 m side length inside a high-µ rectangular tube of 3.2 m side length is plotted in Fig. 4 as a function of the observer position. The left axis shows B air and B s while the right axis shows η . The same current produces a larger field in the coil center when the coil is inside a shield ( η > 1), as expected due to the reduction of the magnetic resonance. For z > 2 m, η becomes less than 1 because the surrounding high-permeability tube attracts most of the flux lines, leading to a strong decrease between 2 m and 6 m. For z > 6 m, η is already less than 0.1.
Once η is calculated by FEM for enough sample points, η can be interpolated to any observation point without using FEM again. Since the tube is infinitely long, any coil of the Step 1 Step 2 FEM same shape in the xy-plane obeys the same η behavior when the parameter z of the η curve is replaced by relative distance z , the distance between z of the observation point and the z coordinate of the coil center. This means that for calculating the field from the original coil and the mirror coils, we only need to calculate η once by FEM. This is also true for a coil set with identical windings on different z planes, like a Helmholtz coil. For different coil shapes the η curve must be recalculated.
Multiplying the magnetic field generated by each coil in Eq. (14) with the corresponding reaction factor η yields the final expression for the magnetic flux density B z within the shielded volume observed at (x, y, z) where z i = z − z i . Note that the expression for the reaction factor η and the current ratio τ * i are given here only for the B z component. For the B x and B y components, they must be calculated accordingly.

A. Comparison to FEM calculations
Here we demonstrate the three-step model for the case of two square coil loops of side length L coil = 2.75 m enclosed in a cubic single-layer MSR of side length L MSR = 3.2 m, a thickness ∆ = 4 mm and a constant µ r = 30000. The square coils are in the xy-plane, as illustrated in the inset of Fig. 5. The distance between the two coils is set as d coil = 0.5445L coil to compose a square Helmholtz coil in air 38 . Fig. 5 shows B z for the single steps of the three-step model and the result for the respective FEM calculation in depen- dence of the z-coordinate, i.e. moving along the coil central axis, where the magnetic field B z of the coil has the strongest change. The infinite sum over all contributing images in Eq. (16) was truncated at i = 5, i.e. we stopped after the fifth mirror image on both sides. The red curve in Fig. 6 shows that the contribution of the next term (i = 6) for this coil arrangement is already less than 0.05% for z = 0 m. We checked that even for z up to 0.3 m, the difference to the z = 0 curve is less than 5%. The FEM calculation (done with COMSOL v5.4) is regarded as the benchmark, even though it has a finite uncertainty because it is a mathematical approximation which depends on the chosen boundary condition and the mesh size. Using optimized configurations, the uncertainty is expected to be < 0.1%.
Step 1 is the exact solution of the coil in air (Eq. (5)), which underestimates the field inside the shield by 33.2%. Adding the influence of the two walls parallel to the coils as infinite plates with finite thickness leads to a field 1.2% below the FEM result. The complete three-step model follows well the increasing FEM trend for larger z values and underestimates the magnetic field by only about 0.1%, which is already close to the expected uncertainty of the FEM used as benchmark. To demonstrate the influence of introducing the finitethickness current ratio (Eq. (13)) instead of classical current ratios (Eq. (7)) in step 2, we also calculated this case, marked with classical mirror method in Fig. 5. The assumption of an infinite thickness leads to an expected field overestimation. The difference of 0.3% to the FEM benchmark is small for this case but increases significantly when ∆ or µ r are smaller. For example, for the same setup with µ r = 10000 and 2000, the differences are 1.0% and 6.1%.
The final deviation of the three-step model is a sum of various effects. Effects leading to an overestimation are: infinite width assumption for the parallel plates, infinite height assumption for the side tube and neglected coupling between the parallel plates and the side tube. In contrast truncation TABLE I. The relative error of the field derived from the three-step model compared to the FEM results at the origin with varied parameters values. The setup is the same as in Fig. 5 with default values µ r = 30000, d coil = 1.5 m and L MSR = 3.2 m. In each column, only one parameter value is altered and the others are left at their default value. Since we would like the coils to be attached to inside walls of the shield, the length of the coil is changed when the length of the MSR is changed according to L coil = 0.86L MSR . of the sum over all contributions leads to an underestimation. The use of a circular coil as an approximation for the square coil to calculate τ * and the error of the reaction factor η derived via FEM simulations may lead to a positive or a negative error depending on the given geometry. The relative error of the three-step model compared to the FEM results over wide ranges of parameter values for the same setup as in Fig. 5 is calculated and summarized in Table I. The error never exceeds 1.28% for all considered cases. µ r has the greatest influence, but for a realistic value of µ r = 30000 (permalloy) the error is less than 0.32% if only the coil distance d coil or the length of the MSR L MSR is varied. The maximal influence was found for µ r = 2000 which is a suitable value for ferrites.

B. Comparison to magnetic field measurements
A measurement to further validate the three-step model was conducted. A movable circular 3-axes Helmholtz coil was used to generate a uniform magnetic field inside BMSR-2. BMSR-2 has a passive shielding factor greater than 7 × 10 5 at 0.01 Hz and an effective relative permeability of 17500 31 . The innermost 4-mm permalloy layer surrounds a volume of 3.2 × 3.2 × 3.216 m 3 . For the three-step model and FEM calculations we included only the innermost shielding layer and neglected all others because the field enhancing effect should be dominated by the innermost layer. Due to this assumption and in order to take care of the additional field increase caused by the remanence of the shielding material, which is not included in the model, we set µ r = 30000.
The Helmholtz coil used was oriented along the y-axis of BMSR-2 and has a diameter of 1.6 m. The setup is shown in Fig. 7(a) and (b). The coil was powered by a low-noise current source from Magnicon 39 . The current was measured to be 38.6 mA for each coil with 60 turns in series, generating a 2.6 µT magnetic field in the center measured in a previous nuclear magnetic resonance experiment using 3 He spin polarizations. The magnetic flux densities were measured with a 304-channel vector Superconducting Quantum Interference Device (SQUID) system, mounted in a two-axis ball-bearing slideway. The overall absolute positioning error of the SQUID system is ±1 cm. The possible movement range of the SQUID system along the y axis inside the 3-axes Helmholtz coil is from -4 cm to 5 cm. The position of the SQUIDs was moved in steps of 1.0 ± 0.1 cm. The SQUID's position was locked by a pneumatic brake after each movement. In each position, we measured the magnetic field for 10 seconds to average the noise. Although there are numerous available SQUID channels at different locations, we only used the SQUID channel with the minimum distance to the coil center. In our case, y7 was used.
In the beginning of the field measurements, the y7 SQUID sensor was aligned with the coils axis (y-axis) by tilting and rotating the SQUID system to detect the maximum field signal. The remanent field of BMSR-2 in the nT range was measured after degaussing without the coil field and subtracted in data post-processing. Since the DC offset of the SQUID system is unknown, only the change of the magnetic field is plotted in Fig. 8. The change of the created magnetic field from the field minimum (at -0.4 cm) to a position 5 cm further (at 4.6 cm) is around 98 pT, corresponding to a relative change of 38 ppm. The calculated field change of this Helmholtz coil in air is more uniform, shown as a grey dotted line in Fig. 8. This demonstrates how the surrounding high permeability walls distort the magnetic field. The results calculated with the three-step model and with the FEM are in acceptable agreement with the measured values. The slight difference may result from the uncertainty of the position of the SQUID chip, and the fact that µ r is not a constant 40 as was assumed both for FEM calculations and for the three-step model. Nevertheless, both models successfully predict the maximum of the field at y = 6 cm caused by the coupling to the shielding material. The reason for the This convinced us that the threestep model could be used for the design and optimization of shield-coupled coils. A successful practical example is presented in the next section.

VII. APPLICATIONS
Some of the authors have been using BMSR-2 within a collaboration to measure the 129 Xe EDM using the free spin pre- cession of Xe gas 41 . This experiment has stringent requirements on the uniformity of the magnetic field in the region of the sample cell with a size of about 5 cm × 5 cm × 5 cm.
According to the requirement of the EDM experiment, we defined the magnetic field uniformity in the region of interest as where B y (y) and σ (B y (y)) are the average value and the standard deviation of the magnetic field series B y (y). Here we used a range for y from 0 (coil center) to 5 cm in steps of 1 cm, meaning 6 points in B y (y). U is a unitless quantity independent of the absolute field strength in the center. We built two large square y-axis coils. They were attached to the MSR walls as shown in Fig. 7(c) in order to test if the coils could be permanently hidden behind the internal lining in future. This kind of built-in coil could provide an improved magnetic field homogeneity even though they are located at only a few cm distance to the shielding material. Compared to the Helmholtz coil used before, this built-in coil set would expand the applicable space of the chamber and save time spent in assembling and disassembling the Helmholtz coil system every time measurements with additional magnetic fields are performed. The side length of the square coils was L coil = 2.75 m and the distance d between the two coils is the parameter to be optimized. A single coil consists of 8 turns in series, and in total the coil set had a resistance of 18.8 Ω.
In order to find the optimal distance d of the coil we first applied the three-step method to obtain the uniformity U as a function of the distance d in the range of 1.25 m to 1.45 m, shown as the light blue curve in Fig. 9(a). Within 15 minutes of CPU time, we obtained an optimal value of d = 1.35 m. A following FEM simulation was carried out in a narrower range (1.32 m, 1.39 m), shown as the dark blue dashed curve in Fig. 9(a). After 20 hours of CPU time the optimal distance using repetitive FEM calculations differs by less than 1 mm from the value acquired by the three-step model. The threestep model is not able to reproduce the sharp resonance-like increase of U obtained from the FEM calculation. The main reason is that the reaction factor curve η in such a small region of 5 cm is uncertain as η is obtained in a FEM model with a tube height over 30 m, leading to a coarse mesh. This uncertainty is also responsible for the unsmooth curve obtained by the three-step model in Fig. 8. Fig. 9(a) shows that the position of the coil is critical and should be adjustable around the optimum position.
The coil pair for the experimental confirmation was set up with a distance of d = 1.35 m ± 0.02 m. A current of 0.42 A through the coil pair generated a magnetic field of 2.7 µT in the center, measured with a fluxgate MAG-03 from Bartington Instruments. Before the current was switched on, we measured the remanent field along the y-axis, shown as the grey dotted line in Fig. 9(b). This measured remanent field with a gradient around 3 pT/cm was subtracted from the measured field when the current to the coils was on. Before each measurement, the BMSR-2 was degaussed in order to reduce the remanent field of the shield 42 , providing reproducible starting conditions after the change of the spacing distance between the coils.
The measured residual-field-corrected curve of the coil in a 9 cm range is shown for three coil distances d in Fig. 9(b). The curve with d = 1.35 m (blue) corresponds to the expected optimum. To validate the expected strong dependency on the distance, the measurement was repeated for d = 1.32 m and d = 1.38 m. The mechanical accuracy in the alignment of the provisional setup of the coil leads to the observed shift of the field extremum away from y = 0. This shift is different each time the coil distance is changed. To compare the measurements always a 5 cm regions starting from the extremum was chosen. The limited mechanical accuracy of the temporary setup also worsens also the homogeneity so that we expect to underestimate the possible improvement. Even for the final coil setup with a higher mechanical accuracy an  adjustment of the experiment to the location of the most homogenous magnetic field is possible. For the blue curve, the most uniform 5 cm wide region is from -1.6 cm (where the field is strongest) to 3.4 cm with a field change of only 18 pT. This is a relative change of 6 ppm in relation to the maximal field strength of 2.7 µT, which is 5.4 times smaller than that of the previously used circular Helmholtz coil. Comparing the results for the three distances, Fig. 9(b) confirms that the field homogeneity around the optimal value close to d = 1.35 m deteriorates fast. A shift of 3 cm to 1.32 m or 1.38 m increases the field change in a 5 cm region starting from the extremum from 18 pT to 114 pT or 190 pT. The change from a maximum field strength for d = 1.38 m to a local minimum field strength for d = 1.32 m is a known behavior for such a coil arrangement. In our temporary setup it was difficult to move the coils and to align them in parallel, leading to an additional shift of the position of the field extremum as seen in Fig. 9(b). The gradient of the d = 1.35 m curve is already close to the gradient of the residual field (remanence). It is unclear if the residual field contribution is unchanged on this level when the coil is turned on and the inner layer of the shield is degaussed. We subtracted the residual field even though the raw data show a more homogeneous field. By designing a larger and optimized coil set, we significantly improved the uniformity of the coil field, compared to the circular Helmholtz coil. Another advantage of the square coil prototype is that a noise peak at 4 Hz due to vibrations disappeared because the coils were attached to the walls. Note that the distance between the cables and the MSR walls is only a few cm. Close to the cables the shielding material is exposed to strong fields. It even could reach saturation. Therefore, µ r will decrease near the cables, leading to a reduction of the shielding factor. It was measured that when the coil pair is on, the low-frequency shielding factor of the BMSR-2 decreases by around 5%. For the 129 Xe-EDM experiment this is acceptable.

VIII. CONCLUSIONS
The proposed three-step model takes into account the impact of a finite-thickness and finite-permeability magnetic shield on the field generated by internal coils. To deal with the finite-thickness plates parallel to the coil, we developed an improved mirror method as an extension of the classical images method. The three-step model is significantly faster than a pure FEM analysis with still reasonable accuracy, which is important for optimizing coil spacings.
The three-step model was verified with the FEM calculation for a square coil pair inside an MSR. In this case the relative error of the three-step model on the magnetic field strength is in the order of the expected 0.1% accuracy of the FEM calculation used as benchmark. We applied the model to calculate the change of the magnetic field created by a 1.6-m-diameter circular Helmholtz coil inside BMSR-2 and obtained a result consistent with the measurement data.
As an example of coil optimization, the optimal distance between two square coils for the best homogeneity in a defined region along the central axis was calculated. We built a coil prototype of 2.75 m in side length attached to the BMSR-2 walls. The optimal distance was rapidly estimated by the three-step model and verified by FEM calculation as well as by experiment. The experimental result of the prototype coil shows that the optimum distance could be predicted within the geometric uncertainty of the setup. The maximal change of the magnetic field is 18 pT over a 5 cm region near the center at the field strength of 2.7 µT, corresponding to a relative field change of 6 ppm. This is more than a factor of 5 smaller than the value measured for the previously used Helmholtz coil. Due to this result we have already optimized a built-in coil set with four square coils with the three-step method. This coil system will be installed inside BMSR-2 during a planned upgrade with a new innermost shielding layer.