MMathematical models for magnetic particle imaging
Tobias Kluth ∗ March 21, 2018
Abstract
Magnetic particle imaging (MPI) is a relatively new imaging modality. The nonlinear magnetizationbehavior of nanoparticles in an applied magnetic field is employed to reconstruct an image of the concen-tration of nanoparticles. Finding a sufficiently accurate model for the particle behavior is still an openproblem. For this reason the reconstruction is still computed using a measured forward operator whichis obtained in a time-consuming calibration process. The state of the art model used for the imagingmethodology and first model-based reconstructions relies on strong model simplifications which turnedout to cause too large modeling errors. Neglecting particle-particle interactions, the forward operatorcan be expressed by a Fredholm integral operator of the first kind describing the inverse problem. In thisarticle we give an overview of relevant mathematical models which have not been investigated theoreti-cally in the context of inverse problems yet. We consider deterministic models which are based on thephysical behavior including relaxation mechanisms affecting the particle magnetization. The behavior ofthe models is illustrated with numerical simulations for monodisperse as well as polydisperse tracer. Wefurther motivate linear and nonlinear problems beyond the solely concentration reconstruction related toapplications. This model survey complements a recent topical review on MPI [30] and builds the basisfor upcoming theoretical as well as empirical investigations.
Magnetic particle imaging (MPI) is a relatively new imaging modality [13] which relies on the behavior ofsuperparamagnetic iron oxide nanoparticles. The main goal is to reconstruct the spatially dependent concen-tration of the particles. Measurements are obtained from multiple receive coils where a potential is inducedby the particle’s nonlinear response to the applied dynamic magnetic field. These potential measurementsare used for image reconstruction. A high temporal resolution and a potentially high spatial resolutionmake MPI suitable for several in-vivo applications without the need for harmful radiation. Applications alsobenefit from the fast data acquisition of MPI.The number of potential medical applications is still increasing. First proposed medical applications arevascular imaging and medical instrument tracking. The potential of imaging blood flow was shown first in invivo experiments using a mouse [60]. The usability of a circulating tracer for long term monitoring has beeninvestigated recently [21]. Another medical application which benefits from the high temporal resolution istracking medical instruments [18]. Recently it was shown that MPI is also suitable for tracking and guidinginstruments for angioplasty [49]. Further promising applications of MPI can be found in cancer detection[61] and cancer treatment by hyperthermia [40].The relationship between particle concentration and measured potential is modeled by a Fredholm integralequation of the first kind which is motivated by the suppression of particle interactions due to the nonmagneticcoating. Determining the concentration is thus a linear inverse problem. The integral kernel is defined bythe properties of the receive coils, the applied magnetic fields and the dynamic behavior of the tracermaterials. MPI imaging methodologies are characterized by the applied magnetic fields which are generatedby moving a field free point (FFP) [13] or a field free line (FFL) [58] along a given trajectory rapidly. Themost prominent FFP trajectories are Lissajous and Cartesian trajectories but also other trajectories wereinvestigated in simulation studies [25]. Cartesian sequences [7, 34, 62] allow modeling the MPI signal by a ∗ Center for Industrial Mathematics, University of Bremen, Bibliothekstr. 5, 28357 Bremen, Germany( [email protected] ) a r X i v : . [ phy s i c s . a pp - ph ] M a r patial convolution. In this case deconvolution methods are commonly used to obtain image reconstructions,also known as x-space reconstruction [14, 15]. In contrast, using a Lissajous trajectory, measurements can beobtained faster but the corresponding system matrix obeys a complex structure. In the absence of suitablemodels, the system matrix is usually measured [33, 17] for different scanner setups and tracer materials in atime-consuming measurement process where a “delta” probe is moved through the field of view. Besides thetime-consuming measurement process the methodology suffers from high memory requirements particularlyfor three-dimensional imaging.The problem of modeling MPI, respectively finding the correct integral kernel, also known as systemfunction, is an unsolved problem. Existing model-based reconstructions incorporate particle behavior basedon the theory of paramagnetism [33, 27, 14, 39, 24]. Methods based on ideal magnetic fields [43, 14] andon realistic magnetic fields [27] are promising but they are not yet able to reach the quality of measuredsystem functions. One possible reason may be a nonlinear dependence on the concentration which have beenreported for large concentrations [36]. But as the affected concentrations are larger than the concentrationscommonly used in MPI experiments, these effects are not included in existing models used for imaging sofar. Another possible reason which arouses increasing interest in MPI is the particle relaxation which islikely to emerge for rapidly changing applied magnetic fields in MPI. Simplified models given by ordinarydifferential equations are combined with the superparamagnetic particle behavior [7] to deal with this issue.The authors report a fitted relaxation time of 2 . µ s for Resovist in an applied field of 22 . This sections aims at technical definitions required for the formulation of MPI models. For a more detailedintroduction to the principles of MPI the interested reader is referred to [28]. A recent review of MPI tracerdevelopment can be found in [2] and for a detailed chronological overview on the technical developments inMPI during the last decade we refer to [30].
The inherent nature of the problem is three-dimensional which is the reason why the relevant vector valuedfunctions remain three dimensional even if the domain of the spatial variable is a subset of a d -dimensionalaffine subspace E d ⊂ R . Let Ω ⊂ E d , d = 1 , ,
3, be a bounded domain in E d and let S = { x ∈ R |(cid:107) x (cid:107) = 1 } be the unit sphere. Furthermore, let T > I = [0 , T ] the timeinterval during which the measurement process takes place. The temporal partial derivative of a function g : I → R k , k ∈ N , is denoted by ˙ g , i.e., ˙ g = ∂∂t g . In the following we define the problem and subsequently describe how the lower dimensional case is con-structed. The signal v k : I → R , k = 1 , . . . , L , which is obtained from the L ∈ N receive coils, is givenby v k ( t ) = − (cid:90) I (cid:90) Ω c ( x )˜ a k ( t − τ ) s k ( x, t ) dx dτ − (cid:90) I (cid:90) R ˜ a k ( t − τ ) p k ( x ) T ˙ B app ( x, t ) dx dτ (cid:124) (cid:123)(cid:122) (cid:125) = v E ,k ( t ) (1)where c : Ω → R + ∪ { } is the concentration of the magnetic nanoparticles and s k : Ω × I → R , k = 1 , . . . , L ,denote the system functions characterizing the behavior of the nanoparticles. a k : I → R , k = 1 , . . . , L , arethe periodic kernel functions of the analog filter in the signal acquisition chain and ˜ a k denotes its periodiccontinuation. p k : R → R , k = 1 , . . . , L , denotes the vector field which characterizes the sensitivityprofile of the receive coils. In the remainder of the article it is assumed that the magnetic flux density B app : R × I → R of the applied magnetic field and the kernel functions a k , k = 1 , . . . , L , are chosen in away such that v E ,k = 0 holds for all excitation signals v E ,k : I → R , k = 1 , . . . , L , as defined in (1). Remark 2.1.
Common choices for the analog filters a k , k = 1 , . . . , L , are band stop filters adapted to thefrequencies of sinusodial excitations used in the subsequently described drive field. The assumption regardingthe excitation signals v E ,k , k = 1 , . . . , L , is commonly made when the structure of the system functions isstudied but it is not fulfilled in MPI applications in general [55]. Efforts are made to remove the excitationsignal or reduce its influence on the concentration reconstruction [32, 57, 55, 24]. The applied magnetic fields used in MPI can be characterized ideally by a spatially dependent magneticfield g : R → R and a time-dependent homogeneous magnetic field h : I → R . The applied magnetic fieldis then given by their superposition, i.e., B app ( x, t ) = g ( x ) + h ( t ). The ideal case relies on the assumptionthat the sensitivity profile of the field generating coils is homogeneous, respectively linear, in the field ofview. The field g , also known as selection field , guarantees that a field-free-region is generated. Ideally, g is assumed to be linear such that it can be represented by its transformation matrix G ∈ R × . Here twocases for the methodology are distinguished, namely whether a FFP is generated (rank( G ) = 3) or a FFL is3sed (rank( G ) = 2). The field h , also known as drive field , then moves the field-free-region along a certaintrajectory. Remark 2.2.
In the literature it was also proposed for the FFL approach to rotate the selection field over timesuch that the FFL is rotated [58, 29]. In [29] the authors also show a relation between the Radon transformand the FFL approach combined with the particle model based on the Langevin function (subsequently termedequilibrium model). The selection field is then given by g : R × I → R with g ( x, t ) = P ( t ) T GP ( t ) x where P : I → R × is a rotation matrix for all t ∈ I . Remark 2.3.
In common MPI applications the drive field is defined by weighted sine or cosine functionswith different frequencies in each component of h . But also other kinds of excitation signals were investigatedin simulations [25]. We now can formulate the most general version of the MPI problem. It is then to find the concentration c which fulfills v k ( t ) = − (cid:90) I (cid:90) Ω c ( x )˜ a k ( t − t (cid:48) ) s k ( x, t (cid:48) ) dx dt (cid:48) s k = µ p Tk ddt ¯ m (2)for k = 1 , . . . , L and where ¯ m : Ω × I → R is the mean magnetic moment vector of the nanoparticles.The 1st equation in (2) models the analog filter process by a convolution with respect to t (cid:48) . The spatialintegration describes the induction of a potential in the receive coil and is obtained from Faraday’s law ofinduction combined with the law of reciprocity to obtain the sensitivity profile of the coils [28]. The 2ndequation comprises the sensitivity of the receive coil and the particle behavior in the applied magnetic field.The d -dimensional case for d < δ -distributionwith respect to the orthogonal complement of the affine subspace E d . Then the d -dimensional problem isconstructed by assuming c ( x ) = ˜ c ( x ) δ ( x ) where x = x + x with x ∈ E d , x ∈ E ⊥ d , and ˜ c : Ω ⊂ E d → R + ∪ { } . The parametrization of Ω ⊂ E d then allows reformulating the spatial integral in (2) as an integralover an integration domain Ω d ⊂ R d . Given the affine linear parametrization Γ : Ω d → Ω we can considerthe problem with respect to c d : Ω d → R + ∪ { } , c d ( x ) = ˜ c (Γ( x )). All other spatially dependent functionsare then treated analogously. Remark 2.4.
The linear dependence on the concentration c is based on the assumption that the particle-particle interactions can be neglected. There is increasing evidence that demagnetization effects relying onthese interactions can significantly influence the particle signal [36]. The nonlinear dependence is not con-sidered in this article and remains to be explored in future work. Now we are able to formulate various MPI models for different particle behavior and tracer compositions.For the subsequent considerations we assume single-domain particles which is reasonable for sufficientlysmall diameters of the ferromagnetic core [6]. In this case each particle has a uniform magnetization for anyapplied magnetic field. Here the magnetic moment is then considered for the whole single-domain particleinstead of single atoms. In general the single-domain particles are not isotropic due to the particle’s shape,internal stress, and the internal structure [6]. In the following it is assumed that the single-domain particleshave an uniaxial anisotropy.
For monodisperse tracers it is assumed that the tracer material consists of one single kind of nanoparticlesonly, and all nanoparticles have the same characteristic behavior. Two possible mechanism of relaxationare taken into account for MPI following [48, 45]. Based on the previous assumptions on the particlesit is assumed that each particle is equipped with its magnetic moment vector which follows the externalmagnetic field. The particle can change the orientation of its magnetic moment vector by a rotation of thewhole particle, known as Brownian rotation [5], and by an internal rotation, known as N´eel relaxation [41, 3].4irst, we formulate the simplified model based on the Langevin function, termed equilibrium model . Thenthe general MPI problem is specified with respect to the relaxation mechanisms resulting in MPI modelswhich have not been studied with respect to the imaging problem yet.
One of the most extensively studied models in MPI is based on the Langevin function. It is the only modelwhich has been studied with respect to the imaging problem so far. The model is motivated by the assump-tions that the applied magnetic fields are static and the particles are in equilibrium. These assumptions alsomotivate the term equilibrium model which is used in the following. Using these assumptions, it is assumedthat the mean magnetic moment vector of the nanoparticles immediately follows the magnetic field, i.e.,¯ m ( x, t ) = L β ( (cid:107) B app ( x, t ) (cid:107) ) B app ( x, t ) (cid:107) B app ( x, t ) (cid:107) (3)where L β : R → R is given in terms of the Langevin function by L β ( z ) = m (cid:18) coth( βz ) − βz (cid:19) (4)for m , β >
0. The Langevin function determines the length of the mean magnetic moment vector inequilibrium and can be derived from the Brownian rotation model stated below by assuming a static magneticfield. A more detailed description of this relationship can be found in B. The final problem based on theLangevin function then is to obtain the concentration c from the following system of equations: v k ( t ) = − (cid:90) I (cid:90) Ω c ( x )˜ a k ( t − t (cid:48) ) s k ( x, t (cid:48) ) dx dt (cid:48) s k = µ p Tk (cid:18)(cid:20)(cid:18) L (cid:48) β ( (cid:107) B app (cid:107) ) (cid:107) B app (cid:107) − L β ( (cid:107) B app (cid:107) ) (cid:107) B app (cid:107) (cid:19) B app B T app + L β ( (cid:107) B app (cid:107) ) (cid:107) B app (cid:107) I (cid:21) ˙ B app (cid:19) (5)for k = 1 , . . . , L and the identity matrix I ∈ R × . Remark 3.1. m and β are determined by the saturation magnetization M C of the core material, thevolume of the particle’s core V C , the temperature T B , and the Boltzmann constant κ B , i.e., m = M C V C and β = m / ( κ B T B ) . Assuming spherical particles, the influence of the particle diameter D is given by V C = 1 / πD . Note that β also depends on D as it depends on m . For example, particles consistingof magnetite with a typical diameter of nm (20 nm ) at room temperature K are characterized by β ≈ . /µ × − (0 . /µ × − ) . Remark 3.2.
Under certain assumptions on B app , this problem can be formulated in terms of a spatialconvolution evaluated along the trajectory of the field free point with a temporally changing convolutionkernel [39, 24]. A first theoretical investigation related to this model can be found in [39]. Here we assume that each particle has a magnetic moment vector which changes its direction due to therotation of the whole particle [5, 45]. Including the Brownian rotation and thermal noise results in a Langevinequation from which a deterministic system of differential equations is derived to model the mean magneticmoment vector. The dynamics of a particle’s magnetic moment vector ˜ m : Ω × I → m S with (cid:107) ˜ m (cid:107) = m isgiven by ∂∂t ˜ m = νm (cid:16) ˜ m × B app + ˜ D Γ (cid:17) × ˜ m (6)5ith physical parameters ν, ˜ D > (cid:104) Γ i ( t ) (cid:105) = 0, (cid:104) Γ i ( t )Γ j ( t ) (cid:105) = δ ij δ ( t − t ), for all t, t , t > i, j = 1 , , (cid:104)·(cid:105) denotes theexpectation value of a random variable. In case of zero noise the magnetic moment vector is directly dampedinto the direction of the magnetic field. The Fokker-Planck equation for the probability density function f : S × Ω × I → R + ∪ { } is used to derive the mean magnetic moment vector. The problem then becomes v k ( t ) = − (cid:90) I (cid:90) Ω c ( x )˜ a k ( t − t (cid:48) ) s k ( x, t (cid:48) ) dx dt (cid:48) k = 1 , . . . , Ls k = µ p Tk ∂ ¯ m∂t k = 1 , . . . , L ¯ m ( x, t ) = m (cid:90) S mf ( m, x, t ) dm in Ω × I∂f∂t = − div m (cid:18) ν ( B app − m ( m · B app )) f − τ ∇ m f (cid:19) in S × Ω × I (cid:90) S f ( m, x, t ) dm = 1 in Ω × If ( · , x,
0) = f in Ω (7)for m , ν, τ > f : S → R + ∪ { } with (cid:82) S f dm = 1 is the initial distribution function. The1st and 2nd equations of (7) describe the general problem as stated in (2). The 3rd equation defines themean of the magnetic moment vector in terms of the probability density function f . The 4th equation is thedifferential equation for the probability density function which takes into account the Brownian rotation.The derivation of the Fokker-Planck equation starting from the Langevin equation can be found in A. The5th equation guarantees that f is a probability density function and the sixth equation is the initial condition. Remark 3.3. m , ν , and τ are determined by the saturation magnetization M C of the core material, thevolume of the core V C , the hydrodynamic volume V H of the particles, the dynamic viscosity η , the temperature T B , and the Boltzmann constant κ B , i.e., m = M C V C , ν = m V H η , and τ = 3 V H η/ ( κ B T B ) . τ is known as therelaxation time. Assuming spherical particles, the particle diameter D influences V C as well as V H . Thus m , ν , and τ depend on the particle diameter. An overview of common particle parameters which are obtainedfrom [8] can be found in Table 1. Remark 3.4.
Possible further applications in MPI are motivated by multi-color MPI [42] which also dis-tinguishes different kinds of tracer by their relaxation behavior [56]. In this case the characteristic behaviorencoded in multiple system matrices allow the distinction. In the Brownian rotation model it corresponds tothe simultaneous reconstruction of c and a spatially dependent viscosity η . As a result a nonlinear inverseproblem must be solved. Assume that the particles are spatially blocked such that Brownian rotation is suppressed. Each particlehas a magnetic moment vector which changes its direction due to the change of internal electric states[41, 3, 45]. Including the applied magnetic field, the particle’s anisotropy, and thermal noise results in aLangevin equation based on the Landau-Lifshitz-Gilbert equation for the particle’s magnetic moment vector˜ m : Ω × I → m S with (cid:107) ˜ m (cid:107) = m given by ∂∂t ˜ m = ˜ γ (cid:18) ( B eff + ˜ D Γ) × ˜ m + αm ( ˜ m × ( B eff + ˜ D Γ)) × ˜ m (cid:19) (8)with physical parameters ˜ γ, α, ˜ D > (cid:104) Γ i ( t ) (cid:105) = 0, (cid:104) Γ i ( t )Γ j ( t ) (cid:105) = δ ij δ ( t − t ), for all t, t , t > i, j = 1 , ,
3. In contrast to Brownian rotation the magnetic momentvector moves on a precessional trajectory while it is damped in direction of the magnetic field. From thisequation a deterministic system of differential equations is derived to model the mean magnetic moment6ector ¯ m . In contrast to the previous models, an effective magnetic field B eff : S × Ω × I → R is consid-ered which consists not only of the applied magnetic field. Here the uniaxial anisotropy of the particles ismodeled by a field B anis : S → R with B anis ( m ) = 2 K anis V C m ( m · n ) n for a given anisotropy direction n ∈ S and anisotropy constant K anis ∈ R . This field is obtained from the Stoner-Wohlfarth model for uniaxialanisotropic particles [52, 54]. Here K anis > K anis < B eff = B app + B anis . The Fokker-Planck equation for the probabilitydensity function f : S × Ω × I → R + ∪ { } is used to derive the mean magnetic moment vector. The problemthen becomes v k ( t ) = − (cid:90) I (cid:90) Ω c ( x )˜ a k ( t − t (cid:48) ) s k ( x, t (cid:48) ) dx dt (cid:48) k = 1 , . . . , Ls k = µ p Tk ∂ ¯ m∂t k = 1 , . . . , L ¯ m ( x, t ) = m (cid:90) S mf ( m, x, t ) dm in Ω × I∂f∂t = − div m (cid:18) ˜ γ ( B eff × m + α ( B eff − ( m · B eff ) m )) f − τ ∇ m f (cid:19) in S × Ω × I (cid:90) S f ( m, x, t ) dm = 1 in Ω × If ( · , x,
0) = f in Ω (9)for m , ˜ γ, α, τ > f : S → R + ∪ { } with (cid:82) S f dm = 1 is the initial distribution function. The1st and 2nd equation of (9) describe the general problem as stated in (2). The 3rd equation defines the meanof the magnetic moment vector with respect to the probability density function f . The 4th equation is thedifferential equation for the probability density function which takes into account the N´eel relaxation. Theeffective magnetic field is B eff = B app + 2 K anis V C m ( m · n ) n , K anis , V C > n ∈ S . Depending on the choiceof n ∈ S , the anisotropy field may counteract the applied magnetic field. The larger the applied magneticfield strength the weaker is the influence of the anisotropy field. The derivation of the Fokker-Planck equationstarting from the corresponding Langevin equation can be found in C. The 5th equation guarantees that f is a probability density function and the sixth equation is the initial condition. Remark 3.5. m , ˜ γ , and τ are determined by the saturation magnetization M C of the core material, thevolume of the core V C , the gyromagnetic ratio γ , the damping parameter α , the temperature T B , and theBoltzmann constant κ B , i.e., m = M C V C , ˜ γ = γ α , and relaxation time τ = m α ˜ γκ B T B . Assuming sphericalparticles, the particle diameter D influences V C . Thus m and τ depend on the particle diameter. Note that B anis does not depend on the diameter as the volume of the core V C cancels out. An overview of commonparticle parameters which are obtained from [8] can be found in Table 1. Remark 3.6.
The temperature influences the behavior of the nanoparticles. In the context of multi-colorMPI it was also suggested to determine the temperature while reconstructing the concentration [51]. Theauthors motivate the simultaneous reconstruction by real time monitoring in hyperthermia applications. Theproblem of finding a spatially dependent temperature T B has a nonlinear nature requiring a deeper analysis ofthe problem. Simultaneous usage of different kinds of particles differing in their physical properties result inanother possible nonlinear problem similar to the case of multi-color MPI with viscosity mapping. Structuraldifferences might also be found in the relaxation times τ indicating different environmental structures, e.g.,whether particles are blocked or non-blocked. In MPI, polydisperse tracer has been proposed in context of the equilibrium model for reconstruction so far.The investigations of tracer materials commonly used in MPI motivate the introduction of a distribution ofthe core diameter [33] which is mainly motivated by related physical investigations [23, 10, 9, 22]. The tracermaterial is then modeled by a distribution of particles differing in their diameter
D >
0. Assuming that theparticle’s diameter distribution is given by a density function ρ : R + → R + ∪ { } with (cid:107) ρ (cid:107) L ( R + ) = 1, the7xtended problems are obtained. Below we state these models to highlight the parameter dependence onthe particle’s core diameter. For the same reasons as for the monodisperse problems with Brownian rotationand N´eel relaxation their polydisperse extensions have not been considered for imaging yet. Remark 3.7.
The distribution function of the particle diameters for polydisperse tracers can be approximatedby a log-normal distribution [23, 10, 9, 22]. For example, for Resovist a mean µ = ln(13 × − ) anda standard deviation of σ = 0 . of a log-normal distribution were reported for the diameter [10], i.e. ρ ( D ) = Dσ √ π e − (ln( D ) − µ ) / (2 σ ) . The equilibrium model stated in (5) can be extended to polydisperse tracers by adapting the functiondefining the length of the mean magnetic moment vector in (4). The extended problem is then given by v k ( t ) = − (cid:90) I (cid:90) Ω c ( x )˜ a k ( t − t (cid:48) ) s k ( x, t (cid:48) ) dx dt (cid:48) s k = µ p Tk (cid:18)(cid:20)(cid:18) L (cid:48) β,ρ ( (cid:107) B app (cid:107) ) (cid:107) B app (cid:107) − L β,ρ ( (cid:107) B app (cid:107) ) (cid:107) B app (cid:107) (cid:19) B app B T app + L β,ρ ( (cid:107) B app (cid:107) ) (cid:107) B app (cid:107) I (cid:21) ˙ B app (cid:19) (10)where L β,ρ : R → R is given in terms of the Langevin function by L β,ρ ( z ) = (cid:90) R + ρ ( D ) m ( D ) (cid:18) coth( β ( D ) z ) − β ( D ) z (cid:19) dD (11)for m , β : R + → R + describing the influence of the particle diameter on the volume of the core, respectivelythe magnetic moment. For further details on the physical paramters and their relationship to the particlediameter, we refer to Remark 3.1.The Brownian rotation model stated in (7) can be extended for polydisperse tracers by extending thedomain of the magnetic moment vector’s probability density function f : S × Ω × I × R + → R + ∪ { } . The Brownian rotation model then becomes v k ( t ) = − (cid:90) I (cid:90) Ω c ( x )˜ a k ( t − t (cid:48) ) s k ( x, t (cid:48) ) dx dt (cid:48) k = 1 , . . . , Ls k = µ p Tk ∂ ¯ m∂t k = 1 , . . . , L ¯ m ( x, t ) = (cid:90) R + ρ ( D ) m ( D ) (cid:90) S mf ( m, x, t, D ) dm dD in Ω × I∂f∂t = − div m (cid:18) ν ( D ) ( B app − m ( m · B app )) f − τ ( D ) ∇ m f (cid:19) in S × Ω × I × R + (cid:90) S f ( m, x, t, D ) dm = 1 in Ω × I × R + f ( · , x, , D ) = f in Ω × R + (12)for m , ν, τ : R + → R + and where f : S → R + ∪ { } with (cid:82) S f dm = 1. In contrast to the monodispersemodel, the 3rd equation of (12) comprises the mean of the mean magnetic moment vector and the mean overthe particle diameter. The differential equation for the probability density function in the 4th, 5th, and 6thequation are extended to take the dependence on the particle diameter D into account. For further detailson the physical parameters and their relation to the particle diameter we refer to Remark 3.3.The N´eel relaxation model is extended analogously by considering f : S × Ω × I × R + → R + ∪ { } . It is8hus given by v k ( t ) = − (cid:90) I (cid:90) Ω c ( x )˜ a k ( t − t (cid:48) ) s k ( x, t (cid:48) ) dx dt (cid:48) k = 1 , . . . , Ls k = µ p Tk ∂ ¯ m∂t k = 1 , . . . , L ¯ m ( x, t ) = (cid:90) R + ρ ( D ) m ( D ) (cid:90) S mf ( m, x, t, D ) dm dD in Ω × I∂f∂t = − div m (cid:18) ˜ γ ( B eff × m + α ( B eff − ( m · B eff ) m )) f − τ ( D ) ∇ m f (cid:19) in S × Ω × I × R + (cid:90) S f ( m, x, t, D ) dm = 1 in Ω × I × R + f ( · , x, , D ) = f in Ω × R + (13)for m , τ : R + → R + , ˜ γ, α > f : S → R + ∪ { } with (cid:82) S f dm = 1, cf. Remark 3.5 forfurther details on the physical parameters and their dependence on the diameter D . The 3rd equation of(13) comprises the mean of the mean magnetic moment vector and the mean over the particle diameter andthe 4th, 5th, and 6th equation are extended to take the dependence on the particle diameter into account. Remark 3.8.
The polydisperse models are parametric equations which are similar to the monodispersemodels based on the stochastic differential equations.
To illustrate the behavior of different models we computed the temporal derivative of the mean magneticmoment vector ¯ m in one special case where the probability density function is circular symmetric with respectto the e -axis. On the e -axis the third component of the mean magnetic moment vector is nonzero only.The numerical solution is based on the formulation of a system of ordinary differential equations obtainedby approximating the probability density functions for Brownian and N´eel relaxation by a finite number ofLegendre polynomials. Further details about the numerical solution which follows [8] can be found in D.In the subsequent simulations the first 50 Legendre polynomials were used. The initial distribution f isassumed to be uniform.The simulation setup is as follows. We assume an excitation in the direction of e only and considerthe case that Ω ⊂ E d = { qe | q ∈ R } . The drive field is given by h ( t ) = A cos(2 πf t ) e where A > f > G ∈ R × .For the simulations we use physical parameters typical for MPI applications. An overview of the param-eters can be found in Table 1. The remaining parameters are computed according to Remarks 3.1, 3.3, and3.5.The simulations for the monodisperse models are illustrated in Figure 1 for a diameter of 20 nm. Thegraphs of the equilibrium model are point symmetric to the point (1 . ,
0) which is a direct result of thecosine excitation. The extremal points in the equilibrium model are close to the field free point in time.The Brownian rotation model shows damped and skewed behavior in time direction when it is comparedto the equilibrium model. This is due to the viscous rotation of the particles. Neglecting the shift intime, the N´eel relaxation shows a similar qualitative behavior in terms of temporal symmetry at the origin( z = 0 mm) compared to the equilibrium model. The shift in time may be related to the particle anisotropywith preferred e -direction. A certain strength of the magnetic field is required before the magnetizationchanges. More structural differences can be observed at the remaining z -positions where one extremal pointis more damped than the other one. The more it is damped the less is the applied magnetic field able tocounteract the anisotropy in the preferred direction of anisotropy.The simulations for the polydisperse models are illustrated in Figure 2. The polydisperse equilibriummodel is smaller in amplitude but the qualitative behavior remains similar to the monodisperse version with9 arameter Value Magnetic permeability µ π × − H/mBoltzmann constant κ B . × − J/K
Scanner
Excitation frequency f A x .
012 TExcitation repetition time T R . × − sGradient strength G , Particle
Temperature T B
300 KSat. magnetization M C /TParticle core diameter D × − mParticle core volume V C / πD Particle hydrodynamic volume V H V C Dynamic viscosity (water) η . × − Pa sGyromagnetic ratio γ . × rad/sDamping parameter α . n e Anisotropy constant K anis Tracer
Mean lognormal distribution µ ln(13 × − )Standard deviation lognormal distribution σ We summarized several models relevant for magnetic particle imaging and showed simulations for one par-ticular special case to illustrate the different behavior of the models. These differences might be one of thereasons why model-based reconstructions with the equilibrium model are not of the same quality comparedto reconstructions with a measured linear operator. The relaxation effects are considered independent of eachother as it is commonly assumed that one cause for the relaxation dominates. However, combined modelsconsidering Brownian rotation and N´eel relaxation simultaneously are desirable and were investigated byusing their Langevin equations [46].Imaging quality can suffer when neglecting particle relaxation as it modifies the measured time signal.A loss of quality is then due to the spatial encoding in the time-dependent signal. A certain delay like inBrownian rotation and N´eel relaxation may cause shifted reconstructions in space when taking Cartesiantrajectories into account. The damping and smoothing observed in the Brownian rotation model may causean underestimated concentration and a spatially blurred reconstruction. The ill-posed nature of the problemallows for a certain degree of model errors only. In case of multidimensional trajectories, like Lissajoustrajectories, the rotation of the applied magnetic field vector has to be taken into account. Due to the lossof circular symmetry, the probability density function then needs to be approximated on the whole surfaceof the sphere. The higher dimensionality also increases the computational costs such that more efficient10 a) z = − . z = 0 mm(c) z = 2 . Figure 1: Simulated mean magnetic moment vector in e -direction considering the monodisperse models withparameters and particularly particle diameter specified in Table 1. The dotted vertical lines (red) highlightthe point in time with zero applied magnetic field. 11 a) z = − . z = 0 mm(c) z = 2 . Figure 2: Simulated mean magnetic moment vector in e -direction considering the polydisperse models withparameters and particularly the lognormal particle diameter distribution specified in Table 1. The dottedvertical lines (red) highlight the point in time with zero applied magnetic field.12umerical solutions are required to compute a solution of the probability density function. Furthermore,finding a direct and more efficient solution to compute the mean magnetic moment vector is highly desirable.Given the numerical solutions for the multidimensional case of the Brownian rotation and N´eel relaxation,these models still need to be physically validated for applied magnetic fields in MPI. For this purpose arecently proposed magnetic particle spectrometer (MPS) [4] can provide the required measurements. Theadvantage of the proposed MPS is that it allows applying a drive field with a three-dimensional excitationpattern and a constant offset field simulating the selection field at a fixed position.This overview about mathematical models for magnetic particle imaging builds the basis for severaldirections of future research. The numerical treatment of all models in the multidimensional case requiresfurther analysis and the development of efficient algorithms. A first step into the direction of theoreticalinvestigations was made by formulating a different problem setting motivated by the equilibrium model[39]. But by neglecting the temporal dependencies in the methodology, the equilibrium model as definedin this work is not directly covered. The mathematical models for MPI summarized in this work have notbeen investigated analytically including the related inverse problems. Besides the linear inverse problem ofreconstructing the concentration, several nonlinear inverse problems are motivated by applications and theparticle behavior itself. For example, joint concentration reconstructions combined with a spatial viscosityor temperature distribution were already motivated in the context of multi-color MPI [42, 51]. Interactionsbetween particles cause a nonlinear concentration dependence which becomes of interest [36]. Prior theconsideration of combined problems, a series of analytical works regarding the presented models is required. Acknowledgements
T. Kluth is supported by the Deutsche Forschungsgemeinschaft (DFG) within the framework of GRK 2224/1“Pi : Parameter Identification - Analysis, Algorithms, Applications”. The author also thanks Bangti Jin,University College London, for his helpful comments on this manuscript and fruitful discussions on MPI. References [1] C. Bathke, T. Kluth, C. Brandt, and P. Maass. Improved image reconstruction in magnetic particleimaging using structural a priori information.
International Journal on Magnetic Particle Imaging ,3(1):ID 1703015, 10 pages, 2017.[2] L. M. Bauer, S. F. Situ, M. A. Griswold, and A. C. Samia. Magnetic particle imaging tracers: state-of-the-art and future directions.
The Journal of Physical Chemistry Letters , 6(13), 2015.[3] W. F. Brown. Thermal fluctuations of a singledomain particle.
Journal of Applied Physics , 34(4):1319–1320, 1963.[4] X. Chen, M. Graeser, A. Behrends, A. von Gladiss, and T. M. Buzug. First measured result of the3d magnetic particle spectrometer. In L. Sefc, T. M. Buzug, and T. Knopp, editors,
InternationalWorkshop on Magnetic Particle Imaging , page 123. Infinite Science Publishing, 2017.[5] W. T. Coffey, P. J. Cregg, and Y. U. P. Kalmykov.
On the Theory of Debye and N´eel Relaxation ofSingle Domain Ferromagnetic Particles , pages 263–464. John Wiley & Sons, Inc., 1992.[6] W. T. Coffey and Y. P. Kalmykov. Thermal fluctuations of magnetic nanoparticles: Fifty years afterBrown.
Journal of Applied Physics , 112(12):121301, 2012.[7] L. R. Croft, P. W. Goodwill, and S. M. Conolly. Relaxation in x-space magnetic particle imaging.
IEEETransactions on Medical Imaging , 31(12):2335–2342, 2012.[8] R. J. Deissler, Y. Wu, and M. A. Martens. Dependence of Brownian and N´eel relaxation times onmagnetic field strength.
Medical Physics , 41(1):012301, 1–12, 2014.139] Y. Du, P. T. Lai, C. H. Leung, and P. W. T. Pong. Design of superparamagnetic nanoparticles formagnetic particle imaging (mpi).
International Journal of Molecular Sciences , 14(9):18682–18710, 2013.[10] D. Eberbeck, F. Wiekhorst, S. Wagner, and L. Trahms. How the size distribution of magnetic nanopar-ticles determines their magnetic particle imaging performance.
Applied Physics Letters , 98(18):182502,2011.[11] M. Gehre, T. Kluth, A. Lipponen, B. Jin, A. Sepp¨anen, J. P. Kaipio, and P. Maass. Sparsity recon-struction in electrical impedance tomography: an experimental evaluation.
Journal of Computationaland Applied Mathematics , 236(8):2126–2136, 2012.[12] M. Gehre, T. Kluth, C. Sebu, and P. Maass. Sparse 3D reconstructions in electrical impedance tomog-raphy using real data.
Inverse Problems in Science and Engineering , 22(1):31–44, 2014.[13] B. Gleich and J. Weizenecker. Tomographic imaging using the nonlinear response of magnetic particles.
Nature , 435(7046):1214–1217, 2005.[14] P. W. Goodwill and S. M. Conolly. The x-space formulation of the magnetic particle imaging process:1-D signal, resolution, bandwidth, SNR, SAR, and magnetostimulation.
IEEE Transactions on MedicalImaging , 29(11):1851–1859, 2010.[15] P. W. Goodwill and S. M. Conolly. Multidimensional X-space magnetic particle imaging.
IEEE Trans.Med. Imaging , 30:1581–1590, 2011.[16] M. Graeser, K. Bente, A. Neumann, and T. M. Buzug. Trajectory dependent particle response foranisotropic mono domain particles in magnetic particle imaging.
Journal of Physics D: Applied Physics ,49(4):045007, 2015.[17] M. Gr¨uttner, T. Knopp, J. Franke, M. Heidenreich, J. Rahmer, A. Halkola, C. Kaethner, J. Borgert,and T. M. Buzug. On the formulation of the image reconstruction problem in magnetic particle imaging.
Biomedical Engineering , 58(6):583–591, 2013.[18] J. Haegele, J. Rahmer, B. Gleich, J. Borgert, H. Wojtczyk, N. Panagiotopoulos, T. Buzug,J. Barkhausen, and F. Vogt. Magnetic particle imaging: visualization of instruments for cardiovas-cular intervention.
Radiology , 265(3):933–938, 2012.[19] B. Jin and P. Maass. Sparsity regularization for parameter identification problems.
Inverse Problems ,28(12):123001, 2012.[20] M. G. Kaul, O. Weber, U. Heinen, A. Reitmeier, T. Mummert, C. Jung, N. Raabe, T. Knopp, H. Ittrich,and G. Adam. Combined preclinical magnetic particle imaging and magnetic resonance imaging: Initialresults in mice.
Fortschr R¨ontgenstr , 187(05):347–352, 2015.[21] A. Khandhar, P. Keselman, S. Kemp, R. Ferguson, P. Goodwill, S. Conolly, and K. Krishnan. Evaluationof peg-coated iron oxide nanoparticles as blood pool tracers for preclinical magnetic particle imaging.
Nanoscale , 9(3):1299–1306, 2017.[22] A. P. Khandhar, R. M. Ferguson, H. Arami, and K. M. Krishnan. Monodisperse magnetite nanoparticletracers for invivo magnetic particle imaging.
Biomaterials , 34(15):3837 – 3845, 2013.[23] L. B. Kiss, J. Sderlund, G. A. Niklasson, and C. G. Granqvist. New approach to the origin of lognormalsize distributions of nanoparticles.
Nanotechnology , 10(1):25, 1999.[24] T. Kluth and P. Maass. Model uncertainty in magnetic particle imaging: Nonlinear problem formulationand model-based sparse reconstruction.
International Journal on Magnetic Particle Imaging , 3(2):ID1707004, 10 pages, 2017.[25] T. Knopp, S. Biederer, T. Sattel, J. Weizenecker, B. Gleich, J. Borgert, and T. M. Buzug. Trajectoryanalysis for magnetic particle imaging.
Physics in Medicine and Biology , 54(2):385–397, 2009.1426] T. Knopp, S. Biederer, T. F. Sattel, and T. M. Buzug. Singular value analysis for magnetic particleimaging. In
IEEE Nuclear Science Symposium Conference Record 2008 , pages 4525–4529, 2008.[27] T. Knopp, S. Biederer, T. F. Sattel, J. Rahmer, J. Weizenecker, B. Gleich, J. Borgert, and T. M. Buzug.2D model-based reconstruction for magnetic particle imaging.
Medical Physics , 37(2):485–491, 2010.[28] T. Knopp and T. M. Buzug.
Magnetic Particle Imaging: An Introduction to Imaging Principles andScanner Instrumentation . Springer, Berlin/Heidelberg, 2012.[29] T. Knopp, M. Erbe, T. F. Sattel, S. Biederer, and T. M. Buzug. A Fourier slice theorem for magneticparticle imaging using a field-free line.
Inverse Problems , 27(9):095004, 2011.[30] T. Knopp, N. Gdaniec, and M. M¨oddel. Magnetic particle imaging: From proof of principle to preclinicalapplications.
Physics in Medicine and Biology , 62(14):R124, 2017.[31] T. Knopp and M. Hofmann. Online reconstruction of 3D magnetic particle imaging data.
Physics inMedicine and Biology , 61(11):N257–67, 2016.[32] T. Knopp, J. Rahmer, T. F. Sattel, S. Biederer, J. Weizenecker, B. Gleich, J. Borgert, and T. M.Buzug. Weighted iterative reconstruction for magnetic particle imaging.
Physics in Medicine andBiology , 55(6):1577–1589, 2010.[33] T. Knopp, T. F. Sattel, S. Biederer, J. Rahmer, J. Weizenecker, B. Gleich, J. Borgert, and T. M.Buzug. Model-based reconstruction for magnetic particle imaging.
IEEE Transactions on MedicalImaging , 29(1):12–18, 2010.[34] J. Konkle, P. Goodwill, D. Hensley, R. Orendorff, M. Lustig, and S. Conolly. A convex formulation formagnetic particle imaging x-space reconstruction.
PLoS ONE , 10(10):e0140137, 2015.[35] J. Lampe, C. Bassoy, J. Rahmer, J. Weizenecker, H. Voss, B. Gleich, and J. Borgert. Fast reconstructionin magnetic particle imaging.
Physics in Medicine and Biology , 57(4):1113–1134, 2012.[36] N. L¨owa, P. Radon, O. Kosch, and F. Wiekhorst. Concentration dependent mpi tracer performance.
International Journal on Magnetic Particle Imaging , 2(1), 2016.[37] F. Ludwig, D. Eberbeck, N. L¨owa, U. Steinhoff, T. Wawrzik, M. Schilling, and L. Trahms. Character-ization of magnetic nanoparticle systems with respect to their magnetic particle imaging performance.
Biomedizinische Technik/Biomedical Engineering , 58(6):535–545, 2013.[38] M. Martens, R. Deissler, Y. Wu, L. Bauer, Z. Yao, R. Brown, and M. Griswold. Modeling the Brownianrelaxation of nanoparticle ferrofluids: Comparison with experiment.
Medical Physics , 40(2), 2013.[39] T. M¨arz and A. Weinmann. Model-based reconstruction for magnetic particle imaging in 2D and 3D.
Inverse Problems and Imaging , 10(4):1087–1110, 2016.[40] K. Murase, M. Aoki, N. Banura, K. Nishimoto, A. Mimura, T. Kuboyabu, and I. Yabata. Usefulness ofmagnetic particle imaging for predicting the therapeutic effect of magnetic hyperthermia.
Open Journalof Medical Imaging , 5(02):85, 2015.[41] L. N´eel. Thermoremanent magnetization of fine powders.
Rev. Mod. Phys. , 25:293–295, 1953.[42] J. Rahmer, A. Halkola, B. Gleich, I. Schmale, and J. Borgert. First experimental evidence of thefeasibility of multi-color magnetic particle imaging.
Physics in Medicine and Biology , 60(5):1775, 2015.[43] J. Rahmer, J. Weizenecker, B. Gleich, and J. Borgert. Signal encoding in magnetic particle imaging:properties of the system function.
BMC Medical Imaging , 9(4):1–21, 2009.[44] J. Rahmer, J. Weizenecker, B. Gleich, and J. Borgert. Analysis of a 3-D system function measured formagnetic particle imaging.
IEEE Transactions on Medical Imaging , 31(6):1289–1299, 2012.1545] D. B. Reeves and J. B. Weaver. Approaches for modeling magnetic nanoparticle dynamics.
CriticalReview in Biomedical Engineering , 42(1):85–93, 2014.[46] D. B. Reeves and J. B. Weaver. Combined N´eel and Brown rotational Langevin dynamics in magneticparticle imaging, sensing, and therapy.
Appl. Phys. Lett. , 107(22):223106, 2015.[47] H. Risken. Fokker-planck equation. In
The Fokker-Planck Equation , pages 63–95. Springer, 1996.[48] H. Rogge, M. Erbe, T. M. Buzug, and K. L¨udtke-Buzug. Simulation of the magnetization dynamics ofdiluted ferrofluids in medical applications.
Biomedizinische Technik/Biomedical Engineering , 58(6):601–609, 2013.[49] J. Salamon, M. Hofmann, C. Jung, M. G. Kaul, F. Werner, K. Them, R. Reimer, P. Nielsen, A. vomScheidt, G. Adam, T. Knopp, and H. Ittrich. Magnetic particle/magnetic resonance imaging: In-vitroMPI-guided real time catheter tracking and 4D angioplasty using a road map and blood pool tracerapproach.
PloS ONE , 11(6):e0156899–14, 2016.[50] S. A. Shah, D. B. Reeves, R. M. Ferguson, J. B. Weaver, and K. M. Krishnan. Mixed brownian alignmentand n´eel rotations in superparamagnetic iron oxide nanoparticle suspensions driven by an ac field.
Phys.Rev. B , 92:094438, Sep 2015.[51] C. Stehning, B. Gleich, and J. Rahmer. Simultaneous magnetic particle imaging (MPI) and temperaturemapping using multi-color MPI.
International Journal on Magnetic Particle Imaging , 2(2):ID 1612001,6 pages, 2016.[52] E. Stoner and E. Wohlfarth. A mechanism of magnetic hysteresis in heterogeneous alloys.
IEEETransactions on Magnetics , 27(4):3475–3518, 1991.[53] M. Storath, C. Brandt, M. Hofmann, T. Knopp, J. Salamon, A. Weber, and A. Weinmann. Edgepreserving and noise reducing reconstruction for magnetic particle imaging.
IEEE Transactions onMedical Imaging , 36(1):74–85, 2017.[54] C. Tannous and J. Gieraltowski. The Stoner-Wohlfarth model of ferromagnetism.
European Journal ofPhysics , 29(3):475–487, 2008.[55] K. Them, M. G. Kaul, C. Jung, M. Hofmann, T. Mummert, F. Werner, and T. Knopp. Sensitivityenhancement in magnetic particle imaging by background subtraction.
IEEE Trans Med Imaging ,35(3):893–900, 2016.[56] M. Utkur, Y. Muslu, and E. U. Saritas. Relaxation-based viscosity mapping for magnetic particleimaging.
Physics in Medicine and Biology , 62(9):3422, 2017.[57] A. Weber, F. Werner, J. Weizenecker, T. M. Buzug, and T. Knopp. Artifact free reconstruction with thesystem matrix approach by overscanning the field-free-point trajectory in magnetic particle imaging.
Physics in Medicine and Biology , 61(2):475–487, 2015.[58] J. Weizenecker, B. Gleich, and J. Borgert. Magnetic particle imaging using a field free line.
Journal ofPhysics D: Applied Physics , 41(10):105009, 2008.[59] J. Weizenecker, B. Gleich, J. Rahmer, and J. Borgert. Micro-magnetic simulation study on the magneticparticle imaging performance of anisotropic mono-domain particles.
Physics in Medicine and Biology ,57(22):7317, 2012.[60] J. Weizenecker, B. Gleich, J. Rahmer, H. Dahnke, and J. Borgert. Three-dimensional real-time in vivomagnetic particle imaging.
Physics in medicine and biology , 54(5):L1, 2009.[61] E. Y. Yu, M. Bishop, B. Zheng, R. M. Ferguson, A. P. Khandhar, S. J. Kemp, K. M. Krishnan, P. W.Goodwill, and S. M. Conolly. Magnetic particle imaging: A novel in vivo imaging platform for cancerdetection.
Nano Letters , 17(3):1648–1654, 2017.1662] B. Zheng, T. Vazin, P. W. Goodwill, A. Conway, A. Verma, E. Ulku Saritas, D. Schaffer, and S. M.Conolly. Magnetic particle imaging tracks the long-term fate of in vivo neural cell implants with highimage contrast.
Scientific Reports , 5:14055, 9 pages, 2015.
A Fokker-Planck equation for Brownian rotation
We derive (7) from the Langevin equation of the particle dynamics. As the spatial dependence in theprobability density function f solely results implicitly from the spatially dependent applied magnetic field B app , we consider the problem of determining f for one fixed x ∈ Ω, and further omit the spatial variable x .The Langevin equation for Brownian rotation of a single particle with magnetic moment vector ˜ m reads ∂∂t ˜ m = 16 V H η (cid:16) ˜ m × B eff + ˜ D Γ (cid:17) × ˜ m, (14)with ˜ D > (cid:104) Γ i ( t ) (cid:105) = 0, (cid:104) Γ i ( t )Γ j ( t ) (cid:105) = δ ij δ ( t − t ), for all t, t , t > i, j = 1 , , δ i,j is the Kronecker delta and δ is the Dirac deltadistribution. Substituting m = ˜ m/m , m = | ˜ m | yields ∂∂t m = 16 V H η (cid:16) m × m B eff + ˜ D Γ (cid:17) × m = 16 V H η (cid:16) m ( m × B eff ) × m + ˜ D Λ( m )Γ (cid:17) (15)where Λ( m ) = − m m m − m − m m . (16)The Fokker Planck equation for the probability density function f : S × (0 , T ) → R + thus becomes [47] ∂∂t f ( m, t ) = − div m (cid:18) a ( m, t ) f ( m, t ) + 12 B ( m, t ) (cid:16) div m ( B ( j ) ( m, t ) f ( m, t )) (cid:17) j =1 , , (cid:19) (17)where a : S × (0 , T ) → R and B : S × (0 , T ) → R × with B ( j ) being the j -th column of B . From theLangevin equation it follows a ( m, t ) = m V H η ( B eff − ( m · B eff ) m ) (18)and B ( m, t ) = ˜ D V H η Λ( m ) T . (19)By using the fact that m · ∇ m f = 0 (as the gradient is tangent to the unit sphere) we obtain by usingΛ( m ) T Λ( m ) ∇ m f = ( m × ∇ m f ) × m = ∇ m f the desired Focker Planck equation ∂∂t f ( m, t ) = − V H η div m (cid:32) m ( B eff − ( m · B eff ) m ) f ( m, t ) − ˜ D V H η ∇ m f ( m, t ) (cid:33) (20)The diffusion coefficient ˜ D is determined by considering the equilibrium case [48], i.e. let t ∈ I be the pointin time such that ∂∂t f ( m, t ) = 0. We further assume that f : S → R + ∪{ } with f ( · ) = f ( · , t ) correspondsto the Boltzmann distribution for Brownian rotation in the equilibrium, i.e., f ( m ) = ke − β H ( m,t ) , where H is the Hamiltonian with ∇ m H = − B eff and β = m κ B T B . Assume B eff does not depend on m , which holds if B eff = B app . We use H ( m, t ) = − B app · m and ∇ m f · m = 0 to obtain0 = div m (cid:32) m ( B app − ( m · B app ) m ) e m κ B T B B app · m − ˜ D V H η m κ B T B B eff e m κ B T B B app · m (cid:33) = m κ B T B | B app | − ˜ D V H η (cid:18) m κ B T B (cid:19) | B app | . (21)17rom this it follows ˜ D = √ V H ηκ B T B . Defining τ = V H ηκ B T B yields ∂∂t f ( m, t ) = − τ div m (cid:18) m κ B T B ( B app − ( m · B app ) m ) f ( m, t ) − ∇ m f ( m, t ) (cid:19) . (22) B Derivation of Langevin function
In the following we give an example how the equilibrium model and particularly using the Langevin function ismotivated and related to Brownian rotation. As can be seen in A, the Fokker-Planck equation is parametrizedsuch that the probability density function of the magnetic moment vector in equilibrium at time t ∈ I is ofthe form f ( m ) = ke − β H ( m,t ) . (23)Considering the mean magnetic moment vector yields¯ m ( t ) = m (cid:90) S mf ( m ) dm = m k (cid:90) S me β ( R T m ) · ( R T B app ( t )) dm = m kR (cid:90) S ye y βb dy = m kR (cid:90) π (cid:90) π sin( θ ) cos( φ )sin( θ ) sin( φ )cos( θ ) | sin( θ ) | e cos( θ ) βb dφ dθ = 2 πm kRe (cid:90) π cos( θ ) sin( θ ) e cos( θ ) βb dθ = − πm kRe (cid:90) − xe xβb dx = − πm kRe (cid:32)(cid:20) βb xe xβb (cid:21) − − (cid:90) − βb e xβb dx (cid:33) = − Re m (cid:18) e βb + e − βb e βb − e − βb − βb (cid:19) (24)where R ∈ R × is a rotation matrix sucht that R T B app ( t ) = e b , where b = | B app ( t ) | , and where k isobtained by similar transformations such that it is given by1 k = (cid:90) S e − β H ( m,t ) dm = (cid:90) S e y βb dy = 2 π (cid:90) − e xβb dx = 2 πβb ( e βb − e − βb ) . (25)The length of the mean magnetic moment is thus given by | ¯ m ( t ) | = L β ( | B app ( t ) | ) with L β given by (4). C Fokker-Planck equation for N´eel relaxation
We derive (9) from the Langevin equation based on the Landau-Lifshitz-Gilbert equation for the particledynamics. As the spatial dependence in the probability density function f solely results implicitly from thespatial dependent applied magnetic field B app we consider the problem of determining f for one fixed x ∈ Ω,and further omit the spatial variable x . The Langevin equation for N´eel relaxation of a single particle withmagnetic moment vector ˜ m is obtained from the Landau-Lifshitz-Gilbert equation and reads ∂∂t ˜ m = ˜ γ (cid:18) ( B eff + ˜ D Γ) × ˜ m + αm ( ˜ m × ( B eff + ˜ D Γ)) × ˜ m (cid:19) , (26)with ˜ D > γ = γ/ (1 + α ) and where Γ is a white noise component with (cid:104) Γ i ( t ) (cid:105) = 0, (cid:104) Γ i ( t )Γ j ( t ) (cid:105) = δ ij δ ( t − t ), for all t, t , t > i, j = 1 , , δ i,j is the Kronecker delta and18 is the Dirac delta distribution. Substituting m = ˜ m/m , m = | ˜ m | yields ∂∂t m = ˜ γ (cid:16) ( B eff + ˜ D Γ) × m + α ( m × ( B eff + ˜ D Γ)) × m (cid:17) = ˜ γ (cid:16) B eff × m + α ( m × B eff ) × m + ˜ D (Λ( m ) T + α Λ( m ) T Λ( m ))Γ (cid:17) (27)where Λ( m ) = − m m m − m − m m (28) ∧ Λ( m ) T Λ( m ) = − m + m − m m − m m − m m m + m − m m − m m − m m m + m . (29)The Fokker-Planck equation for the probability density function f : S × (0 , T ) → R + ∪ { } thus becomes [47] ∂∂t f ( m, t ) = − div m (cid:18) a ( m, t ) f ( m, t ) + 12 B ( m, t ) (cid:16) div m ( B ( j ) ( m, t ) f ( m, t )) (cid:17) j =1 , , (cid:19) (30)where a : S × (0 , T ) → R and B : S × (0 , T ) → R × with B ( j ) being the j -th column of B . From theLangevin equation it follows a ( m, t ) = ˜ γ ( B eff × m + α ( m × B eff ) × m ) (31)and B ( m, t ) = ˜ D ˜ γ (cid:0) Λ( m ) T + α Λ( m ) T Λ( m ) (cid:1) . (32)By using the fact that m · ∇ m f = 0 (as the gradient is tangent to the unit sphere) we obtain by usingΛ( m ) T Λ( m ) ∇ m f = ( m × ∇ m f ) × m = ∇ m f the desired Focker Planck equation ∂∂t f ( m, t )= − ˜ γ div m (cid:32) ( B eff × m + α ( B eff − ( m · B eff ) m )) f ( m, t ) − ˜ D ˜ γ (1 + α )2 ∇ m f ( m, t ) (cid:33) (33)The diffusion coefficient ˜ D is determined by considering the equilibrium case [48], i.e. let t ∈ I be the pointin time such that ∂∂t f ( m, t ) = 0. We further assume that f : S → R + ∪{ } with f ( · ) = f ( · , t ) correspondsto the Boltzmann distribution for Brownian rotation in the equilibrium, i.e., f = ke − β H ( m,t ) , where H isthe Hamiltonian with ∇ m H = − B eff and β = m κ B T B . We use ∇ m f = βB eff f and ∇ m f · m = 0 to obtain0 = div m (cid:32) ( B eff × m + α ( B eff − ( m · B eff ) m )) f − ˜ D ˜ γ (1 + α )2 ∇ m f (cid:33) = div m (( B eff × m ) f ) + α div m ( B eff f ) − ˜ D ˜ γ (1 + α )2 ∆ m f = div m ( B eff × m ) f + (cid:32) αβ − ˜ D ˜ γ (1 + α )2 (cid:33) ∆ m f = (cid:32) αβ − ˜ D ˜ γ (1 + α )2 (cid:33) ∆ m f (34)where the last inequality holds as div m ( B eff × m ) = 0 for B eff = B app + B anis . It follows ˜ D = (cid:113) ακ B T B ˜ γ (1+ α ) m .Defining τ = m α ˜ γκ B T B yields ∂∂t f ( m, t ) = − div m (cid:18) ˜ γ ( B eff × m + α ( B eff − ( m · B eff ) m )) f ( m, t ) − τ ∇ m f ( m, t ) (cid:19) . (35)19 Numerical solution for 1D excitation
The following numerical solution follows [8]. We compute the mean magnetic moment vector ¯ m for anapplied magnetic field of the form B app ( x, t ) = b ( x, t ) e with a given field b : Ω × I → R . We additionallyassume that Ω ⊂ { te | t ∈ R } .Under these assumptions the Fokker-Planck equations for Brownian rotation in (7) and N´eel relaxationin (9) can be reformulated in terms of the angle θ between m and e . Let x ∈ Ω be fixed. By assuming z = cos( θ ) the resulting ordinary differential equations for ˜ f x : [ − , × I → R + ∪ { } then become forBrownian rotation 2 τ B ∂∂t ˜ f x ( z, t ) = ∂∂z (cid:18) (1 − z ) (cid:18) ∂∂z ˜ f x ( z, t ) − βb ( x, t ) ˜ f x ( z, t ) (cid:19)(cid:19) (36)with physical parameters τ B = 3 V H η/ ( κ B T B ) and β = m κ B T B . For N´eel relaxation we obtain2 τ N ∂∂t ˜ f x ( z, t ) = ∂∂z (cid:18) (1 − z ) (cid:18) ∂∂z ˜ f x ( z, t ) − βb ( x, t ) ˜ f x ( z, t ) − ˜ kz ˜ f x ( z, t ) (cid:19)(cid:19) (37)with physical parameters τ N = m α ˜ γκ B T B , β = m κ B T B , and ˜ k = K anis m M C κ B T B . The mean magnetic moment vector¯ m is then obtained by ¯ m ( x, t ) = m (cid:90) − z ˜ f x ( z, t ) dze . (38)We expand the probability density function ˜ f in Legendre polynomials { P n } n ∈ N and a n : I → R , n ∈ N , i.e.,˜ f ( z, t ) = ∞ (cid:88) n =0 a n ( t ) P n ( z ) . (39)Using this relation and the constraint (cid:82) − ˜ f ( z, t ) dz = 1, t ∈ I , we obtain the system of differential equationsfor Brownian rotation a ( t ) = 122 τ B n ( n + 1) ddt a n ( t ) = − a n ( t ) + βb ( x, t ) (cid:18) a n − ( t )2 n − − a n +1 n + 3 (cid:19) n ≥ a ( t ) = 122 τ N n ( n + 1) ddt a n ( t ) = − a n ( t ) + βb ( x, t ) (cid:18) a n − ( t )2 n − − a n +1 n + 3 (cid:19) + ˜ K (cid:18) ( n − a n − ( t )(2 n − n −
1) + na n ( t )(2 n − n + 1) − ( n + 1) a n ( t )(2 n + 1)(2 n + 3) − ( n + 2) a n +2 ( t )(2 n + 3)(2 n + 5) (cid:19) n ≥ a − ( t ) = 0. Assuming a uniform distribution as initial condition results in the initialvalues a (0) = 1 / a n (0) = 0, n ≥
1, for both systems. Using the relation (cid:82) − z ˜ f x ( z, t ) dz = 2 / a ( t )determines the mean magnetic moment vector. The time derivative of the mean magnetic moment vector isthus given by the respective differential equation for n = 1.For the simulation we use an approximated probability density function for N ∈ N given by˜ f ( z, t ) = N (cid:88) n =0 a n ( t ) P n ( z ) . (42)The applied magnetic field is given by b ( x, t ) = A cos(2 πf t ) + G , x . (43)where A, f > G3
1, for both systems. Using the relation (cid:82) − z ˜ f x ( z, t ) dz = 2 / a ( t )determines the mean magnetic moment vector. The time derivative of the mean magnetic moment vector isthus given by the respective differential equation for n = 1.For the simulation we use an approximated probability density function for N ∈ N given by˜ f ( z, t ) = N (cid:88) n =0 a n ( t ) P n ( z ) . (42)The applied magnetic field is given by b ( x, t ) = A cos(2 πf t ) + G , x . (43)where A, f > G3 ,3