Generally, the distortion of the magnetic field is connected to the susceptibility difference Δχ = χ_i − χ_e between the magnetized object and its surrounding medium and leads to a variation in the local Larmor frequency (1)ω(r)=δω∂2∂z2∫Gd3r′|r−r′| , with frequency shift δω = γB₀Δχ and gyromagnetic ratio γ = 2.675 × 10⁸s⁻¹T⁻¹ [34]. Let p(r, r₀, t) be the transition probability of a proton spin at point r and time t that started at position r₀ and time t = 0. Since, presumably, spin diffusion inside the dephasing volume encounters no barriers, transition probability p(r, r₀, t) is nothing less than the Green's function of the diffusion equation (2)∂∂tp(r,r0,t)=DΔp(r,r0,t), where Δ denotes the Laplace operator and D is the diffusion coefficient. Thus, (3)p(r,r0,t)=etDΔδ(r−r0), and it is evident that p(r, r₀, t) only depends on the spatial difference r − r₀.

Then, to investigate stochastic field fluctuations around object G, a two-point correlation function K(t) can be introduced (4)K(t)=∫Vd3r∫Vd3r0ω(r)p(r,r0,t)p(r0)ω(r0) (5)=1V∫Vd3rω(r)etDΔω(r), where we made use of Equation (3) in the last step. The probability density of the equilibrium p(r₀) = 1/V due to the assumption of homogeneous spin density. The frequency autocorrelation function K(t) is linked to the measured MR-signal S(t) and can be determined through measurement of asymmetric spin echos (see also Figure 1B and [18]). An asymmetric spin echo S(t)|_tp is a spin echo that is measured after a 180° refocusing pulse given at a time t_p ≠ t/2 where a refocusing pulse given at time t/2 corresponds to a symmetric spin echo S(t)|_t/2. When |t_p − t/2| = t_s, the following equality in our notation holds for low orders of t_s according to Jensen et al. [18]: (6)S(t)|ts + t/2=S(t)|t/2exp(−2ts2K(t2)). Usually, Equation (6) is used in practice to be fitted to MRI data to obtain the correlation function.

In this work, two kinds of magnetic perturbers are considered: spheres and cylinders (see Figures 2A,B, respectively). For instance, an ensemble of microscopic spherical magnetic objects can represent the intracellular accumulation of iron particles in neurodegenerative disease [35] or the accumulation of ultrasmall paramagnetic iron oxide as a contrast agent to track macrophages in nerve injury or atherosclerotic disease [4, 5]. Collectives of cylindrical magnetic perturbers can be represented through a microvascular network in skeletal or cardiac muscle tissue [36, 37] where capillaries contain paramagnetic deoxygenated hemoglobine [7].

The magnetic field around a spherical magnetic perturber is that of a magnetic dipole, i.e., with the external magnetic field B₀ being parallel to the z-axis of the coordinate system and utilization of spherical coordinates (r, θ, φ) (see also Figure 2A), the Larmor frequency reads (7)ω(r)=δωRi33cos2(θ)−1r3. Here, δω=13γΔχB0 represents the frequency shift on the equatorial surface of the sphere with radius R_i. Its dephasing volume V ranges from the sphere surfaces to the surface of the mean relaxation volume. The latter depends on the form of the magnetized object and, therefore, the outer surface of V can be considered as a sphere with radius R. Consequently, the volume fraction of the magnetic material η=Ri3/R3.

For a cylinder, the angle θ between its axis and the external magnetic field B₀ has to be acknowledged in the frequency shift to obtain the local frequency (see Figure 2B) (8)ω(r)=δωRi2cos(2φ)r2, where δω=12γΔχB0sin2(θ) is the external frequency shift. In analogy to spheres, the local volume fraction η=Ri2/R2 for a outer cylindrical surface of the dephasing volume with radius R. In both cases, the rationale of the approach has already been discussed in detail elsewhere [13, 38, 39].

For ensembles of several (microscopic) magnetic objects, it is helpful to avail oneself of the principle of supply areas as in the Krogh capillary model [40]. If the volume hosting the magnetic perturbers can be split into smaller volumes of similar geometry and if the local resonance frequency around the magnetic perturber is point symmetric, i.e., ω(r) = ω(−r) which is valid for both spheres and cylinders, it is possible to think of spin trajectories that leave their initial supply area or dephasing volume as being reflected at its boundary. Also, the same applies for spins whose trajectories hit the surface of the magnetic perturber that is assumed impermeable. Then, the Krogh model allows reducing the problem to one single object, i.e., diffusion is not considered in the whole tissue but, instead, for one single magnetic perturber only and reflective boundary conditions are assumed (9)∂∂rp(r,r0,t)|r = Ri=0=∂∂rp(r,r0,t)|r = R. The case of unrestricted diffusion for a dephasing volume V extending from the surface of the magnetic perturber to infinity has already been analyzed in detail for spheres [26, 27, 41] and cylinders [27] and shall not be the subject of this investigation. It corresponds to very small volume fractions η ≈ G/V_Voxel where spin diffusion around the magnetic object does not reach the outer surface of V during the measurement. A detailed account in a similar notation can be found in Ziener et al. [25].

The diffusion equation (Equation 2) can be solved with the help of an eigenfunction expansion of the transition probability p(r, r₀, t): (10)p(r,r0,t)=∑κn≥0∞e−κn2DtRi2ϕn(r)ϕn*(r0)             =1V+∑κn>0∞e−κn2DtRi2ϕn(r)ϕn*(r0). The eigenfunctions ϕ_n are solutions of the eigenvalue equation (11)Δϕn(r)=−κn2Ri2ϕn(r), with eigenvalues κ_n, and, they fulfill the orthogonality condition (12)∫d3rϕn(r)ϕn′*(r)=δnn′, where R_i is the radius of the inner sphere or cylinder, respectively. Thereby, the lowest eigenvalue κ₀ = 0 corresponds to the lowest eigenfunction ϕ0=1/V [42]. The influence of diffusion is taken into account by the correlation time τ (13)τ=Ri2D. By substituting the spectral expansion of the transition probability p(r, r₀, t) into the definition of the autocorrelation function (Equation 4), we obtain (14)K(t)=∑n = 1∞Fn2e−κn2tτ with expansion coefficients (15)Fn=1V∫d3rω(r)ϕn(r). Since ϕ0(r)=1/V, the expansion coefficient F₀ = 0 does not contribute to the correlation function.

Explicit expressions can be given for both spherical and cylindrical geometries; in the case of spheres, Fn2 reads (16)Fn2=2165κn2ηη−1[j2′(κnη3)−η43j2′(κn)]2[κn2−6] [j2′(κnη3)]2+η13[6η23−κn2] [j2′(κn)]2, with κ_n obeying the eigenvalue equation (17)j2′(κn)y2′ (κnη3)=j2′ (κnη3) y2′(κn). Functions j2′ and y2′ are the first derivatives of the spherical Bessel functions of the first and second kind, respectively [25]. An approximation for the first zero κ₁ is provided in Appendix A.

For cylinders, we have with cylindrical Bessel functions J₂ and Y₂ of the first and second kind, respectively, (18)Fn2=8κn2ηη−1[J2′(κnη)−η32J2′(κn)]2[κn2−4] [J2′(κnη)]2+[4η−κn2] [J2′(κn)]2, and eigenvalues κ_n are the solution of (19)J2′(κn)Y2′ (κnη)=J2′ (κnη) Y2′(κn). The long time behavior of K(t) can be determined by considering the first eigenvalue κ₁ only: (20)KL(t)=F12e−κ12tτ. Evidently, K_L(t) converges asymptotically against zero for large t or K(∞) = 0. On the other hand, the value of the correlation function at t = 0 can be obtained by making use of Equation (3) with the initial condition p(r, r₀, 0) = δ(r − r₀): (21)K(0)=〈ω2(r)〉=1V∫Vd3r ω2(r) (22)=∑n = 1∞Fn2. The last relation is useful for numerical computation and corresponds to Parseval's theorem in Fourier analysis.

The generalized moment expansion method or generalized moment approximation (GMA) has been used recently as an effective algorithm to approximate the time dependence of observables that are connected to reactive and non-reactive processes involving Brownian motion [29, 30, 43, 44]. It can be applied for monitoring diffusive redistributions and barrier crossing problems as well as to approximate autocorrelation functions [30]. In our case, the generalized moment analysis will be used to derive analytical expressions for the correlation function for biexponentially approximated decay forms that will prove to be more accurate than monoexponential approximations. The following subsections will introduce generalized moments and the Padé approximation and results from monoexponential approximations, as detailed in Ziener et al. [25], will be briefly recapitulated to give a buildup to biexponential approximations.

In order to solve differential equation (Equation 30) for spherical magnetic objects, it is helpful to start with the ansatz (41)μ−n(r)=δωRi3μ−n(r)[3 cos2(θ)−1] that separates the angle dependent part of the local Larmor frequency given in Equation (7). Also, the radial part of the Laplace operator can be written as (42)Δd,r=1rd−1ddrrd−1ddr with d = 3 for spheres and d = 2 for infinitely long cylinders. This leads to an iterative set of equations for the functions μ_−n(r) of the form (43)[Δ3,r−6r2] μ−n(r)=−μ−(n−1)(r), or, equivalently, (44)μ−n(r)=[6r2−Δ3,r]−n1r3. The differential equation (Equation 43) can be transformed to an Euler type differential equation. Using the general expression (Equation 41) for the solution of the differential equation (Equation 30), where the local Larmor frequency corresponds to Equation (7), the low-frequency moments given in Equation (28) can be found as (45)μ−n=125δω2Dnη1−ηRi3∫RiRdrμ−n(r)r. The evaluation of the integral yields the general form (46)μn=δω2τnln(η). For spherical geometries, the functions μ_n(r) can be found as the solution of Equation (43). The recursion starts with μ₀(r) which results from Equation (44). The functions are given by

We have used the abbreviation (48)di=∑j=0iηj3 and the constants c₁, …, c₆ that are given by (49)c1=-d2d4,c2=1+c1,c3=25126η13d12d42+1136d4,c4=c3-328-25126d2d4,c5=61588[η-13-d3d4]-62515876η-2d1d43-531752331d2+250d1+232d42,c6=c5+61588+62515876η13d12d42+531752385-54d2d4, where η=Ri3/R3. The functions l_n(η) are necessary to determine the moments μ_n in Equation (46) for negative indices from Equation (45). They are given by

It is convenient to define a ratio τ_n of successive moments as (51)τn=|μnμn+1|. As in Ziener et al. [25], we make use of Koenig's theorem [48, 49] on the zeros of Padé approximants [48, 49]: the theorem states that the series 1/τ_n converges to the lowest eigenvalue of Equation (11) and [39], since then (52)κ12=τlimn→∞μ-nμ-(n+1), where τ is the correlation time from Equation (13). The lowest eigenvalue κ₁ is the solution to the geometry-dependent eigenvalue equation, i.e., Equation (17) for spheres and Equation (19) for cylinders. In Figure 3, a sequence of successive moment ratios τ_n ranging from n = −3, …, 0 is compared with the exact numerical solution for κ₁ for spherical magnetic objects and Koenig's theorem is nicely exemplified for η → 1. The asymptotic value of the first eigenvalue for large volume fractions could be determined as limη→1κ12(η)=6 for spheres and limη→1κ12(η)=4 for cylinders (from Equations [A14] and [A5] in Ziener et al. [25], respectively), simply stating that, in this limit, the correlation function decays monoexponentially with time constant κ12. An approximation for κ₁ is provided in Appendix A. Figure 3 also validifies the concept of generalized moments to define a search interval for numerical computation of the exact solution from Equations (17) and (19). However, the successive moments curves diverge from the exact solution for small volume fractions.

Different cases of the monoexponential functions K_i(t) (see Appendix B and Equation B3) are being shown in Figure 4A to demonstrate the dependence of the monoexponential approximation on the chosen moments μ_n. Evidently, the initial value of K_i(t), i ≥ 1, decreases with increasing i while the exponential cutoff begins at longer times for increasing i. Also, the approximants K_i(t) converge to the long time behavior given in Equation (20), as expected from Koenig's theorem. To evaluate the accuracy of the biexponential approximation, a double logarithmic plot has been chosen in Figure 4B. The (2, 2)-GMA with the two lowest high and low frequency moments approximates the short as well as the long time limit sufficiently well. In contrast, the (1, 3)-GMA with the three lowest low and the first high frequency moment approximates the long time limit significantly better than the short time limit. An application to biological tissue is provided in Figure 4C where mono- and bi-exponential approximations in the GMA of the correlation function are visualized for values that are typical for susceptibility gradients based on nonheme iron-storage in brain tissue (through iron-rich oligodendrocytes in the Globus pallidus that approximately have radii R_i of 3–4 μm, a frequency shift of δω = 184 s⁻¹, a tissue density of η = 0.04 and D = 0.76·10⁻⁹m²s⁻¹, leading to τ = 16.12 ms; see [50, 51] for more details). It can be seen that the biexponential approximations K_{(1, 3)}(t) and K_{(2, 2)}(t) are in excellent agreement with the actual correlation function K(t), while the monoexponential functions K_i(t) only approximate K(t) sufficiently for t ≳ 30 ms. However, changes in the size R_i of the spherical perturbers affect the correlation time τ as well according to Equation (13). The correlation function in dependence on correlation time τ is therefore shown in Figure 4D at the characteristic time t = 10 ms for the same set of mono- and bi-exponential approximations as in Figure 4C. In this case, the biexponential approximations show an equally well agreement with the numerically obtained correlation function K(τ) for all correlation times, whereas the monoexponential approximations exhibit significant deviations from K(τ) for correlation times that range from 40 ms to 5 s. For very short and very long correlation times, however, the monoexponential approximations K₀(τ) and K₁(τ) coincide with K(τ). Mono- and bi-exponential approximations as compared to K(t) for nanoparticle experiments based on parameters in Muller et al. [52], and for lung tissue based on parameters given in Stone et al. [53] and Baete et al. [54], are shown in Appendix C.

In analogy to spherical magnetic objects, the correlation functions for spins that diffuse between two concentric cylinders can be obtained within the GMA by employing a similar ansatz as the one for a spherical geometry and by using Equation (42) for d = 2 as shown in Ziener et al. [25]. Then, generalized moments are determined by the differential equation (53)μ-n(r)=[4r2-Δ2,r]-n1r2. Again, for n = 0, 1, 2, 3, the functions μ_−n(r) for cylindrical geometries that obey Equation (53) can be determined as

and the constants c₁, …, c₆ are given by (55)c1=12ln(η)1-η2,c2=c1-12,c3=-3c1,c4=67[c1-1],c5=11+η[2c1+116η],c6=η1+η[116[1+1η+1η2]-2c1]. The functions l_n(η) are obtained in analogy to the spherical case and are given as (56)l1(η)=2η[1+η],l0(η)=η2,l-1(η)=18ηln(η)η-1,l-2(η)=132,l-3(η)=1256η1-η[341-η2η2+[ln(η)]21-η2]. As in the case for spheres, the low-frequency moments μ_−n for cylinders from Equation (28) can be expressed as (57)μ-n=δω2Dnη1-ηRi2∫RiRdrμ-n(r)r, that correspond to Equation (77) in Ziener et al. [25]. Evaluation of the integral leads to Equation (46), but with the respective values of δω, τ, and l_n(η) for cylinders (see Equation 56).

In analogy to Figure 4A, the time evolution of the functions K_i(t) and their proximity to the long time approximation K_L(t) from Equation (20) is shown in Figure 5A. As in the spherical case, the initial value of K_i(t), i ≥ 1, decreases with increasing i and successive moments converge to the long time approximation K_L(t) for d = 2.

In comparison, the biexponential approximation is shown in Figure 5B: for the (2, 2)- and the (1, 3)-GMA it can be seen that the first approximation describes short and long time limit sufficiently well whereas the second approximation is more accurate in the long time limit. An application to biological tissue is provided in Figure 5C where mono- and bi-exponential approximations in the GMA of the correlation function are visualized for values that are typical for susceptibility gradients between capillaries, that contain deoxygenated hemoglobine, and the surrounding tissue [with approximate values R_i = 4 μm [55], δω = 200 s⁻¹ at 3 Tesla [12, 56], a tissue density of η = 0.03 [55], and D = 1·10⁻⁹m²s⁻¹ [57], leading to τ = 16 ms]. It can be seen that the biexponential approximations K_{(1, 3)}(t) and K_{(2, 2)}(t) are in excellent agreement with the actual correlation function K(t) for t ≳ 20 ms, while only the monoexponential functions K₂(t) and K₃(t) sufficiently approximate K(t) for t ≳ 80 ms. Still, the deviation of the biexponential approximations to K(t) for t < 30 ms are visibly smaller than for the monoexponential approximations. However, changes in the size R_i of the cylindrical perturbers affect the correlation time τ as well according to Equation (13). The correlation function in dependence on correlation time τ is therefore shown in Figure 5D at the characteristic time t = 10 ms for the same set of mono- and bi-exponential approximations as in Figure 5C. In this case, the biexponential approximations show an equally well agreement with the numerically obtained correlation function K(τ), especially for correlation times τ ≲ 50 ms, whereas the monoexponential approximations exhibit significant deviations from K(τ) for correlation times that range from 1 ms to 5 s. For very short and very long correlation times, however, the monoexponential approximations K₀(τ) and K₁(τ) coincide with K(τ), as in the case of spherical perturbers.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.