Magnetic-field modeling with surface currents: Physical and computational principles of bfieldtools

Surface currents provide a general way to model static magnetic fields in source-free volumes. To facilitate the use of surface currents in magneto-quasistatic problems, we have implemented a set of computational tools in a Python package named bfieldtools. In this work, we describe the physical and computational principles of this toolset. To be able to work with surface currents of arbitrary shape, we discretize the currents on triangle meshes using piecewise-linear stream functions. We apply analytical discretizations of integral equations to obtain the magnetic field and potentials associated with the discrete stream function. In addition, we describe the computation of the spherical multipole expansion and a novel surface-harmonic expansion for surface currents, both of which are useful for representing the magnetic field in source-free volumes with a small number of parameters. Last, we share examples related to magnetic shielding and surface-coil design using the presented tools.


Introduction
Modeling magnetic phenomena with surface currents has various applications in physics and engineering. One large field of applications is surface-coil design, where continuous surface currents are used to design coil winding patterns. This design method has most prominently been utilized in plasma physics [1,2,3] and in magnetic resonance imaging (MRI) [4,5,6,7], but it has also been used in transcranial magnetic stimulation (TMS) [8,9], magnetic particle imaging (MPI) [10], and in zero-field magnetometry [11]. Further applications of surface-coil design include, .e.g., field control in physics experiments [12,13] and pickup coils of magnetic sensors [14,15].
The methods used in coil design are also involved in modeling eddy current patterns induced in thin conductive sheets [16,17,18] and field fluctuations due to thermal noise currents [19,20,21]. In addition, surface currents could be applied as equivalent models in magnetic shielding with highpermeability materials [22,23]. Modeling the magnetic field in free space using equivalent current densities on the volume boundary could also have various other applications. This method can be used directly for modeling the field of uniformly magnetized bodies [24], but could also be used as an equivalent model when interpolating magnetic-field data in e.g., in geomagnetism [25] and biomagnetism [26,27]. Additionally, such a field model could be applied, for example, when modeling magnetic fields for interference rejection [28].
Although, surface currents are used in the theory of static magnetic fields, their application in field modeling has been Email address: antti.makinen@aalto.fi (Antti J Mäkinen) limited because of a lack of general computational tools. To facilitate the application of surface-current densities in magneto-static problems, we introduce a novel Python software package bfieldtools (available at https://bfieldtools.github.io).
This package provides tools for representing currents on arbitrarily-shaped surfaces, and calculating the associated magnetic field and potentials. Further, tools for designing current patterns that generate desired magnetic fields are included. The whole software package is described in two papers. In this paper, we present physical and computational principles of the software and applications that showcase the capability of the presented tools. The accompanying paper [29] describes the Python-based implementation in detail and gives examples of the software use in different applications from the user perspective.
In bfieldtools, surface currents are represented by scalar stream functions [4,16,30,17], which effectively constitute the degrees of freedom required for modeling divergencefree currents on surfaces. We discretize the stream functions on a triangle mesh using piecewise-linear interpolation functions equivalent to piecewise-constant surface current as, e.g., in Refs. [4,30,5,31]. Compared to analytical methods (e.g., [32,33,34]) that require certain symmetries for the current distributions, discretizing the stream function on a freely shaped triangle mesh allows studying currents and magnetic fields in a wide range of geometries.
In this work, we first review the physics of the stream function. As an additional feature to previous works, we relate the stream function to harmonic potential theory used, e.g., in bioelectric volume-conduction problems [35]. By introducing the magnetic scalar potential to the computational framework, analogies to other fields utilizing potential theory can be exploited, facilitating the formulation and solution of magneto-quasistatic problems.
The main objective of this work is to describe the field calculations and their discretization as they are implemented in bfieldtools. Based on previous studies that utilize the linear discretization of the field source [36,37,38,4], we obtain a consistent analytical discretization of the integral equations involved in field calculations. The same principles can also be used to obtain discrete differential operators on a surface [39,40], which we utilize in the integral formulas.
We have also implemented computations for series representations of the magnetic field in a source-free volume. First, we review the multipole expansion in terms of spherical harmonics, which is the conventional way of describing such a field. We adapt the multipole expansion of 3D current densities [41,28,42] to obtain the expansion for the field from a surface current on a mesh. In addition, we introduce a novel field representation based on expanding the stream function with the eigenfunctions of the surface-Laplacian [43,44], which can be seen as a generalization of the multipole expansion.
Finally, we share a few examples demonstrating the capability of these tools in coil design and magnetic shielding. More applications are described in the accompanying paper [29], including references to the software implementation.

Divergence-free surface currents
A divergence-free current density j ( r ) on a surface can be expressed with a scalar stream function ψ on the surface [16,30,17] j where r is the position on the surface,n is the unit surface normal, ∇ the tangential gradient operator, i.e., the 3D gradient projected to a tangent plane on the surface: ∇ = ∇ −n(n · ∇) [45]. As ∇ ψ( r ) ×n( r ) is perpendicular to ∇ ψ( r ), the streamlines of the current correspond exactly to the isocontours of ψ( r ). For convenience, we define the operator ∇ (·) ×n as the rotated gradient. By taking a path integral of ψ from r 0 to r on the surface, we find that the difference in the stream function between the two ends of the path equals the flux of surface current j ( r ) passing the curve [16] where d l ×n is a path differential perpendicular to the direction of the path. In consequence, a path integral from a reference point r 0 determines the stream function uniquely on the surface. From another point of view, the stream function can be interpreted as a surface density of magnetic dipoles normal to the surface (see, e.g., [46,47] and Appendix A) This formula corresponds to the continuous limit of dividing the surface currents into smaller and smaller loops that eventually correspond to infinitesimal magnetic dipoles [48,49]. This interpretation of the stream function enables analogies to dipole layers involved in, e.g., volume conductor problems and the calculation of magnetic scalar potentials for divergence-free surface currents. As the surface current density is assumed divergence-free everywhere, the flux of current through any boundary on the surface must be zero. Applying Eq. (2) on a boundary, we can deduce that this condition is equivalent to ψ( r ) being constant on the boundary. With only one boundary, the constant can be set to zero since an additional constant in ψ( r ) does not affect j ( r ). When the surface contains holes, the hole boundaries can have their own constants. These matters are further discussed in Sec. 3 when discretizing the stream function.

Stream functions and the magnetic scalar potential
The magnetic field B originating from sources outside the volume of interest can be expressed as the gradient of a scalar potential U : B = −µ 0 ∇U . As the magnetic field is divergencefree, the magnetic scalar potential U is harmonic, i.e., it satisfies Laplace's equation ∇ 2 U = 0. From the theory of harmonic potentials [48], we know that U can be determined uniquely in the volume (up to a constant) when either the potential or the normal derivative of the potential is specified on the boundary enclosing the volume. Thus, any external source distribution whose potential reproduces the boundary conditions of a given U , can be used to generate U in the volume.
In particular, the boundary condition can be satisfied by the potential of a dipole density ψ( r )n on the same surface [50,51]. In potential theory, this source distribution is known as a double layer, equivalent to two parallel layers of opposite charge. In magnetostatic calculations, as discussed above, such a layer of magnetic dipoles corresponds to a surface-current density ∇ ψ( r )×n. Any magnetic field within a source-free volume can thus be expressed with a stream function on the boundary of the volume.
As the discussion above applies only to a closed surface, a stream function on a surface with openings cannot generally represent all possible field patterns in the volume. This must be taken into account in coil designs where the current may only be placed in restricted regions as well as in fieldinterpolation tasks with equivalent surface currents. However, the dipole-layer analogy still applies to a stream function on an open surface: the stream function always corresponds to a discontinuity in the scalar potential [17] similar to a dipole layer [50].

Integral equations
In the following, we layout the integral equations for calculating the quasistatic magnetic field and magnetic potentials from a stream function. The integrations are discretized in Sec. 3.
In current-free volumes, the magnetic field has both vector and scalar potential representations [48] The vector potential of a surface current density can be written as an integral over the surface S, where the stream function is defined: The vector potential can be equivalently written in terms of a magnetic dipole layer m = ψn (see Appendix A), which is the more convenient form to express the the magnetic scalar potential Similar dipole-layer potentials are used for the electric field in volume conductor problems [35]. Finally, the Biot-Savart formula for the magnetic field is obtained as the curl of the vector potential In computations, it is useful to expand the stream function with a set of basis functions ψ k ( r ) as The coefficients s k parametrize the stream function, enabling linear-algebraic techniques for processing it. Furthermore, the basis functions ψ k ( r ) can be made to satisfy possible boundary conditions so that any combination of them satisfies the same conditions. In some geometries, ψ k ( r ) can be chosen as, e.g., sinusoids or spherical harmonics [17,46,13]. The rotated gradients of the basis functions provides a basis set of vector functions j k ( r ) = ∇ ψ k ( r ) ×n( r ) that expand the current density. The basis function coefficients s k , forming a column vector s, can be used to write the inductive energy and resistive dissipation power of a surface current as quadratic forms of s [46,17,10]. The inductive energy, i.e., the energy stored in the magnetic field, can be written as s M s/2, where the matrix M consists of the mutual inductances of the current patterns, which can be calculated as [48] where A l and B l are the magnetic vector potential and the magnetic field generated by the current pattern j l ( r ), respectively. The power dissipation due to resistive heating can be written as s R s, where is the mutual resistance associated with the two current patterns. Here, e k is the electric field associated with j k and σ s ( r ) = σ( r )d ( r ) is the surface conductivity defined by material conductivity σ and surface thickness d . Further, assuming constant surface conductivity and using stream functions to describe j k and j l , we get Partial integration (Gauss theorem) was used to get the last identity, where ∇ 2 = ∇ · ∇ is the surface Laplacian or the Laplace-Beltrami operator [44,43]. The possible boundary terms vanish similar to derivation in Appendix A. The relationship between the mutual resistance and the Laplacian is utilized further in the next section.

Piecewise-linear stream function
In bfieldtools, surface-current densities are represented by stream functions on triangle meshes. A triangle mesh consists of an ordered collection of vertices r 1 , ..., r V , forming a point cloud in a 3D space, and of a set of triangular faces ∆ f , each defined by a triplet (i , j , k) of vertex indices. We discretize the integral and differential equations described in the previous section by approximating the stream function as linear on each face of the triangle mesh similar to, e.g., Refs. [5,30]. Such piecewise-linear functions can be conveniently expressed as in Eq. (8) by choosing the basis functions ψ k ( r ) to be so-called hat functions h i ( r ), where the index i corresponds to the i th vertex of the mesh. The hat function h i ( r ) has the value one at vertex i and zero at all other vertices. Within triangles, the value is interpolated linearly (see Fig. 1A). As in Eq. (8), the stream function can be written as a sum of the basis functions where s i , the weight for the vertex i , is equal to the current circulating around the vertex on the neighbouring triangles. We obtain current-density basis functions by taking the rotated gradient of the hat function j i ( r ) = ∇ h i ( r ) ×n , which corresponds to an eddy current circulating around vertex i as illustrated using black arrows in Fig. 1A.
In each neighbouring triangle, the gradient and rotated gradient of a hat function are constant and can be expressed using the local geometry [39,40] as where A f is the area of the neighbouring triangle f ,n f the triangle normal and e i the edge opposing the vertex i in the triangle. The constant condition on the outer mesh boundary can be implemented by setting the boundary-vertex values to zero. However, each hole boundary can float at an arbitrary value. To satisfy the constant boundary condition on the hole boundaries, we construct a combined basis function for each hole boundary C k as These functions are constant along the hole boundaries and can be conceptualized as a current flowing around the hole within the triangles neighbouring the hole vertices (Fig. 1A).
The stream function can now be expressed as where the first part sums over the inner vertices of the mesh and the latter sums over the holes of the mesh. With this vertex-wise discretization of the stream function, we can represent physical quantities using operators acting on the vertex values s i . Stacking the weights s i into a column vector s, linear operators (e.g. the surface-Laplacian) acting on ψ discretize to matrices that can be used in linear mappings s → As or in quadratic forms s → s As. Additionally, fields originating from the discretized current can be expressed as B ( r ) = n B n ( r ) s n = B ( r ) s, where B ( r ) is a column vector of the magnetic field contributions at r from each vertex in the mesh.

Differential operators
With the hat-function discretization, the gradient of any scalar function can be calculated on the faces of the mesh from the neighbouring vertex values [40,39]. As this calculation is linear with respect to the vertex values, we define a discrete gradient operator G as a map from scalar values at the vertices (i ) to Euclidean vectors at the faces ( f ). Since there are only three non-zero hat functions on each triangle, the result of this operation can expressed as where i , j , and k are the vertices of ∆ f . One element of the operator is obtained directly from the gradient of the basis function Eq. (13) as The elements of the rotated-gradient operator are defined as Using the hat functions to discretize the surface-Laplacian operator leads to the so-called cotan formula derived and applied in the context of partial differential equations as well as in geometry and graphics processing [52,53,54,55]. As second derivatives are ill-defined for hat functions (zero on the faces, infinite on vertices and edges), the discrete Laplacian operator L is understood as the weak (integrated) form of the surface Laplacian: Using Eq. (18), the non-zero off-diagonal elements of L can be expressed as where i and j correspond to two neighbouring vertices, angles α i j and β i j are the angles opposing the edge connecting Figure 2: Geometric quantities in the discretization of inductance, resistance and Laplacian operators. Laplacian and resistance matrices are sparse: only the vertices whose neighbouring triangles overlap contribute to to the matrix elements. The inductance matrix is dense with elements describing the coupling between the currents circulating the vertices.
the vertices, and indices 1 and 2 correspond to the two triangles that share the edge, as illustrated in Fig. 2. The diagonal elements are then given by so that each row (or column) of the matrix sums to zero. When the surface has boundaries (outer edges or holes), the Laplacian has to be modified. Elements that correspond to the zero-valued boundary can be left out of the matrix. Using Eqs. (15) and (16), it can be further shown that the elements that correspond to the basis function of a hole boundary can be obtained by summing the rows and columns associated with the vertices on the boundary.

Analytical integrals
The analytical integral formulas introduced in this section are the basic building blocks of bfieldtools mesh operators for the magnetic field and magnetic potentials (Sec. 2.3). As the integrals needed to compute the mesh operators involve singular quantities, the analytical formulas behave more smoothly in the proximity of the source mesh compared to numerical quadratures. These formulas have been derived in the literature related to boundary-element methods in bioelectromagnetism [36,37,38] and in antenna modeling [56,57,58]. To introduce concepts and to unify notation, we review the analytical formulas, which can also be seen as potentials of simple charge or dipole configurations visualized in Fig. 3A.
The first building block, used in all the field calculations, is the solid angle subtended at r by triangle ∆ f consisting of vertices i , j , and k. The integral corresponds to the scalar potential of unit (magnetic) dipole density on a triangle and can be computed as [36] where the denominator for the two-argument inverse tangent function atan2 defined as in most standard programming languages. Here, d i = r − r i is a vector pointing from the vertex i to the evaluation point r (see Fig. 1B). The second building block is the potential of a line charge. Following Refs. [37,38], we define the potential of a unit line charge on edge e i as where r j and r k are the two ends of edge e i , which is opposite to vertex i . In Refs. [37,38], this potential is known as γ 0 i to distinguish it from the potential of linearly varying line charge density (γ 1 i ). With the two integrals above, we can express the potential of a unit charge density on a triangle ∆ f as [38] Here, d f =n f · d l +1 is the signed distance from the triangle plane along the plane normal and x l = G f ,l · d l +1 is the normalized signed distance from the line extended from the edge e l along G f ,l such that x l ( r l ) = 1. The geometry related to these distances can be found in Fig. 1C. In addition to the derivation of Eq. (25) in Ref. [38], we provide an alternative, more compact derivation of the formula in Appendix B using the notation of this work. As the last building block, we present the potential of a lin-early varying dipolar density h i ( r ) on a triangle ∆ f [37]: where c i ,l = e i · e l /(2A f ). A concise derivation of Eq. (26) is provided in Appendix B.

Magnetic field and magnetic potentials
We express the magnetic field and potentials using mesh operators such that, e.g., the magnetic field at r is where the sum is taken over the vertices of the mesh and B i ( r ) is the magnetic field corresponfing to hat function h i ( r ). When a set of field evaluation points { r j } is given, the operators can be expressed as coupling matrices, whose elements equal the coupling between the hat-function currents and the field components at the evaluation points.
The magnetic field B i ( r ) of a constant current density in a triangle is derived in Appendix B and, using that result, the magnetic field due to a single hat-function current becomes where N i denotes the set of triangles neighbouring vertex i as shown in Fig. 2. A corresponding formula expressed in local coordinates of a triangle can be found in Ref. [4]. The vector and scalar potentials for the current of a hat stream function h i can be obtained in a straightforward manner using the integrals in Sec. 3.3. The vector potential [Eq. (5)] can be expressed using the discrete rotated-gradient [Eq. (19)] and the integral φ f ( r ) [Eq. (25)] [4,9]: The scalar potential [Eq. (6)] of h i involves only the potentials of linearly varying dipole densitiesΩ f ,i ( r ) (26) The magnetic field and potentials due to a single hat-function current are illustrated in Fig. 3B.

Mutual inductance and resistance
The mutual inductance and mutual resistance for hat functions can also be expressed in terms of the integral and differential operators described above. The mutual inductance between two hat-function currents (Fig. 2) can be calculated using Eq. (9) as where the second integral is calculated using quadrature points r q with weights w q , as in Ref. [9]. We have implemented this approach in bfieldtools, as it naturally handles the singularity in the double integral when ∆ f = ∆ f . Alternatively, the singularity can be handled with an analytical formula for the self element [59] as in Refs. [30,5,31]. In this basis, the mutual inductance operator M can also be interpreted as a mapping from the discretized stream function s to the magnetic flux (integrated normal component) at the mesh vertices. This can be seen from the second identity in Eq. (9): by replacing m k with the dipole densitynh k , the matrix element M k,l corresponds to the normal magnetic field of current l integrated over the hat function of vertex k.
For mutual resistance we also have to model the surface conductivity σ s . Assuming piecewise-constant surface conductivity on the triangles, and using Eq. (10), we obtain the mutual resistance operator as where σ 1 and σ 2 are the conductivities in triangles 1 and 2 neighbouring the edge from vertex i to vertex j (see Fig. 2). When σ is constant over the surface, the resistance operator is proportional to the discrete surface Laplacian Several studies, e.g., [30,5,47] have utilized formulas equivalent to Eq. (32), although the relation to the Laplacian has not been noted.

Magnetic field representations with source expansions
The magnetic field in free space can be expanded as a series of components, each of which can be interpreted to correspond to a certain type of a source-current pattern. A common series used for the static magnetic field is the spherical multipole expansion. This expansion can represent a spatially smoothly varying magnetic field with a few parameters, which can be helpful, e.g., when designing coils that generate these types of fields [42,60,13]. A disadvantage of the multipole expansion is, however, that the series can diverge in regions where the actual field is wellbehaving. For more general purposes, we introduce a representation of magnetic fields based on a stream function, The magnetic scalar potential corresponding to different components of the spherical harmonic (SPH) and surface harmonic (SUH) inner-source expansions. The scalar potential is depicted by the red-blue colors on the vertical plane and the field source is illustrated by the green-brown color. In both expansion, the first, fifth, 11 th , and 16 th component of the series are shown. Compared to the multipole fields, the fields of the SUH expansion are distributed more uniformly around the surface. The bunny surface mesh was decimated from the original (http://graphics.stanford.edu/data/3Dscanrep/#bunny) and repaired for holes. which can be viewed as an equivalent source of the field. We expand the stream function on a surface as a series of functions that we call surface harmonics. The magnetic field patterns of the surface harmonics then yield a representation of the field similar to the multipole expansion.

Multipole expansion with spherical harmonics
The general solution of Laplace's equation in spherical coordinates ( r , θ, ϕ) is [48] where α l m and β l m are the multipole coefficients, Y l m (θ, ϕ) are spherical harmonic functions with degree l and order m (|m| ≤ l ). The α l m terms involve powers of the inverse distance, representing sources close to origin, whereas the β l m terms represent far-away sources. In bfieldtools, we use the real spherical harmonics, which are orthonormal with respect to integration over the full solid angle, i.e, Ω Y l m Y l m d Ω = δ l l δ mm . The exact definition of these functions can be found in Ref. [61].
We obtain the expansion for the magnetic field by taking the gradient of the scalar potential: where V l m = (−(l + 1)Y l mr + ∇ 1 Y l m )/ (l + 1)(2l + 1) and W l m = (l Y l mr + ∇ 1 Y l m )/ l (2l + 1) are vector spherical harmonics [62,28]. Here, ∇ 1 is the angular part of the gradient on a unit sphere [61]. Both the set of vector spherical harmonics V l m , W l m , and the set of tangential vector spherical harmonics X l m = −r × ∇ 1 Y l m / l (l + 1) are also orthonormal with respect to an inner product Ω f l m · g l m d Ω .
Because the multipole expansion of B is linear with respect to the coefficients, we can express it using linear operators B α and B β as where the expansions coefficients, truncated at a certain degree l , are stacked in the column vectors α and β. When the magnetic field is evaluated at a specific set of evaluation points, the linear operators above can be expressed as coupling matrices that map the given multipole coefficients to field values at the evaluation points.
The coefficients α l m and β l m can be calculated directly from any surface-current distribution j ( r ) with help of the tangential vector spherical harmonics X l m as [41,42] α l m = l (l + 1) Using these equations together with the stream function in Eq. (1) and its discretization in Eq. (12), we define mesh operators (matrices) C α and C β that map the stream-function values to the spherical harmonic coefficients The convergence of the multipole series may be analyzed by inserting, for example, the inner multipole coefficients of Eq. (36) into Eq. (33) which yields terms involving factors (r /r ) (l −1) . If there are any sources with radius r greater than the radius of the field point r , the factors (r /r ) (l −1) approach infinity with growing l and the series fails to converge. A similar analysis can be made for the outer sources. These analyses result in convergence regions for the expansions shown in Fig. 4A. As the figures show, the choice of the origin is crucial for the convergence, but it cannot be chosen in such a way that the series would converge everywhere in the volume where the scalar potential is defined.

Surface-harmonics expansion
With the tools presented in this work, we can generate a field expansion that converges at all points where the magnetic scalar potential is defined [Fig. 4]. The expansion is based on the fact that any potential satisfying the Laplace equation can be written in terms of an equivalent stream function on a boundary of the domain as described in Sec. 2.2. Expanding the stream function as a linear combination of basis functions with increasing order of spatial detail yields a field representation similar to the spherical multipole expansion.
We base the expansion on the eigenfunctions of the (negative) surface Laplacian. These eigenfunctions generalize a sinusoidal function basis, such as the spherical harmonics basis, to an arbitrary surface (see Fig. 5). Generally, these functions are characterized by the eigenvalue equation [43,44] where the eigenvalue k 2 n corresponds to the squared spatial frequency of the n th eigenfunction v n . The higher the order n is, the higher the spatial frequency and the more zero crossings v n ( r ) has (Fig. 5). In relation to spherical harmonics, we call these functions surface harmonics (SUHs). In geometry processing, they are also known as manifold harmonics [63].
For practical computations, we discretize the functions as v n ( r ) = i V i ,n h i ( r ), which leads to a discrete (generalized) eigenvalue equation [44] − Lv n = k 2 n N v n , where L is the Laplace operator in Eq. (21) and N is a mass matrix taking into account the overlap of hat functions, and v n correspond to columns of the matrix V. As both L and N are sparse matrices, the vertex coefficient vectors v n can be solved efficiently with sparse solvers. The resulting eigenfunctions v n ( r ) are orthonormal with respect to integration over the surface, which can be expressed in the discrete form as v n N v m = δ n,m . Substituting the surface-harmonics representation of a stream function ψ( r ) = n a n v n ( r ) to Eq. (6), we can write the magnetic scalar potential as Similarly, the SUH coefficients a can be mapped to the magnetic field as Examples of the scalar potentials of the basis functions v n (SUH) are displayed in Fig. 4B with a comparison to the multipole expansion (SPH). The SUH expansion is not restricted to closed surfaces, but can be applied for stream functions on surfaces with boundaries and any number of holes. Such bases can be used for surface-coil design to decrease the degrees of freedom when optimizing surface currents, as demonstrated in the accompanying paper [29].
Instead of orthogonal stream functions, one may desire orthogonality in their magnetic fields. In that case, the basis functions can be solved from an discrete eigenvalue equation similar to Eq. (41) by replacing the mass matrix N with the inductance matrix M. To enable physical interpretations, L can be replaced by R = L/σ s to get the following eigenvalue equation This equation is related to the independent modes of eddy currents on the conducting surface; τ n is the time constant of the nth mode. These modes can be used, for example, to calculate the time dynamics of eddy-current induced fields [17] or uncoupled current patterns for thermal noise calculations [19,64]. It should be noted, however, that M is now a dense matrix, whereas N was very sparse, disabling the use of sparse eigensolvers and increasing computation time when building the matrix.

Coil design and shielding
In this section, we give examples that utilize the developed tools. As the design of surface coils using distributed currents Figure 6: Two examples of surface currents on two surfaces, generating a desired magnetic field confined in a closed volume. Top: Stream functions of two surface-current configurations with a primary current on the sphere and a shielding current on the rounded cube. The currents are designed so that together they create a homogeneous field (left) and a first-order gradient field (right). Bottom: the magnetic field lines (black) and the corresponding magnetic scalar potential are plotted on the horizontal plane shown in the 3D plots on top.
is probably the most prominent application of these tools, we start by giving a brief overview of the coil-design method. In surface-coil design using triangle meshes, the coil current is expressed with a discretized stream function s on the mesh similarly as in our tools (e.g., [30,5,31]). The stream-function s is optimized by minimizing a cost function while taking into account given constraints for, e.g., the field shape. Finally, the coil wires are placed on the isocontours of the stream function to approximate the continuous current density. Typically, a quadratic form of s such as the magnetic energy s M s or the dissipated power s R s is used as a cost function [30,5,8]. Constraints for the field pattern can be formulated using the mesh operator B ( r ) and they can be incorporated in quadratic programming as demonstrated in the accompanying paper [29].
Here, in the next examples, we take a more theoretical approach to surface-coil design and, in particular, to the design of self-shielded coils. We also share an example of a calculation related to magnetic shielding using the tools described in this work.

Perfect shielding by surface currents on a closed surface
When designing coils for target magnetic fields, it is often also desired to control the field outside the volume of interest, e.g., to shield the external environment from the field of the primary current. For such a situation, a shielding current outside the primary surface can be designed.
Let us consider a closed surface, inside which a desired field pattern is to be designed and a second (outer) surface, the exterior of which is to be shielded from the field. To derive a set of equations with a unique solution, we again discretize the surfaces using triangle meshes. The shielded field pattern can be obtained by designing suitable stream functions s 1 and s 2 on the two surfaces so that their combined field satisfies desired boundary conditions. Based on the discussion in Sec. 3.5, we can write the boundary conditions for the normal component of the field at the surfaces using stacked mutual-inductance matrices M i j as where the first row corresponds to the desired magnetic field b n at the inner surface and the second row corresponds to the zero condition for the outer surface. Because s M s is the total magnetic-field energy of the system, M is a positive semidefinite matrix. As the only zero eigenvalue is the one corresponding to a constant stream function (zero current), the system can be solved by inverting M deflated for the constant vector. Two examples demonstrating the perfect shielding obtained by solving Eq. (45) are shown Fig. 6. The shielding by the outer surface corresponds to the situation where the exterior volume would be a superconductor that expels all fields so that no magnetic field crosses the outer surface. This is also analogous to an electrical volume-conductor problem where the exterior volume is insulating, confining the 3D current.
When the outer current surface contains current-free regions in it, perfect cancellation of the primary field is generally not possible due to a lack of degrees of freedom in the current-pattern design. Next, we demonstrate the computations of the primary and shielding currents applicable also to an open surface geometry.

Self-shielded currents with an open geometry
We now apply the tools presented in this work for the design of self-shielded coils in a more realistic bi-planar geometry. Consider a primary coil with stream function s 1 and a shielding coil with stream function s 2 . In Ref. [6] it is demonstrated that a well-performing shielding coil can be designed by minimizing the magnetic-field energy with respect to s 2 . Using s 1 and s 2 , the field energy can then be expressed as The minimum (for given s 1 ) can be found by equating the gradient of E M with respect to s 2 to zero, which yields a set of linear equations: Based on Eq. (45), we can now explain why this method works: the equation can be interpreted as a condition that the We will now optimize the primary current s 1 for minimal energy with an additional constraint. Instead of a hard equality constraint for the desired field, we modify the cost function with a term that penalizes for the residual in desired multipole moments β: where λ is a trade-off parameter between a perfect multipole fit and a minimal field energy. The coupling matrix C β = C β,1 + C β,2 M −1 22 M 21 is obtained using the constraint in Eq. (47); C β,i are matrices that map s i to the multipole moments. The solution for s 1 that minimizes the cost function can again be found by equating the gradient ∇ s 1 E (s 1 ) to zero: Examples of shielded configurations generated using the method above are shown in Fig. 7. Compared to the fields in Fig. 6, which are solved for a closed geometry, these fields leak in the directions where the shielding surface is missing. Thus, the placement of the shielding surfaces is crucial for the selfshielding performance.

Modeling a high-permeability magnetic shield
Magnetic measurements are usually shielded from the lowfrequency fluctuations of the outside magnetic environment with soft ferromagnetic materials such as µ-metal. The purpose of these materials is to guide the external magnetic field to create a magnetic void inside the shield. A downside is that the shield also distorts the fields generated inside the shield. When the relative permeability of the shield is very large, the effect of the shield can be approximated by a shield with infinite relative permeability. This leads to a boundary condition stating that the magnetic field must be normal to the inner surface of the shield or, equivalently, the inner surface has to be at equipotential in terms of the magnetic scalar potential [48].
Let us consider a primary potential U p ( r ) s p generated by a surface current with (discretized) stream function s p inside the magnetic shield. We can satisfy the equipotential condition on the shield by placing a suitable equivalent surface current s eq on the shield surface. In other words, we require that U p ( r )s p + U eq ( r ) s eq = 0 holds on the shield. To solve for s eq , we apply the condition at collocation points on the shield mesh. Because U eq ( r )s eq is discontinuous across the shield, we apply the condition at collocation points slightly inwards from the shield surface as r j − n j , where r j are the vertex positions of the shield mesh, is a small number compared to the mesh resolution and n j are the vertex normal vectors. After solving the resulting set of linear equations for s eq , we can estimate the effect of the shield by U eq ( r ) s eq at all points inside the shield. Figure 8 shows an example of a magnetically shielded configuration with a surface-current pattern on a bi-planar surface inside a perfect cylindrical magnetic shield. The current patterns on the planes are designed such that the field in the The equivalent surface-current distribution (stream function and surfacecurrent density) representing the induced field source on the µ-metal shield. C: The primary magnetic scalar potential generated by the primary source in A. D: The magnetic scalar potential generated by the equivalent current in B inside the shield surface. E: The combined potential that satisfies the constant-potential boundary condition. The dashed circle represents the volume of interest. target volume, indicated by the dashed circle, is as homogeneous as possible. The secondary field due to the shield amplifies the magnetic field (the gradient of the potential) in the vicinity of the shield outside the target volume. The contribution from the cylinder cap also produces minor inhomogeneity in the magnetic field inside the target volume.

Discussion and outlook
We have introduced a set of tools for static and quasistatic modeling of divergence-free surface currents and their fields. The tools can be used for a variety of tasks from surface-coil design and equivalent-source modeling to eddy-current and thermal-noise calculations. This work has covered the central computations implemented in the software, the structure of which is described in the accompanying paper [29].
The computational and theoretical framework leverages the interpretation of stream functions as magnetic dipole densities normal to the surface. This analogy has been recognized previously [46,47], but has not been fully exploited. In this work, we exploit this interpretation further for the calculation of the magnetic scalar potential of a stream function, which we use, e.g., for visualizing the magnetic field. The scalar potential also enables the application of harmonic potential theory commonly applied in volume-conductor problems in the form of boundary-element methods (BEM) [65,66]. The discretizations implemented in bfieldtools are directly applicable for BEM computations. Namely, the potential of a linearly varying dipole density can be used to calculate the double-layer operator D for linear (hat) basis functions, and the potential of a constant charge density can be used for the single-layer operator S with constant basis functions. Additionally, the mutual inductance operator M is equivalent to a yet another operator called N , which maps a dipole density to the normal field component [50,65].
In the future, it can be fruitful to exploit the analogy between quasistatic magnetic problems and electric volume conductor problems even further. The source of the field in the former can be interpreted as magnetic dipoles, whereas in the latter, the sources are current dipoles. The magnetic scalar potential is analogous to the electric potential (both satisfy Laplace's equation), and if we interpret the permeability µ as the counterpart of electrical conductivity σ in a volume conductor, the magnetic field is perfectly analogous to the volume current density. The vector potential in the magnetic problem further corresponds to the magnetic field in a volume conductor problem. This means that the tools presented in this work could be applied to solve the electric potential in a volume conductor, e.g., for modeling transcranial magnetic stimulation [8,9] or solving the bioelectromagnetic forward problem [66].
The multipole and surface-harmonic expansions implemented in bfieldtools are also suitable for applications in a more general context, e.g., in biomagnetic experiments or geomagnetism. The multipole expansion has been applied to source modeling in bioelectromagnetism [35,67,68,69]. In magnetoencephalography, it has also been applied in signal space separation (SSS) [28], which can be used to design software spatial filters to reject external interference fields. In principle, the surface-harmonic expansion could be used for the same purpose with more general convergence properties.