Active compensation of magnetic ﬁeld distortions based on vector spherical harmonics ﬁeld description

An analytic solution to the magnetostatic inverse problem in the framework of vector spherical harmonic basis functions is presented. This formalism is used for the design of a spherical magnetic ﬁeld compensation system and its performance is compared with an already existing rectangular coil system. The proposed set of spherical coils with 15 degrees of freedom achieves a shielding factor of 1000 or better in a large part of the volume enclosed by the coils for a dipolar type external perturbation. © 2017 Author(s). All article content, except where other-wise


I. INTRODUCTION
Interesting information about complex physical, biological and environmental systems is communicated to the external world by weak and low-frequency magnetic signals reflecting the evolution of the weak magnetization carried by these systems. Similarly, in certain quantum mechanical ensembles such as spin polarized ultra-cold neutrons or atoms, the interesting physical information is superimposed on top of their weak magnetization, which must be strictly controlled in order to observe the desired signal. There exist well established methods to measure DC magnetic fields down to a few pT. Traditional methods to create low magnetic field environments are passive shieldings made of high permeability ferromagnetic alloys or superconducting enclosures. Examples of such systems can be found e.g. in Refs. 1-4. Theoretically, there exists another approach to provide a "magnetic vacuum" in a given volume and frequency range. For a charge-and current-free space, in the quasi-static regime, this is guaranteed by solutions of the Laplace's equation arising from Maxwell's laws of classical electromagnetism. 5 Knowing field vectors at sufficiently numerous points in the volume it is possible to judge the field distribution in a given volume to the required accuracy. This, in turn, allows calculating and applying electric currents in a system of correction coils generating (in the same volume) a field distribution a with the same magnitude and sign opposite to the measured field. The correction signal must appear with a delay negligible for the frequency range in question. Active field compensation utilizing threedimensional Helmholtz-like coil systems were applied e.g. in MiniCLEAN, 6 LIPSION 7 and nEDM 8 experiments.
In practical solutions, especially for a large control volume, one is faced with limitations, the most important being the maximal number of independent correction coils included in the feedback loop. Each coil needs a precise and quickly reacting current source with control electronics which dissipates power (heat) which may cause additional problems in the particular application. The main goal of this study is to estimate the performance of an idealized system with as few degrees of freedom as possible. Such a case will be used in future as a reference for practical solutions.
In this article we propose a novel design for an active compensation coil system based on the field decomposition in the Vector Spherical Harmonics (VSH) described by R. Barrera et al. in Ref. 9, together with a stream function description of currents to obtain an orthogonal system of coils.
We first present the VSH formalism consisting of (i) a description of the field in a given volume created by an external source, (ii) evaluation of the current density distribution generating such a field and (iii) a current density discretization procedure leading to spherical coils. Then we characterize the performance of the spherical coil system in response to various field disturbances and compare it to the performance of a 6 coil system, used by the nEDM experiment at PSI.

II. MAGNETIC FIELD DECOMPOSITION IN TERMS OF VECTOR SPHERICAL HARMONICS
The general situation is illustrated in Fig. 1. We divide the whole space into three distinct regions. The central volume called the experimental volume V exp is assumed to be free of magnetic field sources (electric currents and magnetic material) and is supposed to be free of magnetic fields generated by external sources when the compensation system is active. The experimental volume V exp is surrounded by a shell V C , to which electrical currents are constrained. Outside of V C , in the volume V D , magnetic field disturbances can occur. S C is the surface (see section III), on which the solutions for coil turns will be obtained. We assume that all materials within V exp and nonmagnetic (µ r = 1) which allows us to define the scalar magnetic potential φ( r), fulfilling the Laplace's equation.
Electrical currents in the shell V C generate a magnetic field according to Ampere's law: FIG. 1. Definition of geometrical components. V exp is the volume where the field is supposed to be compensated. V C is the volume where compensating coils are supposed to be located. S C is the surface (see section III), on which the solutions for coil turns will be obtained. V D is the volume where the disturbing field originates from.
where A is a vector potential and J is a current density.

A. Spherical harmonic decomposition of the scalar potential
In the spherical coordinate system, the scalar potential φ( r) -the solution of the Laplace equation -can be expressed as a series of spherical harmonics: where the set of {φ lm } describes uniquely the scalar potential function. For small field disturbances from external sources one expects a fast convergence of expression (2) near the origin. Only the first few terms with small l values will contribute to the field description. An additional advantage of this form results from the orthogonality of spherical harmonics. If the external field changes only a limited number of the terms in Eq. (2) are affected. The reacting compensation system will have to drive currents only to the required corresponding coils without affecting others.

B. Current density distribution
In order to compensate a field disturbance in the controlled volume V exp one has to drive currents in the shell V C such that the resulting magnetic field distribution in V exp has the same amplitude as the disturbance with opposing sign. This problem is often referenced as the "Inverse source problem" 10 and has been subject of many studies. Here we present an analytical solution for a system exhibiting spherical symmetry. To accomplish this we begin with a synopsis of Vector Spherical Harmonics (VSH) essentials. A more detailed discussion can be found in Ref. 9. Vector Spherical Harmonics are related to the scalar spherical harmonics in the following way: where:r is a versor in radial direction, (r, ϑ, ϕ) are spherical coordinates. These equations offer is a complete and orthogonal basis set, meaning that any vector function can be expressed as a series of VSH: where {f Y lm , f Φ lm , f Ψ lm } is a set of expansion factors. A significant advantage of VSH is their well known and simple behavior under action of differential operators. 9 Any vector potential A(r, θ, ϕ) can be expressed as a series of VSH: Making use of the relations for the divergence and gradient of the Y part of VSH expansion 9 we express the magnetic field in the current free experimental volume V exp using scalar potential (Eq. (2)): and in the current shell V C : The current density J depends only on Φ lm (components J Y lm and J Ψ lm of the current density expansion into VSH series vanish for all l and m meaning that currents flow only in directionsˆ θ andˆ ϕ). The expression for the magnetic field generated by currents in the shell V C will simplify significantly, similarly for the equations for currents. This problem was already solved in Ref. 9 This means that a current J ∝ Φ lm flowing on a sphere generates a field described by a scalar potential φ ∝ r l Y lm . Hence, the fields generated by different J ∝ Φ lm will be orthogonal inside the sphere.

C. General current densities
If we relax the constraints of an infinitely thin current density (coefficients J Y lm and J Ψ lm do not vanish anymore) we obtain the general expression: Making use of the continuity relation for current ∇ · J = 0, and Gauss's law for magnetism ∇ · B = 0 as well as exploiting the relations for the divergence of VSH, we arrive at: This means, that the Y and Ψ parts of the series expansion of A and J (see Eqns. (6) and (7)) are dependent and only one of them needs to be considered, e.g. Y. From the curl of VSH (see Ref. 9) we get: where B Φ lm (r) is the function describing the expansion of the magnetic field into the Φ component of VSH. This means that for x ∈ V exp , ( J = 0) also In conclusion to the current density expansion in V C , the terms proportional to Y lm (θ, ϕ) and Ψ lm (θ, ϕ) do not contribute to the field generated in the controlled volume V exp . This means that only currents flowing on the surface of the sphere generate field inside that sphere.

III. COILS -DISCRETIZATION OF THE CURRENT DENSITY DISTRIBUTION
In a practical application the required current density distribution J(x, y, z), where (x, y, z) ∈ V C , will be obtained using a coil system. We assume that the coils consist of thin wire turns distributed AIP Advances 7, 035216 (2017) on a sphere, so that volume V C reduces to a surface. Since all wire turns in a coil are connected in series, the only way to construct a given current density is by proper shaping and distribution of wire turns. To perform this we use the stream function approach from Ref. 11. It is based on the fact that for a static current density function J(x, y, z) the divergence vanishes If currents are constrained to a surface defined by a normal vectorˆ n, then the current density can be written as: where ψ(x, y, z) is the so called stream function defined on the surface containing the currents. This function must be differentiable on that surface. In order to fulfill the relation (16) the stream function must be constant on the boundary of the surface. The currents are allowed to flow along the isolines of the stream function. We can identify the isolines with the coil wire turns. Comparing Eqn. (17) with the Φ part of the VSH series in Eqn. (5) we realize that wire turns of the spherical coils follow the isolines of spherical harmonics.
The necessary number of the coil windings depends on the required field reproduction accuracy. Fig. 2 shows the dependence of the average relative reproduction accuracy |B reproduced |/|B theoretical | as a function of the number of wire turns. The more turns a coil has, the more accurate the field created by it reproduces the field of a continous current distribution. Three l = 1 coils (Fig. 3), the so called cos-theta coils, 12 identical to those obtained by J. E. Everett and J. E. Osemeikhian, 13 produce a uniform magnetic field. Field components of higher order are produced with coils with l > 1. Figures 4 and 5 present the wire distributions for l = 2 and l = 3, respectively.

IV. SPHERICAL AND RECTANGULAR COIL SYSTEMS FOR A COMPARISON
In order to compare the expected performance of the spherical coil system with more 'standard' solutions we describe first two rectangular coil systems. The first one is used in the neutron Electric Dipole Moment spectrometer installed at the Paul Scherrer Institute, Villigen, Switzerland 1,8 ). The system consists of three pairs of rectangular coils as shown in blue in Figure 6. The coils (dimensions -see Table I) have independent power supplies providing 6 degrees of freedom for the compensation system. The second 'standard' used for comparison is based on the Merritt 4-coil system 14 -a candidate for the second generation of nEDM experiment at PSI shown in Fig. 7. It consists of 12 coils, providing 12 degrees of freedom. The dimensions of the coils are given in Table II.
The necessary reaction currents corresponding to the external field changes are calculated using the readouts from a number of vector fluxgate sensors. The algorithm utilizes the fact that the magnetic field strength generated by a coil is proportional to the current driven through it. This leads to a system of linear equations relating the measured field components with the coil currents 3 · · · G z1 m · · · · · · · · · · · · · · · G zn 1 G zn 2 G zn 3 · · · G zn m I 1 The G ij k coefficients are determined by recording the fluxgate magnetometer measurement values while changing the current in one coil at a time. The minimum number of 3-dimensional fluxgate sensors necessary to uniquely establish the coil currents depends on the number of coils in the system. Only two sensors are required in the case of the 6-coil rectangular system. However, such a system would compensate the field changes only in two points (location of fluxgate sensors) ignoring the field changes in other parts of the controlled volume. Moreover, the uncertainties of the measured fields propagate to the calculated correction currents. In order to reduce these drawbacks the system utilizes more fluxgate sensors (12 in the example below) leading to an overdetermined set of equations. The pseudoinverse matrix G ☞1 is calculated using the Singular Value Decomposition method 15 to calculate the optimal currents.  To calculate the coil currents in all systems we adopted corresponding algorithms from the GNU Scientific Library. 16 Magnetic field perturbations generated by a given current loop were calculated by integrating Biot-Savart's law.
FIG. 8. Average ∆ 1 as a function of the radius of the spherical control volume with feedback sensors located on that surface for spherical and rectangular coil systems. and △ mark results for spherical coil systems -8 coil and 15 coil, respectively. denote 6 coil rectangular system, ♦ -12 coil one. Full symbols denote the disturbance current loop located at (10,15,20)  In order to compare the performance of coil systems we assume that the field measurement sensors are located at randomly chosen positions on a sphere with 2 m long radius (centered in the coil system center). The chosen radius for the spherical coils is r = 2.93 m and the number of wire turns in each spherical coil is equal to 100; we assume that the thickness of the wires has negligible influence on the magnetic fields. In this work we use two functions of (x, y, z) describing the performance of the coil systems called ∆ 1 and ∆ 2 and defined in Eqns. 21 and 22.
where B per denotes the field generated by an external perturbation source, B cor the correction calculated by the compensation system, and B per max is the maximum value of the magnetic field perturbation inside the experimental volume.
∆ 1 reflects the absolute value of the (uncompensated) field remnant relative to the local value of the perturbation |B per |, while ∆ 2 relates this remnant to the maximum value of perturbation in the experimental volume. Thus ∆ 1 reflects the distribution of the inverse of the shielding factor in the experimental volume while ∆ 2 describes the size and distribution of the remnant (uncompensated) field in that volume. Both figures-of-merit, ∆ 1 and ∆ 2 , are equivalent once B per is known.
Since the expected external field disturbances at the nEDM experiment in PSI are of dipolar type, we model them with a current loop at a fixed distance, orientation current and diameter. This These maps allow to conclude that, for the spherical system, the magnetic field is fairly well compensated in almost the entire experimental volume, whereas in the case of rectangular coils, the compensation performance expressed in terms of ∆ 1 and ∆ 2 is at least an order of magnitude worse. Regions with compensation errors ∆ 1 less than a given value are much bigger in the case of the spherical coils than for rectangular coils. Adding the third order coils (l = 3) in the spherical systems still improves the shielding factor by almost an order of magnitude.
In order to draw a more detailed conclusion the radius of the spherical experimental volume V exp was varied. The 10 sensors were located at random points on the surface of V exp . The averaged (over the controlled volume) relative differences ∆ 1 and ∆ 2 are presented in Figs. 8 and 9, respectively. This calculation shows that the spherical coil system performs on average 1 to 5 orders of magnitude better than the rectangular version, depending on the size of the experimental volume V exp , independently of the figure-of-merit choice. Figure 10 shows the compensation performance dependence, as a function of the distance of the disturbance field source. For each of the points on this plot, 100 random positions and orientations of the sources at a given distance were averaged to obtain the most universal picture.
The performance of the rectangular coil system is almost independent on the distance from the disturbance source. For spherical systems, the average ∆ 1 decays with the distance to the source as r −(l max ) . The only deviation from this law can be observed for the 15 coil system. For big distances, the average ∆ 1 flattens. This effect can be attributed to errors introduced by using the discrete approximation instead of perfect continuous current distribution. By increasing the coil wire density (from 100 to 400 wires per coil), the lowest value of ∆ 1 gets 10 times smaller than before. For the eight coil system we do not see this effect because the error resulting from the cut-off of the Vector Spherical Harmonic expansion is bigger than the one from the current discretization. The spherical coil system with 15 coils performs 2 to 5 orders of magnitude better than the rectangular system with six coils, depending of the size of the experimental volume relative to the size of the coil system and distance from the disturbance source.
This gain in performance is not caused just by the larger number of degrees of freedom. Adding 6 rectangular coils, resulting in a 12 coil system, does not improve the performance significantly, as can be seen by comparing plots in Figs. 11 and 12. The only noticeable improvement is visible in the larger sizes of controlled volume, where the improved uniformity of the field created by the 12-coil system matters.

V. CONCLUDING REMARKS
In this paper, we described a method for designing active field compensation coil systems by solving the inverse source problem of magnetostatics using vector spherical harmonics (VSH). The perturbation field to be compensated is decomposed in the volume inside the coils in the spherical harmonics and the correction currents are calculated for the individual coils producing the orthogonal correcting field components. We show that only the Φ part of the VSH expansion of currents on the current shell generates fields inside the sphere. Based on the above result, two systems consisting of 8 and later of 15 spherical coils were selected for further evaluation.
The operation of such a system was simulated using disturbances created by a current loop, and was compared to the already working system consisting of 6 coils and a proposed system with 12 rectangular coils. We conclude that a 8-coil spherical system with only two additional degrees of freedom (8 instead of 6) would be a significant improvement compared to the currently used system. Our simulations show that the performance of the proposed magnetic field compensation system expressed as an effective shielding factor and as the (uncompensated) field remnant, might be orders of magnitude better than that for the rectangular coils systems.
Our calculation shows that the optimization of the wire density (current discretization) should be performed only after the range of expected positions of dipolar disturbances and the required compensation accuracy are defined.
In a further study we plan to apply the VSH formalism for different shapes (e.g. cuboid, cylinder, truncated cuboid), which are suitable for a practical realization of a compensation systems to be applied in the neutron Electric Dipole Moment experiment.

ACKNOWLEDGMENTS
This work has been supported by the Swiss National Science foundation under number 200020 149211 and by The National Science Centre, Poland, under the grant No. UMO-2015/18/M/ST2/00056. We acknowledge support by the Foundation for Polish Science -MPD program, co-financed by the European Union within the European Regional Development Fund. The LPC Caen acknowledges the support of the French Agence Nationale de la Recherche (ANR) under reference ANR-09BLAN-0046.