Theory and development of biplanar active shim coils for a permanent NMR analyzer

The general expression for magnetic induction in the z axis direction is derived from magnetic scalar potential, and magnetic induction for biplanar shimming coils (BSCs) is also derived from magnetic vector potentials and Green functions, which simultaneously include Sin and Cos harmonic fields. The relationship between these expressions is discussed, and we show they are partially consistent. Magnetic induction generated Sin and Cos stream functions, which are presented and discussed, and we conclude that the type of stream function determines the type of harmonic field, and that BSCs can not only generate specific harmonic fields directly using Cos stream function, but also generate the rest of the harmonic fields through some specific operations. The detailed design process is presented in the form of a diagram. Subsequently, nine BSCs were calculated using the proposed method and applied to a low field NMR relaxation analyzer. The magnetic field homogeneity after shimming increases significantly, which verifies its practical value.


Introduction
Low field NMR (LF-NMR) relaxation analysis technics has been widely employed in oil exploration, food safety, life sciences, high molecular material, and mineral engineering since it was developed in the 1980s [1][2][3][4][5]. LF-NMR has advantages in running costs, and it is non-invasive, dependable, and highly secure, and it will be expanded to new applications as portability and miniaturization develops. However, it also has disadvantages relative to high field NMR, such as low resolution and low sensitivity due to the magnetic field inhomogeneity, which has shown LF-NMR to be insufficient in detecting small composition content, and it has a short relaxation time.
The primary requirement of any NMR system is that its main magnetic field must be highly uniform. LF-NMR usually adopts rare earth permanent magnets to generate the main magnetic field, which generally has significantly lower homogeneity than that generated by superconducting magnets. Two methods are used to improve field homogeneity: passive and active shimming. Passive shimming is the basic method, with active shimming adjusting the field further to attain improved homogeneity. Passive shimming must be carried out in advance because the initial coil current using just active shimming may be so high that the generated heat will affect the magnet temperature stability. This paper investigates active shim coil design to improve homogeneity based on passive shimming.
Target field is a well-known design method proposed by Turner [6] in 1986. This method has no ill-conditioned problem because Fourier transforms have unique inverses, and it has been used to design cylindrical surface shim coils for superconducting magnets for magnetic resonance imaging (MRI). Forbes and Crozier [7][8][9] advanced the TF method, considering the finite coil length using a Fourier series technique. Wen Tao Liu [10][11][12] proposed an alternative current density expansion, rather than Tikhonov regularization, and Xiao Fei You [13] combined this method with Forbes and Crozier's. Ge Li Hu et al. [14,15] employed L1 norm regularization to replace the original ill-conditioned problem with a nearly well-conditioned problem. There has been considerable TF research, and application of TF in engineering has expanded rapidly. However, shim coils based on TF focus on superconducting magnetic resonance spectroscopy (SMRS) to gain high homogeneity for most research, which generally has a small region of interest (ROI), approximately F × H = 5 mm × 20 mm. This paper proposes a biplanar shim coils (BSCs) approach for LF-NMR to improve homogeneity to 0.5 ppm over a relative large ROI (F × H = 25.4 mm × 40 mm). Most literature has only simulated results to demonstrate the efficiency of their proposed approach, whereas this paper uses physical experiments to validate the proposed approach.
The general expression for magnetic induction from a magnetic scalar potential and the specific magnetic induction for BSCs are derived, and their relationship discussed. We present a more general methodology, valid for spherical harmonics that includes Sin and Cos functions simultaneously. We demonstrate that the BSCs not only directly generate some specific harmonic fields, but also generate all other harmonic fields through some specific operations. The detailed design process is shown diagrammatically, and nine BSC pairs are designed based on the indicated process, then installed in a biplanar permanent MRI. An experiment was performed to verify the shimming effect, with very significant results, which indicates that the theory and development approach is viable.

General formula
The main field, B 0 , is usually in the z axis direction, which is perpendicular to the radio frequency (RF) direction. In diameter spherical volume (DSV), consider a small deviation of B 0 , B' ( B 0 , where the total field is Since RF is perpendicular to B 0 , only B' in the z axis direction affects the NMR signal, so in the DSV through which no current passes for a static magnetic field, so that Δ × B = 0 from Maxwell's equations. It then follows from a vector identity that where C and μ 0 are the magnetic scalar potential and magnetic conductivity in a vacuum, respectively. Maxwell's equations require Δ Á B = 0, which means C satisfies the Laplace equation in air in DSV, i.e., Eq (4) can be solved in spherical polar coordinates rather than Cartesian coordinates, using the usual transformation relations between the coordinate systems, which are presented here explicitly for clarify (Fig 1), and the general solution in spherical coordinates can be expressed as Cðr; y; 0Þ ¼ S 1 n¼0 S n m¼0 r n P m n ðcos yÞðA m n cos m0 þ B m n sin m0Þ; where n and m are the degree and order of C, respectively; P m n is an associated Legendre function; and A n m , B m n are coefficients independent of r,θ, and ϕ. When n = 0, then m = 0 also, from Eq (6), and C is constant.
From the gradient transformation, Substituting Eqs (6) into (7), B z ¼ m 0 fS 1 n¼1 S n m¼0 r nÀ 1 ½ð2n þ 1Þcos yP m n ðcos yÞ À ðn À m þ 1ÞP m nþ1 ðcos yÞ Â ðA m n cos m0 þ B m n sin m0Þg; which can be simplified using the properties of Legendre functions as Since m n for P m n ðxÞ for Legendre functions, then m n −1 for P m nÀ 1 ðcos yÞ in Eq (8), so and for convenience, let α = n − 1 and β = m, then where α and β here are the degree and order of B z . Eq (10) is the general expression derived from C, and is compared to the BSC derived expression in the following section. This indicates that B z is composed of a series of harmonic fields, which is consistent with Anderson's theory [6]. In the ideal case, Eq (10) should contain only the uniform field, B z0 , and the constant A 0 1 with no higher harmonics present when α = β = 0. For given α 6 ¼ 0 and β 6 ¼ 0, r a P b a ðcos yÞðcos b0jsinb0Þ determines the distribution of B z 's higher harmonic. The purpose of shim coils is to eliminate higher harmonics and keep the field uniform. Hence, the target field should be chosen from Table 1.
Biplanar shim coil generated magnetic field from scalar potential According to Maxwell's equations, Δ Á B = 0. Therefore, there is a magnetic vector potential, where J is the current density, and Δ 2 A = −μ 0 J. The general solution was obtained in terms of the Green function 1⁄|r − r'|, and can be expressed in differential form as and the differential form of the magnetic field can be expressed as However, J is planar rather than a volume distribution. Hence, Eq (11a) becomes where ds = ρdρdϕ. Let G = 1⁄|r − r'| where G is the Green function, and define the gradient where G 0 x and G 0 y are the components of G' in the X and Y directions, respectively; and J x and J y are the components of J in the X and Y directions, respectively.
In spherical coordinates, Eq (12) can be expressed as where G 0 r 0 ,G 0 y 0 , and G 0 φ are the components of G' in spherical coordinates; and J ρ and J φ are the components of J in polar coordinates. Since In the interval r a < r' for DSV, the Green function can be expanded [16] (14) into (13), From the recurrence property of associated Legendre functions [17], Let and then Eq (17) becomes From Fig 2, r' = (a + ρ) 1/2 , sinθ' = ρ/r' and cosθ' = a/r'. Therefore, C αβ , C r ab , and C φ ab are all functions of ρ for given α and β.
Eq (18) may be integrated over the whole BSC, where ðC φ ab J φ sinbφ À C r ab J r cosbφÞC ab rdφdr: Eqs (10) and (19) have similar forms, both have the term r a P b a ðcos yÞ cosb0 sinb0 ! ; which indicates they should have the same harmonic fields, i.e., sinusoidal harmonic fields. Therefore, BSC-generated harmonic fields could eliminate the effects of high order harmonic fields from the main field (Eq (10)) by selecting an appropriate current distribution. The effect of current distribution is embodied in A 0b aþ1 and B 0b aþ1 , which are discussed below. The continuous current density, J(ρ, φ), can be divided into a tangential component, J φ , and radial component, J ρ . J(ρ, φ) satisfies continuous flow, i.e., Δ J = 0. Let us introduce the stream function S(ρ, φ), such that Then the stream function may be expressed as where l and k are the degree and order of coils, respectively; "±" represents the upper and lower planes, respectively; and the expansion basis, s q (ρ), may be randomly chosen, e.g. s q (ρ) = sinqπρ⁄ρ max . Eq (21) is a Cos stream function. The current density's tangential component and radial component may be derived from Eqs (20) and (21), : where U q is the current coefficient, j φ,q (ρ) = −@s q (ρ)⁄@ρ, and j ρ,q (ρ) = −ks q (ρ)⁄ρ. Therefore, substituting Eqs (22) into (19), þ½C rþ ab þ ðÀ 1Þ lþk C rÀ ab j r;q sinkφsinbφg C ab rdφdr; ð23Þ ab , and C rÀ ab ¼ ðÀ 1Þ aþb C rþ ab , because cos(π − θ') = −cosθ', sin(π − θ') = sinθ' and P b a ðÀ cosy 0 Þ ¼ P b a ðcosy 0 Þ, Eq (23) becomes where δ kβ is the Kronecker function. For A 0 b aþ1 6 ¼ 0, two conditions must be satisfied, b ¼ k and BSCs can only generate the Cos harmonic fields from Eq (28), i.e., r a P b a ðcos yÞcosb0: Hence, Eq (28) is called the Cos B z function, which is determined by the form of the stream function discussed further below.
If we change the stream function Eq (21) into the other form, S AE ðr; φÞ ¼ ðAE1Þ lþk sin kφ S Q q¼1 U q s q ðrÞ; 0 r r max ; ð30Þ which still satisfies continuous flow. Eq (30) is called a Sin stream function, and the corresponding current density may be expressed as where U q is still the current coefficient, j 0 φ;q ðrÞ ¼ j φ;q ðrÞ ¼ À @s q ðrÞ=@r, and j 0 r;q ðrÞ ¼ À j r;q ðrÞ ¼ ks q ðrÞ=r.
Substituting Eqs (31) into (19), and following the derivation above, When B 0b aþ1 6 ¼ 0, Eq (27) must be satisfied, and similarly Eqs (26) and (33). Therefore, Eq (19) becomes which is called a SinB z function. which can only generate Sin harmonic fields, i.e., Since Eqs (28) and (34) both meet the conditions of Eq (27), when β = k, the conclusion that the stream function type determines the B z function type is evident by comparing Eqs (27) and (34) with Eqs (21) and (30). In other words, Sin stream functions can only generate Sin harmonic fields and the Cos stream functions can only generate Cos harmonic fields. If we want to generate corresponding higher harmonic fields while retaining uniform B z0 as mentioned above, the target harmonic field type must be first determined from Table 1. Eq (30) shows that the stream and harmonic functions must choose Cos function when β = k = 0. Eq (28) shows the results as [12]. However, the methodology in the current paper is more general, since it is valid for spherical harmonics that include Sin and Cos functions simultaneously, whereas [12] only presented Cos harmonic functions, with the generation of Sin harmonic functions only represented by rotating at specific angles. Although [12] conclusions were correct, they lacked a strict proof.
The condition α!β should also be considered along with the two conditions of Eq (27), ( where i = 1,2,3,. . .N; and N represents the number of generated harmonic fields. Eq (36) shows that the shim coil for a given degree l and order k can generate an infinite number of harmonic fields depending on the value of i. Hence, the fields generated by a coil with degree l and order k may be expressed as and from Eq (26), where Eq (38) may be expressed in matrix form as q) can be calculated from Eq (39), which is constant for the corresponding values of i and q.V k a is the coefficient matrix of harmonic fields.
The restricting condition, will make the coil only generate one specific type of harmonic field with l = α(i) while the other N-1 types will not be generated. Values of degree l and order k may be given as targets in advance. The nonzero term's value of V k l may be arbitrarily chosen except zero, because V k l is proportional to the current coefficient matrix, U, when D is constant, and the current coefficient matrix U can be solved. U's magnification only changes the current magnitude, not the distribution. For example, when l = 1, k = 0, then β = 0 and α = 2i − 1 from Eq (36), and there should be N types of harmonic fields generated by the coil. However, from Eqs (41) and (29), only the rP 0 1 ðcos yÞ harmonic field is built, i.e., the Z direction gradient field. S(ρ, φ) may be obtained by substituting the solved matrix U into Eq (21). However, since it is continuous, it cannot be used directly for engineering applications. Therefore, S(ρ, φ) must be processed discretely before application.
S(ρ, φ) may be represented by some current-carrying coils, with the coil current, determined by where I max and I min are the maximum and minimum values of S(ρ, φ), respectively; and M is the order of S(ρ, φ), which is set in advance. Therefore, the discrete S(ρ, φ) can be expressed as where j = 0, 1, 2, 3Á Á ÁM − 1.

Design diagram and results for shim coils
5. The discrete S(ρ, φ) was calculated from Eqs (42) and (43), and is called a Z4 coil.
Following the process of Fig 3, nine pairs of BSCs were calculated that generate X, Y, Z, Z 2 , XZ, YZ, X 2 -Y 2 , Z3 and Z4 harmonic fields (Table 1), as shown in Fig 4. The coils were used for shimming the main magnetic field of the high-performance relaxation analyzer. The gap between magnetic poles is 120 mm, i.e., a = 60mm, and the target field volume is F × H = 25.4 mm × 40 mm. Therefore, the maximum distance from the zero point is less than 60 mm, which meets the requirement of Eq (14).   Table 1.
https://doi.org/10.1371/journal.pone.0181552.g004 by clockwise and anti-clockwise coils, respectively. Fig 5(a) shows the engineering drawing for the shim coils. The pair of printed circuit boards (PCBs) shown in Fig 5(b) and 5(c) are multilayer boards 2 mm thick. Each PCB for the nine shim coils was installed on the two poles of a permanent biplanar MRI in symmetry. A pair of coils generating the same harmonic in two PCBs were connected in series by wires. Thus, there were nine independent groups of constant current sources to power the nine BSC pairs, respectively.

Experiments to test shimming effects
An experiment was conducted to test the effects of the shim coils. The magnetic intensity of the biplanar permanent MRI = 0.5T, hence resonance frequency = 21.3 MHz for pure water. The adjustment order of shim coil current was Z, X, Y, XZ, YZ, Z 2 , X 2 -Y 2 , Z 3 , and Z 4 . The effect of shimming can be measured by comparing the free induced decay (FID) integral area before and after shimming. Once the current in the first BSC pair was adjusted to achieve maximum integral area, then the next BSC pair was adjusted, and so on to through all BSCs. BSC maximum current was 310 mA, and total thermal power < 10 W. Shim coil electric parameters are shown in Table 2, where the resistance terms represent the total resistance of a pair coil in series. Fig 6 shows the test equipment.
Magnetic dephasing caused by field inhomogeneity produces additional suppression of the signal. Hence, damping the magnitude of transverse magnetization is described by replacing  the decay time, T 2 with a modified time, T Ã 2 . Although qualitatively they have the same shape, the new envelope typically decays much faster in time. Decay time T Ã 2 represents a combination of external field induced T 0 2 and thermodynamic T 2 effects, 1 where T 0 2 is both machine and sample dependent [18]. Since the sample used in these experiments was pure water, and the sample's effect on T 0 2 could be ignored, and hence, T 0 2 only depends on inhomogeneities of the main field in this case. Generally, T 0 2 dependence is used in FID measurements to determine if main field homogeneity is adequate. This method was also used in the experiment, and the Fourier transform of FID (data in S1 Table) are shown in Fig 7. The shimming effect is clearly evident. As the FID integral area increases from 28,900 to 5,892,900, an approximately 20-fold increase, the full width at half maximum (FWHM) of the frequency spectrum decreased from 17.90 to 0.81 ppm, and the corresponding full width at 1/10 maximum (FWTM) decreased from 79.62 to 4.87 ppm.
Since thermal deformation of the shim coils will decrease the shimming effect, stability of the shimming effect was also considered, and estimated by measuring FWHM and FWTM every 10 minutes over 24 hours. The average change rate was approximately 0.50%, with a maximum rate of approximately 0.85%.

Discussion and conclusion
The general expression for magnetic induction was derived from magnetic scalar potential and indicates that the field may be expanded as a series of harmonic fields, which is the basis of target field theory. Shimming aims to eliminate the effects of high order harmonic fields through reversed harmonic fields induced by the shim coils.
The harmonic fields contain Sin and Cos harmonic functions, induced by the corresponding Sin and Cos stream functions. The BSCs not only directly generated specific harmonic fields, but also generated other harmonic fields through specific operations.
The design process was shown diagrammatically, and nine BSCs were designed and fabricated based on the theory developed here, then installed in a biplanar permanent MRI. The shimming effect was very significant, which indicates that the theory and development approach is viable.
Supporting information S1

Author Contributions
Data curation: Tian Xia.