Unique compact representation of magnetic fields using truncated solid harmonic expansions

Abstract

Precise knowledge of magnetic fields is crucial in many medical imaging applications such as magnetic resonance imaging (MRI) or magnetic particle imaging (MPI), as they form the foundation of these imaging systems. Mathematical methods are essential for efficiently analysing the magnetic fields in the entire field-of-view. In this work, we propose a compact and unique representation of the magnetic fields using real solid spherical harmonic expansions, which can be obtained by spherical t-designs. To ensure a unique representation, the expansion point is shifted at the level of the expansion coefficients. As an application scenario, these methods are used to acquire and analyse the magnetic fields of an MPI system. Here, the field-free-point of the spatial encoding field serves as the unique expansion point.

1. Introduction

Magnetic fields have been instrumental in the advancement of technology since the invention of the compass. They are fundamental for electric generators or motors, transformers and magnetic storage devices. In the field of medical applications, magnetic fields are the basis of various imaging systems like magnetic resonance imaging (MRI) or magnetic particle imaging (MPI). The precise generation of the magnetic fields has a significant impact on image quality, as even small deviations can lead to image artefacts and misdiagnoses. If the magnetic fields are known, the negative influence of these imperfections can be corrected in most of these applications, like field-related artefacts in MR images [25]. A standard method for magnetic field representation is a spherical harmonic expansion, which can be obtained via a calibration measurement of the magnetic field at several positions on a spherical surface [36]. This offers a robust and compact representation of said fields within a spherical region and allows analysis and solution of related problems [15].

In this paper, we will use MPI as an example to illustrate the strengths of spherical harmonic expansions for the analysis of magnetic fields. However, spherical harmonic expansions are not limited to this particular field. In fact, they have applications in a wide variety of other fields. As they provide a compact representation of magnetic fields, they are used in MRI to effectively design active or passive shimming [35], to determine the magnetic coupling between two electromagnetic sources [51] or to model the earth’s lithospheric magnetic field [33, 47]. Furthermore, spherical harmonics can be used for registration of objects by determination of the object’s orientation [12] or for simulation of high-resolution, full-sky maps of the cosmic microwave background anisotropies [34].

As with MRI, the fundamental building blocks of the recent imaging modality MPI are magnetic fields. It was already shown by Bringout et al. [10, 11] and Weber et al. [53, 52] that the coefficients of a spherical harmonic expansion are suitable for the representation of magnetic fields in MPI. Static magnetic fields spatially encode the MPI signal, while dynamic magnetic fields are used for signal generation [38]. MPI scanners are characterised by the topology of their static signal encoding field, which is either a field-free-point (FFP) or a field-free-line (FFL) [28]. Many reconstruction methods in MPI require some assumptions or knowledge about the magnetic fields. In x-space reconstruction, the position of the FFP is required to grid the measured signal to the FFP positions in the spatial domain and obtain the reconstructed image by a subsequent deconvolution [19]. The system-matrix reconstruction, where an inverse problem formulated with a dedicated system matrix is solved, needs the underlying magnetic fields in order to calculate model-based system matrices [1]. Furthermore, multiple measurement-based methods require the knowledge of the magnetic fields. Enlarging the field-of-view (FOV) is done by multiple shifted measurements, called patches, with dedicated system matrices. To cover the desired extended FOV, the magnetic-field-dependent shifts and expanses of the individual patches must be known. In order to apply a fast implicit reconstruction where a single system matrix is reused for all patches, the centre positions of the shifted patches must be known to avoid severe artefacts [46]. Since the shifted patches experience different magnetic fields, not only the centre position but also the expanse of each patch varies. This information should be incorporated into the reconstruction method to further reduce artefacts [9].

In a typical application, a magnetic field can be represented by three spherical harmonic expansions, one for each spatial component. These are calculated using a quadrature rule for a sphere and magnetic field measurements at the corresponding quadrature nodes. The coefficients of each expansion depend on the size and position of the sphere, and the resulting expansion can be used to predict the field values at certain positions, whereas the coefficients allow for an analysis of global field properties. In order to compare the properties of individual field components or of fields generated by different field generators, it is essential that the corresponding coefficients of the expansion have been recorded on the same reference spheres. Depending on the specific set-up, this is not always feasible.

In such application scenarios, the use of a spherical harmonic expansion that is independent of the geometry of the measurement set-up would be advantageous. The present paper proposes such a unique representation. Specifically, it suggests shifting expansions obtained on different but overlapping reference spheres into a common point by mapping the coefficients corresponding to the original expansion point to the coefficients corresponding to the common expansion point. This mapping is based on the classical shift theorem of spherical harmonics [45], which states how the basis functions of the expansion map under translation. The primary mathematical innovation of this work is the derivation of this mapping. Additionally, we provide an example of how the coefficients of the unique representation can be utilised to characterise imperfections in the magnetic fields of an MPI system.

We have structured the manuscript as follows. To conclude the introduction, we introduce MPI in more detail as this serves as motivation throughout the paper. Afterwards, the theoretical part will start with a review of real solid spherical harmonic expansions as a general solution of Laplace’s equation. This will provide the mathematical background for our proposals. The coefficients of the expansion can be used directly to analyse the spatial characteristics of the magnetic fields at the point of the expansion. To exploit this, we propose a method to shift the reference point of the expansion, which is directly applied to the coefficients. This offers the possibility to obtain the spatial characteristics of the magnetic fields at different positions from one set of coefficients calculated in a measurement-based procedure. For the calculation of the coefficients, we propose to use spherical t-designs as efficient quadrature nodes. An actual measurement at these nodes provides the coefficients of the effective magnetic fields in MPI. Finally, we use the coefficients at different expansion points for the characterisation of static and dynamic fields in MPI.

1.1 Magnetic particle imaging

MPI is a tracer-based imaging modality, which determines the spatial distribution of superparamagnetic iron oxide nanoparticles (SPIOs) using magnetic fields for signal generation and encoding. The spatial encoding field is a static linear field, called selection field, and in the case considered here it has an FFP topology. Only nanoparticles that are located inside a small low-field-region (LFR) around the FFP are unsaturated and able to non-linearly respond to an excitation field, which leads to the spatial encoding of the signal [38]. In our scenario, three orthogonal excitation directions with sinusoidal excitation and rational frequency ratios are chosen, such that the FFP moves on a Lissajous trajectory running through a cuboidal FOV. A detailed mathematical model of the MPI receive signal is described in [26]. Here, we consider the signal under the assumptions of an ideal analogue filter and no feed through. In this case, the magnetisation response of the nanoparticles in the LFR induces the voltage signal into multiple receive coils $k=1,\dots ,K$ given by

(1) \begin{align} u_k(t) = -\mu _0\int _{\Omega } \left \langle \boldsymbol{p}^k_{\textrm {r}}(\boldsymbol{r}) , \frac {\partial \bar {\boldsymbol{m}}}{\partial t}\!\left (\boldsymbol{B}(\boldsymbol{r},t),t\right ) \right \rangle c(\boldsymbol{r}) \textrm {d}\boldsymbol{r}, \end{align}

where $c\,:\,\Omega \,\rightarrow \,\mathbb {R}_{\geq 0}$ describes the particle distribution inside the FOV $\Omega \subseteq \mathbb {R}^3$ , $\bar {\boldsymbol{m}}\,:\,\mathbb {R}^3\times \mathbb {R}\,\rightarrow \,\mathbb {R}^3$ is the mean magnetic moment of the particles, $\boldsymbol{p}^k_{\textrm {r}}\,:\,\mathbb {R}^3\,\rightarrow \,\mathbb {R}^3$ is the coil sensitivity of the receive coils and $\mu _0$ is the vacuum permeability. The magnetic moment of the particles is the response to the applied magnetic field $\boldsymbol{B}\,:\,\mathbb {R}^3\times \mathbb {R}\,\rightarrow \,\mathbb {R}^3$ , which is composed of the static selection and dynamic drive fields.

An exemplary MPI experiment is shown in Figure 1. A mouse is placed in the centre of the scanner bore. During the measurement, a tracer with SPIOs is injected. In a measurement scenario, typically selection fields with gradient strengths between 0.2 and $\text{7 Tm}^{-1}$ [20, 31] in $z$ -direction with half of this value in $x$ - and $y$ -direction as well as drive-field amplitudes of about 6 to 18 mT [20, 54] are used. The drive-field amplitude is limited since higher amplitudes can cause peripheral nerve stimulation [42]. Higher gradient strengths lead to a smaller signal generating LFR such that the resolution of the imaging system increases [38]. However, this comes at the cost of a reduced size of the FOV. For example, using a gradient strength of $\text{2.0 Tm}^{-1}$ with a drive-field amplitude of 12 mT can shift the FFP up to $\pm$ 12 mm in $x$ - and $y$ -direction and $\pm$ 6 mm in $z$ -direction. This yields a FOV of 24 × 24 × 12 mm $^{3}$ , which does not cover larger objects like mice or rats. To this end, a multi-patch approach is used [30]. Additional static magnetic fields, named focus fields, shift the initial FFP such that different patches cover a larger FOV. In Figure 1, an experiment is sketched, where a set of nine different patches is used to cover the mouse.

Figure 1. An MPI measurement is illustrated schematically. MPI scanner and a three-axis robot are controlled by a single computer. Prior to the measurement, a mouse is placed in the centre of the scanner bore using the robot. During the MPI measurement tracer material containing SPIOs is injected into the mouse. As the size of the mouse exceeds the size of a single-patch FOV, multiple patches are used to cover its body. Off-centre patches are warped due to the spatial characteristics of the static and dynamic fields.

Due to field imperfections, the trajectory of each patch is slightly different, which might have multiple negative consequences. If the FFP is not moving along the expected path the spatial encoding changes, which may lead to image artefacts. Moreover, patches might shift to different positions, which may lead to gaps in the sampled FOV like it is illustrated in Figure 1. Lastly, different or spatially dependent drive-field amplitudes can result in incorrect estimations of the tracer concentration.

One main goal of this work is to quantify the severity of the distortions of the underlying magnetic fields. Ideally, only constant and linear fields are present in MPI, which, when represented by the coefficients of a spherical harmonic expansion at the FFP of the selection field, leads to only a few non-zero coefficients as shown in Figure 2. Local imperfections are directly observable in non-zero coefficients of orthogonal field components or coefficients of higher order. For the calculation of the coefficients, we use field measurements at spherical t-design quadrature nodes. The corresponding expansion has the centre of the sphere as its expansion point. However, since the exact position of the FFP is not yet known at the time of measurement, this is likely not the FFP. Additional measurement with the FFP as centre can be avoided by shifting the reference point of the expansion by a linear transformation of the coefficients, which will be introduced in this paper. This can also be exploited to compare the local fields at the centres of the patch positions, which ideally should be identical. This is done by moving the reference point of extension to each centre, respectively.

Figure 2. Spherical harmonic coefficients of two ideal magnetic fields in MPI. On the left, an ideal selection field with gradient strength of 2 $Tm^{-1}$ in $z$ -direction and −1 $Tm^{-1}$ in $x$ - and $y$ -direction is shown. The gradient strength is represented by the linear coefficients ( $l=\textit{1}$ ) of the spherical harmonic expansion of the corresponding field direction. An ideal focus field in $x$ -direction with a 24 mT field strength is visualised on the right. This constant field is represented by the constant coefficient ( $l=\textit{0}$ ) of the expansion in $x$ -direction.

2. Theory

2.1 Unique solution of Laplace’s equation

This section presents the theoretical foundations that are essential for comprehending the concepts that will be discussed subsequently. It does not introduce new ideas. We start with the introduction of solid spherical harmonic expansions as general solution of Laplace’s equation. In order to solve the equation, we use a Dirichlet boundary condition on a sphere, which is a natural choice for solutions expanded with spherical harmonics.

Definition 2.1. Let $f\in \mathcal {C}^2(\mathcal {B}_R( \boldsymbol {\rho } ),\mathbb {R})$ with $\mathcal {B}_R( \boldsymbol {\rho } )\mathrel {\mathop :}=\left \{ \boldsymbol{a}\in \mathbb {R}^3 \,:\, \left \lVert \boldsymbol{a} - \boldsymbol {\rho } \right \rVert _2 \leq R \right \}$ , $ \boldsymbol {\rho } \in \mathbb {R}^3$ , $R\in \mathbb {R}_+$ . Laplace’s equation with Dirichlet boundary condition is given by

(2) \begin{align} \begin{cases} \Delta f(\boldsymbol{a}) = 0\qquad &\forall \boldsymbol{a}\in \overset {\circ }{\mathcal {B}}_R( \boldsymbol {\rho } )\\[3pt] f(\boldsymbol{a}) = \hat {f}(\boldsymbol{a}) \quad &\forall \boldsymbol{a}\in \partial \mathcal {B}_R( \boldsymbol {\rho } ). \end{cases} \end{align}

The boundary condition $\hat {f}\in \mathcal {C}(\partial \mathcal {B}_R( \boldsymbol {\rho } ),\mathbb {R})$ is given on the surface of the ball denoted by $\partial \mathcal {B}_R( \boldsymbol {\rho } )$ , while the interior is denoted by $\overset {\circ }{\mathcal {B}}_R( \boldsymbol {\rho } )$ .

First, we introduce a solution of (2) on the unit ball, which we extend later for an arbitrary radius and centre of the ball.

Proposition 2.2. Let $f\in \mathcal {C}^2(\mathcal {B}_1(\boldsymbol {0}),\mathbb {R})$ and $\hat {f}\,:\,\mathbb {S}^2\,\rightarrow \,\mathbb {R}\in L^2(\mathbb {S}^2)$ fulfill ( 2 ). Then, $f$ can be written as solid spherical harmonic expansion

(3) \begin{align} f(\boldsymbol{a}) = \sum _{l=0}^{\infty }\sum _{m=-l}^l \gamma _{l,m} Z_l^m(\boldsymbol{a})\qquad \forall \boldsymbol{a}\in \mathcal {B}_1(\boldsymbol {0}), \end{align}

where the normalised real solid spherical harmonics $Z_l^m$ as an extension of normalised real spherical harmonics are defined by

\begin{align*} Z_l^m \,:\, \mathbb {R}^3\,\rightarrow \,\mathbb {R}, (r,\vartheta ,\varphi ) \mapsto K_l^{\left \lvert m \right \rvert } r^l P_l^{\left \lvert m \right \rvert }(\cos \vartheta )\cdot \begin{cases} \sqrt {2}\cos (m\varphi ) & m \gt 0\\ \sqrt {2}\sin (\left \lvert m \right \rvert \varphi ) & m \lt 0\\ 1 & m = 0 \end{cases}, \end{align*}

with $ K_l^m = \sqrt {\frac {(l-m)!}{(l+m)!}}$ and the associated Legendre polynomials $P_l^m$ [4]. The solid spherical coefficients $\gamma _{l,m}\in \mathbb {R}$ of the expansion can be calculated by the orthogonal projection

\begin{align*} \gamma _{l,m} = \frac {2l+1}{4\pi }\int _{\mathbb {S}^2}\hat {f}(\boldsymbol{a})Z_l^m(\boldsymbol{a}) \textrm {d}\boldsymbol{a}. \end{align*}

Proof. The proposition holds since the restricted solid spherical harmonics $Z_l^m|_{\mathbb {S}^2}$ form an orthogonal basis of $L^2(\mathbb {S}^2)$ and thus the coefficients can be calculated by the orthogonal projection.

Remark 1. Each solid spherical harmonic expansion $\sum \limits _{l=0}^\infty \sum \limits _{m=-l}^l \gamma _{l,m} Z_l^m$ satisfies (2) [4].

Next, we generalise Proposition 2.2 for arbitrary radius $R$ and centre $ \boldsymbol {\rho }$ of the ball $\mathcal {B}_R( \boldsymbol {\rho } )$ . In this case, $ \boldsymbol {\rho }$ determines the centre of the series expansion. Thus, we obtain an expansion $f^{ \boldsymbol {\rho } }$ , where the centre of the expansion is denoted by the superscript. Analogously, we denote the domain of the expansion by $\mathcal {B}_R^{ \boldsymbol {\rho } }(\boldsymbol {0})$ in its local coordinate system, which corresponds to $\mathcal {B}_R( \boldsymbol {\rho } )$ in the global coordinate system as it is visualised in Figure 3.

Proposition 2.3. Let $f\in \mathcal {C}^2(\mathcal {B}_R( \boldsymbol {\rho } ),\mathbb {R})$ and $\hat {f}\,:\,\partial \mathcal {B}_R( \boldsymbol {\rho } )\,\rightarrow \,\mathbb {R}\in L^2(\partial \mathcal {B}_R( \boldsymbol {\rho } ))$ fulfill equation ( 2 ) for arbitrary $R\in \mathbb {R}_+$ and $ \boldsymbol {\rho } \in \mathbb {R}^3$ . The coefficients $\gamma _{l,m}( \boldsymbol {\rho } , R)$ depending on $ \boldsymbol {\rho }$ and $R$ can be calculated by

(4) \begin{align} \gamma _{l,m}( \boldsymbol {\rho } , R) = \frac {2l+1}{4\pi } \int _{\mathbb {S}^2} \hat {f}{\left (R\boldsymbol{a} + \boldsymbol {\rho } \right )}Z_l^m(\boldsymbol{a}) \textrm {d} \boldsymbol{a}. \end{align}

With $ \gamma _{l,m}( \boldsymbol {\rho } ) = \frac {1}{R^l}\gamma _{l,m}( \boldsymbol {\rho } ,R)$ , the solid harmonic expansion of $f$ can be formulated as

(5) \begin{align} f^{ \boldsymbol {\rho } }(\boldsymbol{a}) = \sum _{l=0}^\infty \sum _{m=-l}^l \gamma _{l,m}( \boldsymbol {\rho } ) Z_l^m(\boldsymbol{a}) \qquad \forall \boldsymbol{a}\in \mathcal {B}_R^{ \boldsymbol {\rho } }(\boldsymbol {0}), \end{align}

with $f^{ \boldsymbol {\rho } } \,:\!=\, \tau _{ \boldsymbol {\rho } }(f)$ where $\tau _{ \boldsymbol {\rho } }(f)\,:\,\mathcal {B}_R(\boldsymbol {0})\,\rightarrow \,\mathbb {R},\boldsymbol{a} \mapsto f(\boldsymbol{a} + \boldsymbol {\rho } )$ denotes the shift operator on $\mathcal {C}^2$ functions. The centre $ \boldsymbol {\rho }$ determines the origin of the underlying coordinate system, which is denoted by a superscript $ \boldsymbol {\rho }$ for $f$ and its domain.

Proof. Let $f^{ \boldsymbol {\rho } ,R} \,:\!=\, \tau _{ \boldsymbol {\rho } }(\sigma _R(f))\,:\,\mathcal {B}_1(\boldsymbol {0})\,\rightarrow \,\mathbb {R},\boldsymbol{a} \mapsto f(R\boldsymbol{a} + \boldsymbol {\rho } )$ and $\hat {f}^{ \boldsymbol {\rho } ,R} \,:\!=\, \tau _{ \boldsymbol {\rho } }(\sigma _R(\hat {f}))$ analogous, where $\sigma _R$ denotes the scaling operator on $\mathcal {C}^2$ functions. Using Proposition 2.2, the coefficients of $f^{ \boldsymbol {\rho } ,R}$ can be calculated by

\begin{align*} \gamma _{l,m}( \boldsymbol {\rho } , R) = \frac {2l+1}{4\pi } \int _{\mathbb {S}^2} \hat {f}^{ \boldsymbol {\rho } ,R}{\left (\boldsymbol{a}\right )}Z_l^m(\boldsymbol{a}) \textrm {d} \boldsymbol{a}, \end{align*}

which yields (4). Thus, we get

\begin{align*} f^{ \boldsymbol {\rho } ,R} = \tau _{ \boldsymbol {\rho } }(\sigma _R(f)) &= \sum _{l=0}^{\infty }\sum _{m=-l}^l \gamma _{l,m}( \boldsymbol {\rho } ,R) Z_l^m\\ \Leftrightarrow \, \tau _{ \boldsymbol {\rho } }(f) = \tau _{ \boldsymbol {\rho } }(\sigma _R(\sigma ^{-1}_R(f))) &= \sum _{l=0}^{\infty }\sum _{m=-l}^l \gamma _{l,m}( \boldsymbol {\rho } ,R)\,\sigma ^{-1}_R(Z_l^m) \\ &= \sum _{l=0}^{\infty }\sum _{m=-l}^l \underbrace {\frac {\gamma _{l,m}( \boldsymbol {\rho } ,R)}{R^l}}_{=:\,\gamma _{l,m}( \boldsymbol {\rho } )} Z_l^m, \end{align*}

since $Z_l^m$ are homogeneous polynomials so that $\sigma ^{-1}_R(Z_l^m(\boldsymbol{a})) = Z_l^m\!{\left (\frac {1}{R}\boldsymbol{a}\right )} = \frac {1}{R^l}Z_l^m(\boldsymbol{a})$ . With $f^{ \boldsymbol {\rho } } \,:\!=\, \tau _{ \boldsymbol {\rho } }(f)$ equation (5) follows.

Note 1. For a better readability, the indices $R$ and $ \boldsymbol {\rho }$ are omitted if $R=1$ and $ \boldsymbol {\rho } = \boldsymbol {0}$ .

2.2 Translation

Figure 3. Different coordinate systems of the coefficients with the domain of the function $f$ . The black coordinate system represents the initial coordinate system of $f^{ \boldsymbol {\rho } }$ at the expansion point $ \boldsymbol {\rho }$ . Using a shift $\boldsymbol{v}$ , the coefficients of $f^{\boldsymbol{q}}$ depend on the shifted blue coordinate system with its origin at $\boldsymbol{q} = \boldsymbol {\rho } + \boldsymbol{v}$ . Both local coordinate systems, $\mathcal {B}_R^{ \boldsymbol {\rho } }(\boldsymbol {0})$ and $\mathcal {B}_R^{\boldsymbol{q}}({-}\boldsymbol{v})$ , are equal to $\mathcal {B}_R( \boldsymbol {\rho } )$ in the global coordinate system (red).

The coefficients $\gamma _{l,m}( \boldsymbol {\rho } )$ correspond to a solid harmonic expansion around the centre $ \boldsymbol {\rho }$ of the domain of the boundary condition. Next, we introduce a translation operator $\hat {\tau }$ that allows to transform them, such that they correspond to a solid harmonic expansion around a new centre point $\boldsymbol{q} \in \mathcal {B}_R( \boldsymbol {\rho } )$ .

Steinborn and Ruedenberg [45] introduced an addition theorem for complex solid spherical harmonics, which states how these functions map under translation. In particular, they state how to calculate $\tilde {Z}_l^m(\boldsymbol{a}_1 + \boldsymbol{a}_2)$ with the sum $\sum _{\lambda ,\mu } \tilde {Z}_{l-\lambda }^{m-\mu }(\boldsymbol{a}_1)\tilde {Z}_{\lambda }^{\mu }(\boldsymbol{a}_2)$ , where $\tilde {Z}$ denotes the complex modified regular solid harmonics defined in their paper and $\boldsymbol{a}_1,\boldsymbol{a}_2 \in \mathbb {R}^3$ . Rico et al. [41] transferred this addition theorem to real solid harmonics. With slight adaptions, this addition theorem can be applied to the normalised real solid harmonics used in this paper.

However, the addition theorem changes the basis functions of the expansion, while we are interested in the coefficients of the shifted expansion. To this end, we transfer the adapted addition theorem from the solid harmonics to the coefficients to obtain a mapping from the coefficients of one expansion point to the coefficients of a new expansion point. Consequently, any shift of the coordinate system is reflected in a transformation of the coefficients of the expansion, rather than the basis functions.

In order to establish the mapping of the coefficients, a linear operator is introduced which describes a truncated solid harmonic expansion. Henceforth, we assume $f\,:\,\mathcal {B}_R( \boldsymbol {\rho } )\,\rightarrow \,\mathbb {R} \in \mathbb {P}^L$ and $\hat {f}\,:\,\partial \mathcal {B}_R( \boldsymbol {\rho } )\,\rightarrow \,\mathbb {R}\in L^2(\partial \mathcal {B}_R( \boldsymbol {\rho } ))$ fulfill (2) for $R\in \mathbb {R}_+$ and $ \boldsymbol {\rho } \in \mathbb {R}^3$ , with $\mathbb {P}^L$ being the space of all harmonic polynomials up to degree $L\in {\mathbb {N}}$ .

Definition 2.4. We define the truncated solid harmonic expansion as a linear operator

\begin{align*} \mathcal {S}_L \, :\, \mathbb {R}^{(L+1)^2} &\rightarrow \mathbb {P}^L\\ \boldsymbol {\gamma } (\boldsymbol{q}) &\mapsto \left (\boldsymbol{a} \mapsto \sum _{l=0}^L\sum _{m=-l}^l \gamma _{l,m}(\boldsymbol{q}) Z_l^m(\boldsymbol{a})\right ), \end{align*}

where $ \boldsymbol {\gamma } (\boldsymbol{q}) = \left (\gamma _{l,m}(\boldsymbol{q})\right )_{\substack {l = 0,\dots ,L \\ m = -l,\dots ,l}} \in \mathbb {R}^{(L+1)^2}$ is a vector containing all coefficients up to $l=L$ at expansion centre $\boldsymbol{q}\in \mathcal {B}_R( \boldsymbol {\rho } )$ .

Remark 2.

1. (i) Since we assume that $f$ is a polynomial of at most degree $L$ , (5) is equivalent to

   \begin{align*} f^{ \boldsymbol {\rho } } = \mathcal {S}_L( \boldsymbol {\gamma } ( \boldsymbol {\rho } )) \end{align*}
   for $ \boldsymbol {\gamma } ( \boldsymbol {\rho } )$ calculated with (4).

2. (ii) A translation of the coordinate system by a shift $\boldsymbol{v} = \boldsymbol{q} - \boldsymbol {\rho }$ to a new centre $\boldsymbol{q}\in \mathcal {B}_R( \boldsymbol {\rho } )$ of the series expansion for $\boldsymbol{a}\in \mathcal {B}_R^{\boldsymbol{q}}({-}\boldsymbol{v})$ is described by

   \begin{align*} f^{\boldsymbol{q}}(\boldsymbol{a}) = \tau _{\boldsymbol{v}}\big (f^{ \boldsymbol {\rho } }(\boldsymbol{a})\big ) &= \tau _{\boldsymbol{v}} \big (\mathcal {S}_L( \boldsymbol {\gamma } ( \boldsymbol {\rho } ))(\boldsymbol{a})\big )\\ &= \sum _{l=0}^L\sum _{m=-l}^l \gamma _{l,m}( \boldsymbol {\rho } ) \,\tau _{\boldsymbol{v}}\big (Z_l^m(\boldsymbol{a})\big ). \end{align*}

The translation of the solid spherical harmonics can be transferred to the coefficients as it is stated in the following theorem.

Theorem 2.5. For any $\boldsymbol{q} \in \mathcal {B}_R( \boldsymbol {\rho } )$ , an operator $\hat {\tau }_{\boldsymbol{v}}\,:\,\mathbb {R}^{(L+1)^2}\,\rightarrow \,\mathbb {R}^{(L+1)^2}$ exists with $\boldsymbol{v} = \boldsymbol{q} - \boldsymbol {\rho }$ such that

commutes, that is,

\begin{align*} \tau _{\boldsymbol{v}} \circ \mathcal {S}_L( \boldsymbol {\gamma } ( \boldsymbol {\rho } )) = \mathcal {S}_L \circ \hat {\tau }_{\boldsymbol{v}}\big ( \boldsymbol {\gamma } ( \boldsymbol {\rho } )\big ). \end{align*}

The actual calculation of $\hat {\tau }_{\boldsymbol{v}}\big ( \boldsymbol {\gamma } ( \boldsymbol {\rho } )\big ) = \boldsymbol {\gamma } (\boldsymbol{q})$ is given by equations ( 22 ) to ( 25 ) in Appendix A.2 .

Proof. The proof of the theorem is given in the appendix. First, we adapt the addition theorem for unnormalised real solid spherical harmonics provided by Rico et al. in their work [41] to our normalised ones in Appendix A.1. Applying the addition theorem to the solid harmonic expansion and reordering of the sums leads to the addition theorem for the solid coefficients and with that to the proof of the theorem as it is shown in the second Appendix A.2.

The domain of the expansion depends on the boundary condition used for the calculation of the coefficients, so the domain of the shifted harmonic polynomial $\mathcal {S}^L( \boldsymbol {\gamma } (\boldsymbol{q}))$ is given by $\mathcal {B}_R( \boldsymbol {\rho } - \boldsymbol{q})$ , respectively, $\mathcal {B}_R^{\boldsymbol{q}}({-}\boldsymbol{v})$ in its local coordinate system. In Figure 3, the domain of the harmonic polynomial is shown with the coordinate systems centred at the original and shifted expansion point.

Remark 3.

1. (i) The effort for calculating $(L+1)^2$ coefficients $\hat {\tau }_{\boldsymbol{v}}\big ( \boldsymbol {\gamma } ( \boldsymbol {\rho } )\big )$ with (22) to (25) is $\mathcal {O}{\left ( L^4 \right )}$ .

2. (ii) Theorem 2.5 can be generalised for a shift between arbitrary points $\boldsymbol{p}, \boldsymbol{q} \in \mathcal {B}_R( \boldsymbol {\rho } )$ .

2.3 Efficient quadrature

In order to obtain the solid coefficients for a polynomial $f\,:\,\mathcal {B}_R( \boldsymbol {\rho } )\,\rightarrow \,\mathbb {R}$ of degree $L\in {\mathbb {N}}_0$ its values on the boundary $\partial \mathcal {B}_R( \boldsymbol {\rho } )$ have to be known. For instance in the magnetic field determination application scenario, these values can be measured. In MPI, for example, the choice of the measurement points is only restricted by the size and shape of the scanner bore and so far a Gauss–Legendre quadrature was used in MPI for the calculation of the coefficients [11]. A more efficient way to choose the measurement points are spherical t-designs [14, 5], which are introduced next.

Definition 2.6. A spherical t-design is a set of nodes $\left \{ \boldsymbol{a}_k \right \}_{k=1,\dots ,N}\subseteq \mathbb {S}^2$ such that

\begin{align*} \int _{\mathbb {S}^2} \mathcal {Y}(\boldsymbol{a}) \textrm {d}\boldsymbol{a} = \frac {4\pi }{N}\sum _{k=1}^{N} \mathcal {Y}(\boldsymbol{a}_k)\quad \forall \mathcal {Y} \in \mathbb {P}^t_{\mathbb {S}}, \end{align*}

with $\mathbb {P}^t_{\mathbb {S}} = \text {span}{\left \{ Z_l^m|_{\mathbb {S}^2} \,:\, l \leq t \right \}}$ being the set of all polynomials up to degree $t\in {\mathbb {N}}_0$ on the unit sphere [5].

Remark 4. The spherical t-design is a very efficient sampling pattern for the quadrature on a spherical surface. More information on the calculation of these t-designs is given in [23], while explicit definitions of various spherical t-designs can be found in [22]. For instance, the smallest known $8$ -design only consists of $36$ points [23]. In comparison, $45$ Gauss–Legendre quadrature nodes are required for the same accuracy [52].

Proposition 2.7. Assume $f\,:\,\mathcal {B}_R( \boldsymbol {\rho } )\,\rightarrow \,\mathbb {R}$ to be a polynomial of degree $L \in {\mathbb {N}}_0$ fulfilling ( 2 ) with $\hat {f}\,:\,\partial \mathcal {B}_R( \boldsymbol {\rho } )\,\rightarrow \,\mathbb {R}\in L^2(\partial \mathcal {B}_R( \boldsymbol {\rho } ))$ . Let $\left \{ \boldsymbol{a}_k \right \}_{k=1,\dots ,N}$ be a $2L$ -design. For $l\leq L$ , it holds that

\begin{align*} \gamma _{l,m}( \boldsymbol {\rho } ) = \frac {2l+1}{N R^l}\sum _{k=1}^N \hat {f}(R\boldsymbol{a}_k+ \boldsymbol {\rho } )Z_l^m(\boldsymbol{a}_k). \end{align*}

Proof. Let $l\leq L$ . Then, it holds that the degree of the product $\hat {f} Z_l^m$ is at most $2L$ . With Proposition 2.3 and Definition 2.6, we get

\begin{align*} \gamma _{l,m}( \boldsymbol {\rho } ) &= \frac {1}{R^l}\frac {2l+1}{4\pi } \int _{\mathbb {S}^2}\hat {f}(R\boldsymbol{a}+ \boldsymbol {\rho } )Z_l^m(\boldsymbol{a}) \textrm {d}\boldsymbol{a} \\ &= \frac {1}{R^l}\frac {2l+1}{4\pi }\frac {4\pi }{N}\sum _{k=1}^N \hat {f}(R\boldsymbol{a}_k+ \boldsymbol {\rho } )Z_l^m(\boldsymbol{a}_k). \end{align*}

Since $\deg (f) \leq L$ , the solid harmonic expansion (5) can be truncated at $l = L$ .

3. Methods

The requisite tools are now available to efficiently measure and analyse the magnetic fields in MPI using spherical harmonic coefficients. We will begin by demonstrating that each component of the magnetic fields satisfies Laplace’s equation. Once this condition has been met, Proposition 2.7 can be applied to obtain solid harmonic expansions representing each component of the fields.

In order to assess the coefficients of the magnetic fields in MPI, we introduce the different fields with their ideal coefficients afterwards. One of the presented fields features a characteristic unique point, which we exploit for a unique representation of the fields by shifting the coefficients into this point.

3.1 Magnetic fields

Magnetic fields in MPI are quasi-static, generated by electric currents outside the scanner bore. As these fields fulfill Laplace’s equation, we are able to apply Proposition 2.3 and expand the field as a solid harmonic expansion.

Lemma 3.1. Let $\boldsymbol{B} = (B_x,B_y,B_z) \in \mathcal {C}^2(\Omega ,\mathbb {R}^3)$ be a quasi-static magnetic field, where $\Omega \subseteq \mathbb {R}^3$ describes the region where no electric current flows. Then, the components $B_i$ for $i=x,y,z$ fulfill Laplace’s equation for all $\boldsymbol{a} \in \Omega$ .

Proof. The static magnetic field fulfills Maxwell’s equations

\begin{align*} \nabla \cdot \boldsymbol{B}(\boldsymbol{a}) &= 0\\ \nabla \times \boldsymbol{B}(\boldsymbol{a}) &= \boldsymbol {0} ,\end{align*}

for all $\boldsymbol{a}\in \Omega$ [24]. The second equation is equal to zero since no electric current flows in $\Omega$ . Thus, it holds that $\boldsymbol{B} = -\nabla \Phi$ where $\Phi \,:\,\Omega \,\rightarrow \,\mathbb {R}$ is the magnetic scalar potential [52, 24]. In [52], it is shown that $-\Delta \Phi = 0$ follows, which implies that $\Delta B_i = 0$ .

Note 2. In the following, the index $i$ denotes the $x$ -, $y$ - and $z$ -component of the magnetic field $\boldsymbol{B}$ .

Proposition 3.2. Assume that the magnetic field $\boldsymbol{B}\,:\,\Omega \,\rightarrow \,\mathbb {R}^3\in \mathbb {P}^L$ can be described by a harmonic polynomial of degree $L$ . Let $R\in \mathbb {R}_+$ and $ \boldsymbol {\rho } \in \mathbb {R}^3$ such that $\mathcal {B}_R( \boldsymbol {\rho } )\subseteq \Omega$ and let $\left \{ \boldsymbol{a}_k \right \}_{k=1,\dots ,N}$ be a $2L$ -design. With a given boundary condition $B_i(\boldsymbol{a}) = \hat {B}_i(\boldsymbol{a}) \, \forall \boldsymbol{a}\in \partial \mathcal {B}_R( \boldsymbol {\rho } )$ with $\hat {B}_i\in L^2{\left (\partial \mathcal {B}_R( \boldsymbol {\rho } )\right )}$ the field can be formulated as a solid harmonic expansion

\begin{align*} \boldsymbol{B}^{ \boldsymbol {\rho } }(\boldsymbol{a}) = \sum _{l=0}^L\sum _{m=-l}^l \boldsymbol {\gamma\!}_{l,m}( \boldsymbol {\rho } ) Z_l^m(\boldsymbol{a}) \qquad \forall \boldsymbol{a}\in \mathcal {B}_R^{ \boldsymbol {\rho } }(\boldsymbol {0}), \end{align*}

with coefficients $ \boldsymbol {\gamma\!}_{l,m} = (\gamma _{l,m}^x, \gamma _{l,m}^y, \gamma _{l,m}^z)$ calculated by

\begin{align*} \gamma _{l,m}^i( \boldsymbol {\rho } ) = \frac {1}{R^l}\frac {2l+1}{N} \sum _{k=1}^N \hat {B_i}{\left (R\boldsymbol{a}_k + \boldsymbol {\rho } \right )}Z_l^m(\boldsymbol{a}_k). \end{align*}

Proof. Due to Lemma 3.1, the magnetic field $\boldsymbol{B}$ and the boundary $\hat {\boldsymbol{B}}$ fulfill all assumptions of Proposition 2.7, which implies both equations.

The expansion of the magnetic field $\boldsymbol{B}$ at the expansion point $ \boldsymbol {\rho }$ with coefficients $ \boldsymbol {\gamma } ^x( \boldsymbol {\rho } )$ , $ \boldsymbol {\gamma } ^y( \boldsymbol {\rho } )$ and $ \boldsymbol {\gamma } ^z( \boldsymbol {\rho } )$ characterises the magnetic field locally in $x$ -, $y$ - and $z$ -direction, respectively. Similar to Taylor series, the expansion can be written as polynomial where the polynomial degree increases with the index $l$ . Thereby, $ \boldsymbol{\gamma\!}_{0,0}$ describes the constant part of the magnetic field, while the coefficients $ \boldsymbol {\gamma\!}_{1,m}$ contain the information about its linear behaviour, $ \boldsymbol {\gamma\!}_{1,-1}$ describes the behaviour in $y$ -direction, $ \boldsymbol {\gamma\!}_{1,0}$ in $z$ -direction and $ \boldsymbol {\gamma\!}_{1,1}$ in $x$ -direction. The coefficients for $l\gt 1$ characterise the non-linear behaviour of the magnetic field.

3.2 Magnetic fields in MPI

In MPI, two main magnetic fields are used for signal encoding and generation: a linear selection field $\boldsymbol{B}_{\textrm {SF}}\,:\,\mathbb {R}^3\,\rightarrow \,\mathbb {R}^3$ and dynamic drive fields $\boldsymbol{B}_{\textrm {DF}}^i\,:\,\mathbb {R}^3\times \mathbb {R}\,\rightarrow \,\mathbb {R}^3$ . Each dynamic drive field can be separated into the constant coil sensitivity $\boldsymbol{p}_{\textrm {DF}}^i\,:\,\mathbb {R}^3\,\rightarrow \,\mathbb {R}^3$ and the sinusoidal current $I^i\,:\,\mathbb {R}\,\rightarrow \,\mathbb {R}$ such that $\boldsymbol{B}_{\textrm {DF}}^i(\boldsymbol{r},t) = I^i(t)\boldsymbol{p}_{\textrm {DF}}^i(\boldsymbol{r})$ . In our set-up, the three orthogonal drive field coils are also the receive coils so that $\boldsymbol{p}_{\textrm {DF}}^k = \boldsymbol{p}^k_{\textrm {r}}$ with $k\in \left \{ 1,2,3 \right \}$ in (1). In the multi-patch setting described in the problem statement, additional patch-wise constant focus fields $\boldsymbol{B}_{\textrm {FF}}^i\,:\,\mathbb {R}^3\,\rightarrow \,\mathbb {R}^3$ are applied to obtain a larger FOV. More general information on the imaging principles of MPI and the set-up of an MPI scanner can be found in [27, 28], while the mathematical background of MPI is described in [26].

The solid coefficients at the expansion point $\boldsymbol {0}$ of ideal selection, focus and drive fields are listed in Table 1. All coefficients not mentioned in the table are zero. Since the selection field is a linear field, only the coefficients for $l=1$ are non-zero and describe the gradient strength of $g\in [0\;\text{Tm}^{-1},2.5\;\text{Tm}^{-1}]$ in $z$ -direction and $-0.5g$ in $x$ - and $y$ -direction. Meanwhile, only the constant coefficients for $l=0$ of the ideal drive and focus field are non-zero. The focus fields are characterised by the shift $f^x,f^y\in [{-}17\; \textrm {mT}, 17\;\textrm {mT}]$ and $f^z\in [{-}42\; \textrm{mT}, 42\; \textrm{mT}]$ , while the drive field is characterised by its amplitude $d^i\in [0\; \textrm{mT},14\; \textrm{mT}]$ with $d^i = \max _t(I^i(t)) \boldsymbol{p}_{\textrm {DF}}^i(\boldsymbol {0})$ .

Table 1. Coefficients of the three different ideal magnetic fields in MPI in Tm $^{-l}$

3.3 Measurement set-up

Measuring magnetic fields can be done using different devices like Hall-effect sensors, SQUID sensors or induction sensors [50]. Gaussmeters with a three-axis Hall sensor are very accurate and therefore widely used for magnetic field measurements [39, 40]. Hence, we use a 3-channel gaussmeter with a three-axis high-sensitivity Hall-effect sensor from Lake Shore (model 460, Westerville, USA) [13] for the measurement of the static fields of the MPI scanner. Its accuracy is the sum of the reading error of about $\pm 0.10\%$ and $\pm 0.005\%$ of the chosen range [13]. The used range of $\pm 0.3$ T results in a maximum reading error of $\pm 300\,\mu$ T and a range error of $\pm 15\,\mu$ T. The coil sensitivity of the dynamic drive field is measured with a three-axis coil sensor, which is connected to an analogue-to-digital converter (ADC) for a digitisation of the induced voltage signal [48]. Each coil has a radius of 2.5 mm and an accuracy of about $\pm 8\%$ . The magnetic fields are measured in our preclinical MPI system 25/20 FF (Bruker BioSpin MRI GmbH, Ettlingen, Germany), which is equipped with a three-axis Cartesian robot (isel Germany AG, Eichenzell, Germany) for an easy and accurate positioning of the measurement devices. The robot has a repetition accuracy in each direction of $\pm 0.02$ mm and an angle error of $\pm 5\%$ for a motor step angle of $1.8^{\circ }$ . Taking the accuracy of the gaussmeter or the coil sensor into account measurement errors due to mislocation of the robot are small and can be neglected. All fields are measured at the $36$ points of a spherical $8$ -design [22], which are approached by the robot. Figure 4 shows the measurement set-up including the chosen spherical $8$ -design. According to the 12 cm diameter of the scanner bore, the points are rescaled to obtain a sphere with radius 42 mm. The spheres were chosen as large as possible while keeping a safety margin that prevents collision of probe and scanner bore. The centre $ \boldsymbol {\rho }$ is chosen near the FFP of the selection field. At every point, all three field directions are measured simultaneously.

For each of the magnetic fields, we perform multiple measurements with different field strengths. Note, that the following field values describe the input parameters at the scanner, which ideally should result in the ideal magnetic fields described in the last section. First, the selection field at 10 different gradient strengths of 0.25 to $\text{2.5 Tm}^{-1}$ with a step size of $\text{0.25 Tm}^{-1}$ is measured with the Hall sensor. The focus fields are measured with the Hall sensor as well with $11$ different shifts of −17 to 17 mT and a step size of 3.4 mT for the fields in $x$ - and $y$ -direction and −42 to 42 mT with a step size of 8.4 mT in $z$ -direction. For the drive field measured with the coil sensor, four different amplitudes of 6 to 12 mT with a step size of 2 mT for all three directions is set, as this range is most commonly used in our experiments.

Figure 4. Measurement set-up. The Hall-effect sensor of the gaussmeter is mounted on a three-axis Cartesian robot, which moves it to the chosen spherical t-design positions inside the scanner bore. Here, the positions of a spherical 8-design are marked in blue, where the lighter blue indicates the positions with a negative sign in $y$ -direction. The voltage sensor of the gaussmeter transfers the measured data to the computer, which controls the robot movements and the settings of the MPI scanner.

3.4 Unique representation

While the single coils of the coil sensor are all centred around the same point, the three orthogonal detectors inside the Hall sensor are slightly shifted away from the centre of the rod. Therefore, each detector measures the field on a slightly shifted sphere as it is shown in Figure 5. All detectors are located 1.8 mm behind the tip of the rod, and the $x$ - and $y$ -detector are additionally shifted by 2.08 mm outward from the centre. This results in three expansions at slightly different expansion points $ \boldsymbol {\rho } _i$ for each direction. Using the translation from Theorem2.5, the coefficients can be shifted into a common coordinate system centred at the tip of the rod $ \boldsymbol {\rho }$ .

Figure 5. The three-axis Hall-effect sensor has three individual sensors for $x$ -, $y$ - and $z$ -direction, respectively. The sensors for $x$ - and $y$ -direction ( $ \boldsymbol {\rho } _x$ and $ \boldsymbol {\rho } _y$ ) are shifted off centre inside the sensor rod (grey square). For each sensor, the corresponding sphere $\mathcal {B}_R( \boldsymbol {\rho } _x)$ and $\mathcal {B}_R( \boldsymbol {\rho } _y)$ on which the magnetic field is measured are shown, as well as the sphere $\mathcal {B}_R( \boldsymbol {\rho } )$ at the tip of the rod. Additionally, the spatial coordinate system of the MPI scanner is shown on the bottom right. Above, the magnetic field coordinate system is displayed as it is given by the detector orientation of the sensor.

In MPI, the main selection field has a unique scanner-specific FFP. It can be determined from the expansions using root-finding methods, such as Newton’s method used in this work. So far, the expansions of the selection, drive and focus fields have a centre that depends on the measurement set-up. That is, the centre of the spheres on which the spherical t-design positions are located. However, by moving all expansions into the FFP, we remove this dependency and end up with a unique expansion of the magnetic fields that depends exclusively on the specific MPI scanner.

3.5 Implementation

All numerical methods described so far are implemented in the programming language Julia (version 1.8) [6] in the open-source software package SphericalHarmonicExpansions.jl (version 0.1) [44]. The package provides methods for storage and handling of the coefficients of spherical or solid expansions, an efficient quadrature based on spherical t-designs to calculate spherical or solid coefficients, a method to translate coefficients to a different expansion point and methods for fast numerical evaluation of the expansions in Cartesian coordinates. Furthermore, a collection of spherical t-designs can be obtained via the MPIFiles.jl package (version 0.12) [29]. The measurements described in Section 3.3 are performed using MPIMeasurements.jl (version 0.3) [21]. An example script, which shows how to obtain the expansion of a $\text{2 Tm}^{-1}$ selection field from the measurements described above, is provided at https://github.com/IBIResearch/SphericalHarmonicExpansionOfMagneticFields (version 1.0).

3.6 Error analysis

The measurement instruments always have a tiny statistical independent measurement error $\delta \gt 0$ , which is propagated to the coefficients and the final magnetic field using the laws of error propagation [7, 16]. Our measurement set-up has two sources of error: uncertainty in magnetic field measurements and errors in positioning the Hall-effect sensor. A systematic positional error may result from a non-ideal robot mount of the measurement sensors or bending of the rods to which those are attached. Such a systematic error effectively results in a shift of the coordinate system, which is compensated by the shift into the FFP and hence can be neglected. Only non-systematic mislocations need to be accounted for, which are so small in our set-up that they are neglected.

Corollary 3.3. Let $ \boldsymbol {\gamma } ^i( \boldsymbol {\rho } )$ be the solid coefficients calculated with Proposition 2.7 and $\boldsymbol{B}^{ \boldsymbol {\rho } }$ the resulting magnetic field calculated with Proposition 3.2 . Additionally, we have the independent observational errors $\delta _k\!\left (R\boldsymbol{a}_k + \boldsymbol {\rho } \right )$ of the measured boundary condition $\hat {B_i}{\left (R\boldsymbol{a}_k + \boldsymbol {\rho } \right )}$ .

1. (i) The standard deviation for the propagated error of the coefficient $\gamma _{l,m}^i( \boldsymbol {\rho } )$ can be obtained by

   \begin{align*} \varepsilon _{l,m}^i( \boldsymbol {\rho } ) = \frac {2l+1}{N R^l} \sqrt {\sum _{k=1}^N \left (\delta _k\!\left (R\boldsymbol{a}_k + \boldsymbol {\rho } \right )Z_l^m(\boldsymbol{a}_k)\right )^2}. \end{align*}

2. (ii) For each component of the magnetic field $B_i^{ \boldsymbol {\rho } }$ , the standard deviation for the propagated error at a position $\boldsymbol{a}\in \mathcal {B}_R( \boldsymbol {\rho } )$ can be calculated by

   \begin{align*} \hat {\varepsilon }_i^{ \boldsymbol {\rho } }(\boldsymbol{a}) = \sqrt { \sum _{k=1}^N \left (\delta _k(R\boldsymbol{a}_k + \boldsymbol {\rho } ) \; \mathcal {S}_L\!\left (\frac {2l+1}{NR^l} Z_l^m(R\boldsymbol{a}_k+ \boldsymbol {\rho } )\right )\!(\boldsymbol{a})\right )^{\negthickspace 2}}. \end{align*}

Remark 5. Propagating the error of the coefficients through the translation mapping can be done analogously since the translation mapping is linear in the coefficients.

For error analysis, we compare the field values provided by the truncated expansion to the measured ones. In our example study, we use the field measurements obtained at the spherical t-design positions which we also used to create the truncated expansion. In a typical application scenario, independent measurements should be used for error estimation. For the selection field, this is done by

(6) \begin{align} \zeta _i(k) = \frac {B_i^{ \boldsymbol {\rho } }(R \boldsymbol{a}_k) - \hat {B}_i(R\boldsymbol{a}_k + \boldsymbol {\rho } )}{R g_i}, \end{align}

which compares the measured field $\hat {B}_i$ at a scaled t-design position $R\boldsymbol{a}_k + \boldsymbol {\rho }$ with the calculated field $B_i^{ \boldsymbol {\rho } }$ at the same position. The difference is normalised to the scaled gradient strength $g_i$ of the considered direction $i$ . The propagated error undergoes the same normalisation $\hat {\zeta }_i(k) = \frac {\hat {\varepsilon }_i^{ \boldsymbol {\rho } }(R\boldsymbol{a}_k+ \boldsymbol {\rho } )}{R g_i}$ , which allows to assess the approximation quality of the truncated expansions.

4. Results

Exemplary, we examine the results of a $\text{2 Tm}^{-1}$ selection field. In the left part of Table 2, the initial coefficients calculated from the measurement without any post-processing are listed, while in the right part, the processed coefficients are listed. Post-processing consists of normalising with the radius of the measured sphere, correcting the shifts of the Hall sensors and shifting into the FFP calculated using the initial coefficients. Due to the shift into the FFP, the coefficients for $[l,m] = [0,0]$ are equal to zero up to floating point precision. Now, the gradient, which slightly deviates from the ideal gradient (cf. Table 1), can be read directly from the coefficients for $l=1$ . All other non-zero coefficients can be attributed to either measurement error or imperfections of the selection field.

Table 2. Comparison of the initial coefficients (in T, left) and the normalised processed coefficients (in $Tm^{-l}$ , right) of a $2\; Tm^{-1}$ selection field

Figure 6. Solid harmonic analysis of the static fields in MPI, that is, the selection and focus fields. The first row shows a selection field with a gradient strength of $\text{2 Tm}^{-1}$ . In the first column, the coefficients at the FFP of the selection field $ \boldsymbol {\xi } _{\textrm {c}}$ are shown, while in the second column the coefficients at another point $ \boldsymbol {\xi } _{\textrm {s}}$ are shown. Both points are marked in the field plot on the right. In the second row, a focus field of −24 mT in $x$ - and 24 mT in $z$ -direction at both positions is shown. This additional field is required to shift the FFP from $ \boldsymbol {\xi } _{\textrm {c}}$ to $ \boldsymbol {\xi } _{\textrm {s}}$ . The combined selection and focus field is shown in the last row. In an ideal MPI system, the coefficients with light blue background would be identical.

4.1 Magnetic fields in MPI

4.1.1 Static fields

Using translation of the coefficients enables comparison of the different magnetic fields applied in MPI. In Figure 6, the static magnetic fields of the central and the lower left patch of Figure 1 are shown. In the first row, the coefficients of the selection field with a gradient strength of $\text{2 Tm}^{-1}$ at its FFP $ \boldsymbol {\xi } _{\textrm {c}}$ and at the centre point of the lower left patch $ \boldsymbol {\xi } _{\textrm {s}}$ are shown. While the coefficients at $ \boldsymbol {\xi } _{\textrm {c}}$ do not feature any significant imperfections, in the shifted coefficients some imperfections for $l \geq 1$ occur. Since $ \boldsymbol {\xi } _{\textrm {s}}$ is not the FFP of the selection field, $ \boldsymbol {\gamma\!}_{0,0}^{\textrm {SF}}( \boldsymbol {\xi } _{\textrm {s}},R_{\textrm {s}})$ are non-zero and contain information about the offset field. Using an additional focus field with the same offset field but opposite sign, the offset field can be cancelled and the FFP is shifted into $ \boldsymbol {\xi } _{\textrm {s}}$ . This focus field is shown in the second row with an offset field of −24 mT in $x$ - and 24 mT in $z$ -direction. At both positions some imperfections occur, but they are slightly higher at the off-centre position $ \boldsymbol {\xi } _{\textrm {s}}$ . The combined selection and focus field is visualised in the last row. The axes of the field plot on the right are shifted due to the translation of the coefficients to the FFP $ \boldsymbol {\xi } _{\textrm {s}}$ . The FFP therefore has the coordinates $(0,0)$ , as it is in the top field plot. Due to the shift into the FFP, the coefficients of the initial selection field and the combined field can be directly compared. It can be observed that the combined field has much more imperfections than the initial selection field, starting already from $l=1$ .

The coefficients do not only enable comparison, and they also allow for calculation of the real gradient strength and focus field shifts differing from the input parameters given in Section 3.3. In case of the set-up of Figure 6, the real gradient strength of the selection field in its FFP is −1.02, −1.01 and $\text{2.03 Tm}^{-1}$ in $x$ -, $y$ - and $z$ -direction. Meanwhile, the real focus field shifts are −22.96 and 23.28 mT in $x$ - and $z$ -direction, respectively.

4.1.2 Dynamic fields

The 12 mT dynamic fields of our MPI scanner are shown in Figure 7. As for the static fields, the coefficients at the FFP of the selection field $ \boldsymbol {\xi } _{\textrm {c}}$ and at the shifted FFP $ \boldsymbol {\xi } _{\textrm {s}}$ are shown in the left columns, while the $x$ -, $y$ - and $z$ -drive fields are shown in the three rows. Overall, the coefficients decrease as $l$ increases, which justifies truncating the expansion at $L=4$ . It can be observed that even in the centre imperfections especially for $l \gt 1$ occur, which are more severe for the $y$ - and $z$ -drive field. The coefficients for $l=0$ show that the drive-field amplitudes 12.35, 12.29 and 12.15 mT for the $x$ -, $y$ - and $z$ -drive field deviate from the 12 mT input parameter. This is because the drive-field coils are located closer to the scanner bore, than the selection- and focus-field coils. Comparing the coefficients at the central FFP and at the shifted FFP, the imperfections of the drive fields increase. Here, imperfections also arise for $l=0$ , which are especially visible for the $z$ -drive field. In the shifted FFP, the constant part of the $z$ -drive field is not perfectly aligned in $z$ -direction but also points slightly in the $x$ -direction. In combination with the imperfections of the selection and focus field, this leads to the distorted trajectory of the lower left patch in Figure 1. The imperfections also manifest in the field plot on the right, where for each drive field a representative plane is shown.

Figure 7. Solid harmonic analysis of the dynamic fields in MPI, that is, the drive fields. Three drive fields in $x$ -, $y$ - and $z$ -direction with 12 mT amplitude are shown in each row. In the left columns, the coefficients up to $L=\textit{3}$ at the FFPs of the selection fields of Figure 6 are visualised, while on the right, the fields in the $xz$ - respectively $xy$ -plane are shown.

4.2 Error analysis

A directional comparison of the field values provided by the truncated expansion to the measured ones at different gradient strengths shows an error (standard deviation) in the range of 4 to 100 $\mu$ T. As shown in Figure 8, this error increases with the gradient strength, is approximately the same in $y$ - and $z$ -direction and is approximately a factor of 2 smaller in $x$ -direction. If we normalise the error as defined in equation (6), one observes errors in the range of $4\cdot 10^{-4}$ to $17\cdot 10^{-4}$ . The largest normalised errors can be observed for the smallest gradient strength, but no clear trend is evident for the remaining gradient strengths. Concerning the spatial dependence, the same observations apply which we just made.

If we propagate the uncertainty of the calibration measurement, we have to expect errors in the range of 5 to 18 $\mu$ T, which increase with the gradient strength, as shown in Figure 8. These are similar in the $x$ - and $y$ -directions and stronger in the $z$ -direction in the range of 20% to 70%, increasing linearly with gradient strength. In direct comparison, the observed error is up to a factor of 6 larger then the propagated one. Hence, the observed error can only partially attributed to uncertainties in the measurements with our Hall-effect sensor.

Figure 8. Standard deviation of the measurement errors (solid lines) and propagated errors (dashed lines).

5. Discussion

In this paper, we have given a review of the real solid harmonic expansions as a general solution of Laplace’s equation and efficient quadrature methods for the calculation of the expansion coefficients via spherical t-designs. Furthermore, we proposed a method to change the reference point of the expansion using spatial shifts. This allows for a unique expansion of the magnetic field that is independent of the measurement set-up. Furthermore, we have illustrated a methodology for the analysis of field imperfections via the polynomial structure of the expansion around the reference point. This is based on the premise that the local properties of the magnetic fields are directly encoded in the coefficients of the solid harmonic expansions. Our methods were evaluated on the signal generating and encoding fields in MPI, where the coefficients provide a compact representation of the fields using the characteristic FFP of the static selection field as unique expansion centre. This uniqueness allows for comparison of the behaviour of magnetic fields of different patches or MPI scanner.

One of the main advantages of using these truncated expansions for the approximation of magnetic fields is the extremely fast acquisition time. Contrary to classical methods where the magnetic field is densely measured on the entire FOV, fewer measurements on the surface of a sphere suffice to obtain the truncated expansion into real solid harmonics, regardless of the polynomial degree of the underlying field. In our set-up, measuring at the positions of a spherical $8$ -design takes about 2 min, which is sufficient to approximate the static and dynamic fields in MPI. In comparison, a Gauss–Legendre quadrature scheme, which has been commonly used in MPI scenarios, has 25% more nodes and would take about 2.5 min.

Our error analysis has shown that the deviations between model (truncated expansion) and measurement on the sphere are in the per mill range. However, only a small part of the observed deviations, 20% in the worst case scenario, can be attributed to the measurement inaccuracy of the Hall-effect sensor used. The error source with the greatest influence must therefore have a different origin. For example, a model error caused by the truncation of the expansion is possible. This hypothesis could be tested for example by choosing an expansion with larger $L$ . For this, of course, a new spherical t-design with $t=2L$ would have to be chosen. However, for the application in the MPI context targeted in this work, the approximation accuracy of the real solid harmonic expansion up to degree $L=4$ achieved in this work is sufficient. Meanwhile, an analysis of the uncertainty of the FFP position is still an open question. For example, Monte Carlo sampling techniques can be used to propagate the error from the field to the FFP position [32, 3].

The coefficients enable straightforward analysis of the different field setting of an MPI system. The local field characteristics are directly encoded in the coefficients, enabling direct comparison of the local behaviour of the magnetic fields without the need for additional field evaluations. In the shifted spatial encoding and excitation fields of the presented MPI scanner, we observe slight imperfections as shown in Figures 6 and 7, respectively. These imperfections are a major cause for imaging artefacts and their precise knowledge is key in their reduction. In some cases, the knowledge about the magnetic fields has already been incorporated into other projects or is future work. First of all, the knowledge can be exploited for MPI measurement planning. Shifting a patch to the correct position is crucial in many scenarios like multi-resolution data acquisition [18] or magnetic actuation [37]. Especially in multi-patch MPI, the distorted shape and position of the shifted patches can lead to uncovered areas inside the FOV as shown in Figure 1 or to imaging artefacts when the imperfections are not included in the patch-wise imaging operator [46]. By incorporating the magnetic fields presented here, the latter can be avoided by dedicated measurements of the operator of each patch for non-negligible field deviations [9], a field-dependent post-processing of the operator [8] or modelling of the imaging operator with integrated field imperfections [1, 2]. In addition, solid harmonic expansions are already being used to analyse the magnetic fields of newly built field generators for MPI [49, 17]. Furthermore, they can be directly incorporated in the reconstruction process [11].

However, using the FFP of the selection field as the unique expansion point does not work for MPI scanners with multiple FFPs or an FFL. For these set-ups, other unique points must be used, such as the robot’s zero position. Although this does not provide coefficients comparable to other scanners, it does allow for comparison and characterisation of different fields of the same scanner. Furthermore, any rotation of the set-ups is not captured by our proposed method. Steinborn and Ruedenberg [45] introduce not only the translation but also the rotation of solid spherical harmonics. A comparable methodology could be employed with regard to the rotation. This would allow for any rotation of the sensors or robot to be addressed by mapping the expansion coefficients to those of a rotated expansion.

Medical imaging set-ups often feature cylindrical gantries, so an expansion of the magnetic fields with cylindrical harmonics would be a natural choice. Nevertheless, an approximation with cylindrical harmonic expansions requires considerably more coefficients, which complicates the analysis of the magnetic field imperfections. Furthermore, we observed that even more coefficients are required to obtain a similar accuracy as with a spherical harmonic expansion since the basis functions of the cylindrical harmonic expansions are not as suitable for the presented magnetic fields. To the best of our knowledge, such a small set of quadrature nodes as the spherical t-design does not exist for quadrature on a cylindrical surface. However, ellipsoidal harmonic expansions are a viable option, as they more closely resemble the cylindrical shape of the gantries than the spherical harmonics. By transferring the spherical t-designs to an ellipsoidal surface, it is possible to calculate the requisite ellipsoidal coefficients while maintaining the original measurement time. A proof of concept for this methodology can be found in the work of Scheffler et al. [43].

The compact representation of the magnetic fields in MPI offer multiple further investigations. Using the presented tools we can deal with the field’s imperfections in various applications. First of all, they are an important parameter for model-based reconstructions. Incorporating the real field parameter into the modelled imaging operator lead to reconstruction results closer to those obtained with a measured operator. Thus, with the provided tools, we can work with the imperfections of the magnetic fields instead of avoiding them at all costs. They can be accounted for in the imaging sequences or for magnetic actuation. Furthermore, the spherical t-design offers sufficient small set of measurement points such that multiple Hall-effect sensors can be used simultaneously to measure a magnetic field in one shot. This can be used for direct feedback for magnetic field calibrations.