Methods: Of Stream Functions and Thin Wires: An Intuitive Approach to Gradient Coil Design

The design of gradient coils is sometimes perceived as complex and counterintuitive. However, a current density is connected to a stream function in fact by a simple relation. Here we present an intuitive open source code collection to derive stream functions from current densities on simple surface geometries. Discrete thin wires, oriented orthogonally to the main magnetic field direction are used to describe a surface current density. An inverse problem is solved and stream functions are derived to find coil designs in the current and stream function domains. The flexibility of the design method is demonstrated by deriving gradient coil designs on several different surface topologies. This collection is primarily intended for teaching, as well as for demonstrating all gradient coil design steps with openly available software tools.


INTRODUCTION
Gradient coils are the key component to enable imaging using nuclear magnetic resonance (NMR). Spatially varying magnetic fields are switched in a time dependent manner to achieve spatial encoding of the NMR signals. Different approaches to gradient coil design have been proposed over the past decades to derive very sophisticated designs compared to initially implemented Maxwell and Golay type coils [1,2].
Most approaches to gradient coil design today are based on boundary element methods [3,4]. These approaches are based on representing the actual electrical current density by its stream function discretized on a surface mesh. Due to the counterintuitive approximation of a surface current density by equidistant iso-contour lines of the stream function the coil design problem is oftentimes perceived as complex. However, the relation between the current density and stream function is oftentimes perceived as complex. For educational purposes the level of complexity of the boundary element method can be reduced to understand and design simple gradient coils. The coil design problem can be derived intuitively with an intermediate step which represents the current density. In addition it can be simplified by reducing its dimensionality if the current carrying surface is oriented in parallel to the main magnetic field. This is done by neglecting currents which flow in parallel to B 0 and therefore do not generate a B z component relevant for encoding. Therefore we propose the most simple framework for coil design based on thin wires oriented orthogonally to the main magnetic field. This paper aims at providing an intuitive access to gradient coil design.
Additionally, Stream functions are derived from an intermediate representation of current densities and by a direct solution. The straight-forward, yet powerful approach presented here is suitable for simple, regular surface geometries such as cylinders or planes.
Most coil designs are derived by deploying commercially available closed source software, e.g. [5,6]. To the knowledge of the authors there is no code openly available which includes the surface mesh definition for gradient coil design, as well. One repository available on GitHub by Bringout et al. [7]. needs a predefined triangular surface mesh. Here, the mesh definition is included within the code that accompanies this paper and the complete sources are available on GitHub.

THEORY AND METHODS
All scripts to demonstrate this method were designed to run with GNU Octave and MatLab (The Mathworks, Nattick, MA, United States). Images in this manuscript were plotted with GNU Octave on a laptop computer.
Gradient coils for spatial encoding are usually operated with frequencies below 10kHz, therefore the coil design may be treated as a magneto-static problem. The relation between the magnetic field of the gradient coil, B(r), and a free current density, J(r), is given by Ampere's law: where μ 0 is the permeability of free space. Conventional magnetic resonance imaging (MRI) usually deploys a main magnetic field which is highly uniform, unidirectional and very strong to achieve a sufficient spin polarization. Due to its historic origin from NMR, direction of the main field B 0 is by convention chosen along the z-axis. The local Larmor precession frequency only depends on the absolute value of the superposition of the main magnetic field and the field induced by the gradient coil. A Taylor expansion of this magnitude of the superposition field yields whereẑ is a unit vector along z. For high field systems the magnitude of the main magnetic field, B 0 , exceeds the magnitude of the gradient encoding fields typically about two orders of magnitude. Therefore, the expansion in Eq. 2 can typically be terminated after the first order term, corresponding to the z-component B z (r), and the transverse field components B xy (r) can usually be neglected.

Thin-Wire Approximation
For simple, regular surface geometries such as planar or cylindrical, a current density, J(r), may be approximated by discrete thin wire segments. As only B z is considered relevant for spatial encoding and due to the curl operation as in Eq. 1 between the current and the generated magnetic field, the orientation of each current carrying element, m, is chosen to be orthogonal to the z-axis.
Wires running parallel to B 0 produce no B z component and can be used in this model to feed the orthogonal current segments. All feeding wires are allowed to overlap in space, as wire thickness is ignored in the basic thin-wire approximation. Each orthogonal thin-wire segment m with a current I m contributes to each target point n, which results in a spatially dependent magnetic field B z,n , at point n, depicted in Figure 1A.
The relation between the current in each thin-wire segment and the field of each point in space is described by the Biot-Savart law: Combining these relations defines the sensitivity matrix S nm . By representing the currents in each thin wire segment as a vector, I m , and extracting only the B z component of the field vector, simple algebraic expressions become possible. The resulting gradient field is defined by the column vector, B z,n , according to B z,n S nm · I m . (4) A corresponding current distribution for a given target field B target z may be calculated by deploying the pseudoinverse S + nm of S nm . However, depending on the number of current elements m and number of target points n, the resulting system of linear equations may be under-determined. In addition, this inverse problem is ill-conditioned, as most inverse problems in magnetostatics. A direct inversion of such problems typically leads to impractical solutions such as depicted in Figure 2A.

Power Optimization
A Tikhonov regularization may be used to find solutions with smaller norms, which penalizes high opposing currents according to the following minimization problem Here, · 2 is the Euclidean norm and Γ describes the Tikhonov matrix which is in this case chosen as the identity matrix multiplied by a regularization parameter λ: The identity matrix in the Tikhonov regularization penalizes the sum of squares of the currents in each of the elements. This is proportional to the power dissipated by all elements and therefore acts to limit the total required electrical power. From a physics point of view the coefficient λ is related to the resistance of the wire segments. An appropriate choice of the weighting factor is balancing efficiency against accuracy of the resulting field. This free design parameter has to be chosen by the coil designer. Different approaches, e.g. the L-curve method [8] may be used to choose practical values. An exemplary regularized current distribution is shown in Figure 2B.

The Stream Function Representation and Discretization
Discrete wire layouts of gradient coil designs are usually derived from a stream function representation of surface currents. The stream function is then used to plot equidistant iso-contours, so Frontiers in Physics | www.frontiersin.org September 2021 | Volume 9 | Article 699468  Frontiers in Physics | www.frontiersin.org September 2021 | Volume 9 | Article 699468 3 called streamlines, representing trajectories of equally spaced steady flow of current. The number of streamlines is chosen according to engineering constraints, e.g. wire thickness, minimum wire spacing or maximum inductance.
The stream function is related to the current density by its integral. An integration of the previously regularized current density is depicted in Figure 2C. In matrix form, the integration may be directly acquired from the current density by calculating the cumulative sum along the z-direction according to Here, it should be noted that the depicted discrete coil layout may not be realized by closed loops, leading to a gradient coil requiring many current terminals on each side. The example shown in Figure 2C would require current terminals and sinks of twice the number of iso-contours on each side, which is impractical.

Realizing Closed Loops
As depicted in Figure 2C, even the regularized current distribution may not be realizable in a straight-forward manner due to a high number of feeding ports. An additional constraint, enforcing the sum of currents along the z-axis to be zero along the surface boundaries, ensures for a design with closed loops. This constraint is given by In this equation the index m is replaced by two indices, corresponding to the position of represented segments in the 2D matrix. Effectively, this is equivalent to Kirchhoff's current law, which states the principle of the electric charge conservation. If this condition is fulfilled, for each current flowing along the direction orthogonal to z there is an opposing current, or currents with the same ij coordinate, present elsewhere along the z axis. Closed loops, including their feeding wires running in parallel to B 0 needed to fulfill Eq. 8 are sketched in Figure 1B. Such interconnections along z might even be realized without altering the resulting target field B z . On continuous surfaces of simple topologies this ensures a design which is realizable with closed current loops on this surface. In the practical implementation Eq. 8 can be converted to a matrix form and appended to Eq. 5 before solving the regularization problem.
Closed loop iso-contours derived from the stream function requires some manual interaction to get realizable gradient coil design, which can be manufactured by a single wire. Manual addition of connections between neighboring iso-contour lines transforms each quadrant into a spiral. Connections are usually chosen along the z-direction while return paths are positioned on top to mitigate effects from the introduced changes [9]. An automated process to interconnect multiple groups of streamlines has been proposed recently [10]. However, such an automatic approach is beyond the scope of this work.
Additional constraints may be added to enforce further design requirements. Exemplary mentioned here is the balance of force or torque which may be expressed by simple matrix operations, as well. Assuming thin wires which host currents I n orthogonal to a strong and unidirectional main magnetic field B z the force on each discrete element, n, is given by To account for field inhomogeneities, the main magnetic field B needs to be defined at each wire location. Accordingly, the excess torque may now be calculated by summation over all moments from force F at position r. The resulting term may be added and minimized during the optimization. A regularization parameter has to be added for an adequate weighting of the constraints. Balancing multiple regularization parameters is the art of coil design and usually requires a comparison of multiple designs within a reasonable parameter space [11].

From Current Distributions to Stream Functions
The thin-wire approximation provides an already powerful approach to find realizable gradient coil designs. However, the coil layout is not directly derived, but requires and intermediate step in the current domain. For calculating stream functions, the boundary element method is usually deployed. Most methods are based on triangular meshes, each triangle describing one boundary element. However, finite element methods (FEMs) in general are not limited to triangular meshes. Boundary elements with a rectangular mesh may be used, as well (e.g. [12,13]). Considering only the z-component of the magnetic field, a simplified rectangular element may be derived from two neighboring thin wires, as used in the thin-wire approach. Since connecting paths along the z-axis do not contribute to B z , they may be neglected. A basic cell from wire elements is sketched in Figure 1C.
Corresponding to the sensitivity matrix S (Eq. 4) a similar sensitivity matrix in the stream function domain S stream may be derived. Size of the original sensitivity matrix S is given by the number of target points, n, the number of circumferential segments, m c , and number of segments along z, m z , S n,m c ,m z . The stream function sensitivity S stream is given by spatial derivative along z of the sensitivity S. In a matrix-like notation the stream function sensitivity is defined as S stream S n,m ( : , 1: end − 1) − S n,m ( : , 2: end).
In analogy to the inverse problem stated in Eq. 5 a corresponding current value may be calculated. For regularizing this equation an effective current has to be used which consists of the combination from currents of two neighboring cells at the same location which can be achieved by an appropriate modification of the matrix Γ.
As previously shown, closed loops were enforced in the current domain by a constraint which was defined in a global manner (see Eq. 8). Due to the definition of the stream-function as spatial difference between two neighboring current carrying Frontiers in Physics | www.frontiersin.org September 2021 | Volume 9 | Article 699468 elements the same constraint does not lead to a meaningful result. Therefore closed loops have to be enforced using a different strategy. The approach chosen here was to demand the peripheral elements to carry zero current. This requires the addition of further terms to the inverse problem. Same currents of all outer elements, in this case zero, ensures closed stream lines. In analogy to the closed loops, further constraints may be added similarly. As stated above, this may include balance of force and torque, maximum current density or others.  An additional level of complexity is added for split gradient coils e.g. for multi-modality or combined MR-Linac systems [14,15]. Such coils in a shielded arrangement consist of four independent current carrying surfaces. Dimensions of the cylindrical shielded coil were used as a basis to simulate such a setup with radii of 0. The strength for each target field was set to 1 [mT/m] and the number of iso-contours was chosen to allow for a good visualization for plotting. It should be noted that the actual number of iso-contours is usually chosen based on constraints given by the overall system design. For larger coils the choice is a trade-off between efficiency and possible switching speed which is limited by the available maximum current and voltage. For smaller coils the number of iso-contours is mostly limited by the available space which is needed for a finite wire thickness.

RESULTS
To demonstrate the versatility of the proposed method, stream functions based on discrete current distributions have been calculated for multiple different numbers of independent current carrying surfaces and different coil geometries. A gradient coil design on a single surface is given by the cylindrical design shown in Figure 2. The original target field along with the deviation of the resulting field is depicted in Figure 3. The regularization chosen here results in a deviation of about ±6%. This deviation may be controlled by the relation between the regularization and the chosen dimensions.
Cylindrical shielded gradient coils are usually used in high-field MRI scanners to suppress eddy currents in the cryostat [16][17][18]. Therefore, a second target region is defined on a cylindrical surface at the position of the cold radiation shield. A wire mesh along with the main and shield target is shown in Figure 4. Compared to the main layer, the shield exhibits opposite polarity which corresponds to opposite direction of currents. A biplanar set of coils deploys two independent current carrying surfaces and is depicted in Figure 5. One should keep in mind that most biplanar setups in the literature are designed for a main magnetic field which is oriented along the vertical axis. However, this is not the case in the example shown here.
An additional level of complexity is added in the design of split coils e.g. for multi-modality or combined MR-Linac systems [14,15]. Such coils in a shielded arrangement consist of four independent current carrying surfaces. Resulting wire layouts are shown in Figure 6. Same dimensions as for the previous shielded design and a gap of 0.2 m was chosen.
One of the main results is the code which was used to generate the rectangular thin-wire mesh and all designs within this manuscript. The full code is available as supporting material and under 1 .

DISCUSSION
One of the main motivations for this work is the educational aspect, introducing gradient coil design in a straight-forward, intuitive and FIGURE 5 | Bi-planar design to demonstrate the feasibility of using flat surface geometries. Direction of the main magnetic field is along the z-axis. Both surfaces (first and second row) along with a 3D visualization are shown for three different gradient fields. Currents on all edges of both surfaces had to be forced to zero for achieving closed loops. easy to understand way. To do so, we approach the gradient coil design from a current-density perspective. If the magnitude of the encoding field is much smaller than the main field strength a simplified representation becomes possible by considering the z-component of the magnetic field, only. Hence one field component instead of all three has to be calculated while significantly reducing the computational effort. Versatility of the proposed method is demonstrated by coil designs on flat and cylindrical surfaces. Coil designs on multiple independent current carrying surfaces are demonstrated by a shielded and a split shielded design.
The coil design methods shown are closely related to previously shown boundary element approaches, e.g. [3,4] which have been maximally reduced. Due to the reduction to B z only, the optimization may not be usable directly for all magnet types. For classical and superconducting electromagnets the current carrying surface(s) are usually orthogonal to the main magnetic field, regardless of field strength. MRI systems based on rare Earth magnets are usually orientated such, that the current carrying surfaces of the gradient coils are not perpendicular to the direction of the main magnetic field. This includes u-shaped and Hallbach array type magnets [19][20][21][22] and others. However, the reduced form presented here may be easily adapted by adding longitudinal thin-wires to consider all field components, including B xy , and closing the boundary element loops. Such a possible addition comes with an increased computational effort.
Due to the simplification using only wires orthogonal to the zaxis, the current method is limited to regular surface geometries. Methods for arbitrary surfaces have been shown e.g. in [5,23]. It should be noted that for manufacturing reasons most gradient coils are realized in fact on simple surfaces. However, with new manufacturing techniques based on 3D printed materials, such as shown in [24] other shapes might be used more often in the future. For such cases a full description of all components of the current density is also possible. For instance in the literature on the FEM methodology, approaches to shear and deform rectangular elements are usually covered, e.g. [12,13]. In the context of gradient coil design non-triangular mesh elements have been used, as well [25,26].
Minimum wire distance is one of the engineering constraints which has to be considered during the design. This has not been considered in the designs presented here and the associated problems are most apparent in the bi-planar designs in Figure 5. One approach is to choose the number of coil turns according to spatial constraints. More sophisticated approaches constrain the gradient of the stream function, the current density, using a minimax design or a p-norm [6,27]. Alternatively, an explicit constraint could be relatively trivially added directly to the thin-wire solver.
As a further step the iso-contour lines of the stream function have to be interconnected. A automatized interconnection algorithm has been demonstrated recently [10]. We hope that a combination of the intuitive approach to coil design shown here with an easy to use interconnection algorithm lowers the boundary for making experimental gradient coils.
A further goal of this publication was to make gradient coil design code which covers main aspects of the coil design problem, including surface mesh definition, available to the public. In addition to the supplemented materials, the code is available on GitHub under 1 . We hope that additional adaptations of the code are available in the future, including extensions to arbitrary orientations of the main magnetic field, e.g. for systems based on Halbach type magnets.