Localization regime in diffusion NMR: theory and experiments

TL;DR

This paper explores the localization regime in diffusion NMR, demonstrating how high-gradient conditions lead to boundary-localized magnetization and a stretched-exponential signal decay, with theoretical and experimental validation.

Contribution

It provides a comprehensive theoretical and experimental analysis of the localization regime in diffusion NMR across various geometries, highlighting its generic nature at high gradients.

Findings

Localization regime causes boundary-localized magnetization.

Stretched-exponential decay observed in experiments.

Agreement between theory, simulations, and experiments.

Abstract

In this work we investigate the emergence of the localization regime for diffusion in various geometries: inside slabs, inside cylinders and around rods arranged on a square array. At high gradients, the transverse magnetization is strongly attenuated in the bulk, whereas the macroscopic signal is formed by the remaining magnetization localized near boundaries of the sample. As a consequence, the signal is particularly sensitive to the microstructure. Our theoretical analysis relies on recent mathematical advances on the study of the Bloch-Torrey equation. Experiments were conducted with hyperpolarized xenon-129 gas in 3D-printed phantoms and show an excellent agreement with numerical simulations and theoretical predictions. Our mathematical arguments and experimental evidence indicate that the localization regime with a stretched-exponential decay of the macroscopic signal is a generic feature of diffusion NMR that can be observed at moderately high gradients in most NMR scanners.

Figures

Tables

Table 1. Table 1: Characteristics of the phantoms (see Fig. 2 ).

#	Geometry	Dimensions
1	slab	$L=3\mathrm{mm}$
2	cylinders	$2R=3.8\mathrm{mm}$
3	cylinders	$2R=2\mathrm{mm}$
4	array of rods	$17\times 10$ rods
		$2R=3.2\mathrm{mm}$; $L=4\mathrm{mm}$
		$b_V=1.48\mathrm{mm}$; $b_H=1.35\mathrm{mm}$
5	array of rods	$20\times 10$ rods
		$2R=3.2\mathrm{mm}$; $L=3.4\mathrm{mm}$
		$b_V=0.78\mathrm{mm}$; $b_H=1.08\mathrm{mm}$

Equations

$$\frac{\partial M}{\partial t}=(D_0{\nabla }^2-i\gamma gx)M,$$

$${\int }_{\Omega }{v}_n(\mathbf{r})v_m(\mathbf{r})\mathrm{d}\mathbf{r}=\{\begin{array}{l}0\text{ if }n\ne m\\ 1\text{ if }n=m\end{array}.$$

$$M(\mathbf{r},t)=\frac{1}{\sqrt{V}}\sum _{n\ge 1}{\mu }_n{v}_n(\mathbf{r})e^{-{\lambda }_{n}t},$$

$${\mu }_n=\frac{1}{\sqrt{V}}{\int }_{\Omega }{v}_n(\mathbf{r})\mathrm{d}\mathbf{r},$$

$$-(D_0{\nabla }^2-i\gamma gx)v_n={\lambda }_n{v}_n.$$

$$M(\mathbf{r},\delta +0)$$

$$=\frac{1}{\sqrt{V}}\sum _{n\ge 1}{\mu }_n^{\ast }{e}^{-{\lambda }_n^{\ast }\delta }\sum _{m\ge 1}{\beta }_{n,m}{v}_m(\mathbf{r}),$$

$${\beta }_{n,m}={\int }_{\Omega }{v}_n^{\ast }(\mathbf{r})v_m(\mathbf{r})\mathrm{d}\mathbf{r}.$$

$$M(\mathbf{r},2\delta )=\frac{1}{\sqrt{V}}\sum _{n,m\ge 1}{\mu }_n^{\ast }{\beta }_{n,m}{v}_m(\mathbf{r})e^{-({\lambda }_n^{\ast }+{\lambda }_m)\delta },$$

$$E={\int }_{\Omega }M(\mathbf{r},2\delta )\mathrm{d}\mathbf{r}=\sum _{n,m\ge 1}{\mu }_n^{\ast }{\beta }_{n,m}{\mu }_m{e}^{-({\lambda }_n^{\ast }+{\lambda }_m)\delta }.$$

$$M(\mathbf{r},2\delta )$$

$$E$$

$$2\mathrm{Re}({\lambda }_1)=|a_1^{\prime }|\frac{D_0}{{\ell }_g^2}+O({\ell }_g^{-3/2}),(g\to \infty )$$

$$-\ln (E)\propto ({\ell }_{\delta }/{\ell }_g{)}^2\propto (b{D}_0{)}^{1/3},(b{D}_0\gg 1)$$

$$2\mathrm{Re}({\lambda }_{k,l})$$

$$-\frac{\sqrt{3}{D}_0}{2|a_k^{\prime }|R_c{\ell }_g}+O({\ell }_g^{-1/2}),$$

$$2\mathrm{Re}({\lambda }_{k,l_1,l_2})$$

$$-\frac{\sqrt{3}{D}_0}{2|a_k^{\prime }|{\ell }_g}\left(\frac{1}{R_{c,1}}+\frac{1}{R_{c,2}}\right)+O({\ell }_g^{-1/2}),$$

$$M($$

$$+\left(c_{2,2}+c_{1,2}{e}^{2i\mathrm{Im}({\lambda }_1)\delta }\right)e^{-2\mathrm{Re}({\lambda }_1)\delta }{v}_2(x),$$

$$E\approx C{e}^{-2\mathrm{Re}({\lambda }_1)\delta },$$

$$C=2(C_{1,1}+\mathrm{Re}(C_{1,2})),$$

$$C_{1,2}={\mu }_1^{\ast }{\beta }_{1,2}{\mu }_2{e}^{2i\mathrm{Im}({\lambda }_1)\delta }.$$

$$2\mathrm{Im}({\lambda }_1)=\pm \left(\gamma gL-\sqrt{3}|a_1^{\prime }|\frac{D_0}{{\ell }_g^2}\right)+O({\ell }_g^{-3/2}),$$

$${\beta }_{n,m}={\int }_{\Omega }(\mathcal{D}{v}_n^{\ast })(\mathbf{r})v_m(\mathbf{r})\mathrm{d}\mathbf{r}.$$

$$E=\exp (-b{D}_0{\sin }^2\theta )E_{\mathrm{slab}}(L,g\cos \theta ),$$

$$D\approx D_0\left(1-\eta \frac{4}{3\sqrt{\pi }}\frac{S}{V}\sqrt{D_{0}T}\right),$$

$$E_{\mathrm{slab}}\approx 2C_{1,1}\exp (-|a_1^{\prime }|{\ell }_{\delta }^2/{\ell }_g^2),$$

$$E_{\mathrm{cyl}}\approx C\exp \left(-|a_1^{\prime }|\frac{{\ell }_{\delta }^2}{{\ell }_g^2}-\frac{{\ell }_{\delta }^2}{R^{1/2}{\ell }_g^{3/2}}+\frac{\sqrt{3}{\ell }_{\delta }^2}{2|a_1^{\prime }|R{\ell }_g}\right),$$

$$E_{\mathrm{rods}}\approx C\exp \left(-|a_1^{\prime }|\frac{{\ell }_{\delta }^2}{{\ell }_g^2}-\frac{{\ell }_{\delta }^2}{R^{1/2}{\ell }_g^{3/2}}-\frac{\sqrt{3}{\ell }_{\delta }^2}{2|a_1^{\prime }|R{\ell }_g}\right),$$

Full text

††journal: Journal of Magnetic Resonance

1 Introduction

Diffusion magnetic resonance imaging (dMRI) is a non-invasive technique which aims at unraveling the microstructural properties of a sample through the measurement of diffusion of spin-bearing nuclei (e.g., protons in water molecules) [1, 2, 3]. One may probe finer-scale details by increasing the magnetic field gradient $g$ or the $b$-value ($b\propto g^{2}$, see below), which characterizes the strength of the diffusion encoding. For example, Wedeen et al employed very large $b$-values, up to $40~{}\rm{ms/\mu m^{2}}$, for the reconstruction of neuronal fiber pathways in the human brain [4]. In contrast, most theoretical approaches rely on small $b$ approximations. In fact, the assumption of small $b$-values allows one to expand the dMRI signal in powers of $b$ (cumulant expansion) and truncate it to its first terms [3, 5, 6]. In particular, the first term of the cumulant expansion yields the classical Gaussian phase approximation for the signal $E=e^{-bD}$, where $D$ is the effective (or apparent) diffusion coefficient of spin-bearing particles in the sample that may differ from their intrinsic diffusion coefficient $D_{0}$. Note that the formula $E=e^{-bD_{0}}$ is exactly valid for free diffusion for any gradient profile and strength.

The validity of the cumulant expansion and related approximations was first questioned theoretically by Stoller et al. [7], who showed that the cumulant expansion has a finite convergence radius and that a new regime emerges at higher gradients. This so-called “localization regime” [8] was observed experimentally by Hürlimann et al. for water in slabs of width $160~{}\rm{\mu m}$ at gradients as small as $20~{}\rm{mT/m}$ and $b$-values of about $10~{}\rm{ms/\mu m^{2}}$ [9]. It features a very strong deviation from the Gaussian phase approximation signal, with $-\ln(E)\propto(bD_{0})^{1/3}$. Moreover, due to the high $b$-values the transverse magnetization far from the boundaries of the medium vanishes and only a thin layer of thickness $\ell_{g}=(\gamma g/D_{0})^{-1/3}$ near the boundaries contributes to the measured signal (with $\gamma$ being the gyromagnetic ratio of the diffusing nuclei). Here, $\ell_{g}$ denotes the gradient length, which can be interpreted as the typical diffusion length traveled by nuclei during the time needed to reach uncorrelated phases. This feature makes the signal at high $b$-values particularly sensitive to the microstructure [10]. The influence of the permeability of the boundaries was later studied in [11, 12] and it was shown that the localization of the magnetization near the boundaries enhances the dependence of the resulting signal on the permeability. The localization regime was shown to emerge when the gradient length $\ell_{g}$ is much smaller than (i) the diffusion length traveled by spins during the gradient sequence and (ii) the typical size of confining domains (e.g., width of the slab, in the experiment by Hürlimann et al).

To our knowledge, the work by Hürlimann et al [9] on diffusion in a slab is the only undeniable experimental evidence of the localization regime. While deviations from the Gaussian approximation were abundantly observed in biological tissues and mineral samples [13, 14, 15, 16, 17, 5, 18], it is usually difficult to identify unambiguously the origin of these deviations. In other words, the observed deviations may be related to the localization regime, but may also originate from distinct populations of water with different effective diffusion coefficients, mixture of restricted and hindered diffusion, etc. The existence of the localization regime in more general domains and, in particular, in unbounded domains (e.g., extracellular diffusion) has not yet been addressed experimentally. In this paper, we present experimental data, numerical simulations and theoretical arguments showing the emergence of the localization regime in two complementary geometries: diffusion inside cylinders and diffusion outside an array of rods. We also treat a slab geometry as a reference case.

The paper is organized as follows. First, we present in Sec. 2 the theoretical ground that we rely on in order to interpret the experimental results. Section 3 describes the experimental setup and the numerical simulations. In Sec. 4, we show that our theoretical analysis and numerical computations are in excellent agreement with the experimental data. The characteristic stretched-exponential decay of the signal in the localization regime is observed at moderately high gradients and in various geometries, including unbounded diffusion with obstacles. We discuss the implications of these results and conclude the paper in Sec. 5.

2 Theory

The evolution of the transverse magnetization $M$ in the simplest case of a constant magnetic field gradient $g$ is described by the Bloch-Torrey equation [19]:

[TABLE]

where we chose the $x$-axis as the gradient direction. The differential operator on the right-hand side of the equation, $D_{0}\nabla^{2}-i\gamma gx$, that we call the “Bloch-Torrey operator” in the following, is usually seen as a perturbation of the classical diffusion (or Laplace) operator $D_{0}\nabla^{2}$. However, the $i\gamma gx$ term makes the whole operator non-Hermitian, which affects deeply its spectral properties [7, 12, 20, 21]: its spectrum may become discrete or even empty in unbounded domains; the eigenvalues are complex; the eigenmodes are no longer orthogonal to each other (in the Hermitian sense) and they do not necessarily form a complete basis [22, 23]. These mathematical facts suggest that theoretical approaches where the gradient term is treated as a small perturbation may be insufficient to fully understand the signal formation.

In the following, we analyze the behavior of the transverse magnetization, hence the signal, at high gradients and relatively long encoding times, which are the conditions of emergence of the localization regime. This section summarizes and further extends the theoretical ground developed in [3, 7, 8, 9, 11, 10, 24]. Throughout this paper, we focus on the Bloch-Torrey equation (1) with the linear gradient, although similar features may be relevant to other spatial profiles of the magnetic field, e.g. the dipole field (see [25, 26] and references therein).

2.1 Eigenmode decomposition

Let us first consider a pulsed-gradient spin-echo sequence with gradient strength $g$, encoding pulse duration $\delta$, and no diffusion step between the two pulses (in other words, the time separation $\Delta$ between two encoding pulses is equal to $\delta$). We first assume that the domain $\Omega$ in which diffusion takes place is bounded.

In bounded domains, the term $i\gamma gx$ is bounded and one can show that the complex-valued eigenmodes $v_{n}$ of the Bloch-Torrey operator form a complete basis, except at some exceptional values of $g$ (that we neglect in the following, see [20]). We stress that these eigenmodes and the corresponding eigenvalues depend on the gradient $g$. Since the Bloch-Torrey operator is symmetric, the eigenmodes are “rectanormal” [7] in the sense that, after appropriate normalization, they satisfy the following orthogonality relation (not to be confused with the Hermitian orthogonality relation, in which $v_{m}({\mathbf{r}})$ is replaced by its complex conjugate $v_{m}^{*}({\mathbf{r}})$):

[TABLE]

This allows one to project functions onto the {$v_{n}\}_{n\geq 1}$ basis by using standard formulas. Let us assume that the initial condition is a uniform magnetization after the exciting $90^{\circ}$ radio frequency (RF) pulse, i.e. $M({\mathbf{r}},t=0)=1/V$, where $V$ is the volume of the domain $\Omega$. Hence we decompose the transverse magnetization as:

[TABLE]

where $\mu_{n}$ are complex coefficients given by

[TABLE]

and $\lambda_{n}$ are the corresponding eigenvalues:

[TABLE]

We set the minus sign in front of the Bloch-Torrey operator to have eigenvalues with a positive real part, thus we sort $\lambda_{n}$ by increasing real part.

At the end of the first gradient pulse (i.e., at $t=\delta$), a refocusing $180^{\circ}$ RF pulse is applied, which is equivalent to applying a complex conjugation to the magnetization $M$:

[TABLE]

where we have decomposed each $v_{n}^{*}$ over the eigenmode basis $(v_{m})_{m\geq 1}$, with complex coefficients $\beta_{n,m}$:

[TABLE]

Finally, after the second gradient pulse, the transverse magnetization is equal to

[TABLE]

and the resulting signal can be computed from

[TABLE]

We stress again that the eigenvalues $\lambda_{n}$ and eigenmodes $v_{n}({\mathbf{r}})$ depend on the gradient $g$, and so do the coefficients $\mu_{n}$ and $\beta_{n,m}$.

2.2 Localization regime and gradient length

In the regime of long encoding times $\delta$, one can truncate the above decompositions to their first terms:

[TABLE]

where the coefficient $C_{1,1}=\lvert\mu_{1}\rvert^{2}\beta_{1,1}$ typically decreases slowly with $g$ compared to the exponential factor. This approximation is valid if the next-order terms are negligible. This may happen if the next eigenvalue is sufficiently far from $\lambda_{1}$, i.e., $\lvert e^{-\lambda_{2}\delta}\rvert\ll\lvert e^{-\lambda_{1}\delta}\rvert$, or if the weights $\mu_{n}^{*}\beta_{n,m}$ and $\mu_{n}^{*}\beta_{n,m}\mu_{m}$ (with $(n,m)\neq(1,1)$) are small enough. Note that the eigenvalues may come in complex conjugated pairs, in which case the above expressions (10) and (11) would involve four terms (two diagonal terms and two cross terms) with similar exponential decays $e^{-2\mathrm{Re}(\lambda_{1})\delta}$. We show in Sec. 2.3 that these additional terms may produce oscillations in the signal on top of the overall decay (11). For clarity, we first focus on Eq. (11) in this subsection.

In this regime, the behavior of the magnetization, and thus the signal, is then dictated by the dependence of $\mathrm{Re}(\lambda_{1})$ on the gradient $g$. In [7, 8, 11, 12, 20], an asymptotic expansion for the $\lambda_{n}$ at high gradients was derived. In particular, the first term of the expansion of the first eigenvalue is universal:

[TABLE]

where $\ell_{g}=(\gamma g/D_{0})^{-1/3}$ is the gradient length and $a^{\prime}_{1}\approx-1.02$ is the first zero of the derivative of the Airy function $\mathrm{Ai}(z)$. This formula is universal in the sense that it does not depend on the geometry of the domain. By combining Eqs. (10) and (12), one gets the stretched-exponential behavior, which is a hallmark of the localization regime,

[TABLE]

where $\ell_{\delta}=\sqrt{D_{0}\delta}$ is the diffusion length during the gradient pulse, and $b=2\gamma^{2}g^{2}\delta^{3}/3$ for this gradient profile.

In order to offer an interpretation for the gradient length $\ell_{g}$, let us consider several nuclei that start at the same position and diffuse freely (i.e., they do not encounter any boundary or obstacle) under the magnetic field gradient $g$. After a time $t$, the nuclei are typically spread over a distance of the order of $\sqrt{D_{0}t}$, and thus have accumulated a random phase difference of the order of $(\gamma gt)\sqrt{D_{0}t}$. The nuclei have practically uncorrelated phases when this difference is comparable with $2\pi$, which yields a decorrelation time: $t_{g}=(\gamma^{2}g^{2}D_{0})^{-1/3}$. Note that by time-reversal, $t_{g}$ can also be interpreted as the time beyond which several diffusing particles arriving at the same position have almost uncorrelated phases. The gradient length $\ell_{g}$ is then the typical distance traveled by these particles during time $t_{g}$: $\ell_{g}=\sqrt{D_{0}t_{g}}$.

If the gradient length $\ell_{g}$ is much smaller than the typical size of the medium $\ell_{s}$ and the diffusion length $\ell_{\delta}$, then nuclei at any point not too close to a boundary have completely uncorrelated phases and the resulting magnetization is negligible. As the motion of the nuclei is restricted near the boundary, the remaining magnetization at the echo time is localized in a thin layer of thickness $\ell_{g}$ close to the boundary: this is the localization regime. In regard to Eq. (10), the localization regime corresponds to the localization of the first eigenmodes of the Bloch-Torrey operator at high gradients over a width $\ell_{g}$ [11, 10, 20]. This localization occurs where the gradient is orthogonal to the boundary, i.e. where the boundary prevents diffusion along the gradient direction the most.

2.3 Dependence on the geometry

Although the first term of the high-gradient expansion (12) of the eigenvalues of the Bloch-Torrey operator is universal, the next-order terms depend on the geometry of the domain, and more precisely on the curvature, permeability, and surface relaxivity of the boundary at the points where the localization occurs [11, 12, 20]. Throughout the article, we discard permeability and surface relaxivity and consider only the effect of curvature. In two dimensions, if one denotes by $R_{c}$ the curvature radius at a localization point, then the eigenvalues associated to the eigenmodes localized near that point can be expanded as [20]

[TABLE]

where $a^{\prime}_{k}$ is the $k$-th zero of the derivative of the Airy function, and $k,l$ are positive integers (the smallest eigenvalue is given by $k=l=1$). This formula applies to diffusion inside cylinders or outside parallel rods, which can be effectively reduced to diffusion in a two-dimensional geometry. A generalization of Eq. (14) to three-dimensional domains is expected to be

[TABLE]

where $R_{c,1}$ and $R_{c,2}$ are the principal radii of curvature of the boundary at the localization point and $k,l_{1},l_{2}$ are positive integers. In the following we focus on the formula (14) for two-dimensional geometries, which corresponds to our experimental setup (see Sec. 3).

In the third term of Eq. (14) the sign of the curvature radius $R_{c}$ should be taken into account: positive for a concave boundary (e.g., interior of a cylinder) and negative for a convex boundary (e.g., exterior of a cylinder). In the above expansion, the first term is related to the variation of the eigenmode along the gradient direction and the second term is related to the variation of the eigenmode in the orthogonal plane.

The real part of the spacing between the first and the next eigenvalue in Eq. (14) is $2D_{0}/(\lvert R_{c}\rvert^{1/2}\ell_{g}^{3/2})$ which means that although the localization regime emerges when $\ell_{\delta}\gg\ell_{g}$, the asymptotic decay (11) is established at a longer time: $\ell_{\delta}\gg\lvert R_{c}\rvert^{1/4}\ell_{g}^{3/4}$ (assuming that $\ell_{g}<\lvert R_{c}\rvert$ in the localization regime, which is always the case for basic geometries). In other words, in the intermediate regime $\ell_{g}\ll\ell_{\delta}\ll\lvert R_{c}\rvert^{1/4}\ell_{g}^{3/4}$, the transverse magnetization is localized at the boundary of the domain, but the truncated formulas (10) and (11) are not valid since one has to take into account other eigenmodes that are not suppressed by the exponential factor in Eq. (8). These eigenmodes have typically the same profile in the gradient direction but different profiles in the orthogonal plane. Since these eigenmodes correspond to different eigenvalues, the signal is a superposition of different decaying functions (as in Eq. (9)) and cannot be reduced to a single stretched-exponential decay.

Another effect of the geometry is the overlapping of localization pockets. In fact, to each localization point (at which the gradient is perpendicular to the boundary) is associated a family of eigenmodes and eigenvalues described by Eq. (14) at high $g$. As a consequence, the smallest eigenvalues from different families can contribute to the signal even at large $\delta$, whereas the overlap between corresponding eigenmodes may lead to oscillating patterns.

To illustrate this effect, let us consider the simple example of a slab of width $L$ (with the gradient direction being orthogonal to the slab). Note that this geometry is reduced to a one-dimensional interval of length $L$. In this setting, the eigenmodes of the Bloch-Torrey operator can be localized at either of two endpoints of the interval, but only if the gradient length $\ell_{g}$ is much smaller than $L$. In fact, as we stated in Sec. 2.2, $\ell_{g}$ represents approximately the width of the localized eigenmodes. For example, Stoller et al computed that the first eigenmodes start to localize at each side of the slab when $\ell_{g}<0.38L$. Due to the left-right symmetry ($x\to-x$), the first two eigenmodes $v_{1}(x)$ and $v_{2}(x)$ satisfy the identity $v_{2}(x)=v_{1}^{*}(-x)$ and $\lambda_{2}=\lambda_{1}^{*}$, and are localized at each endpoint of the interval. Thus the first two eigenvalues have the same real part and, in the limit of large $\delta$, the magnetization may be represented by a superposition of $v_{1}$ and $v_{2}$:

[TABLE]

with $c_{n,m}=\mu_{n}^{*}\beta_{n,m}/\sqrt{V}$, and the signal is given by

[TABLE]

with

[TABLE]

where $C_{1,1}$ was previously introduced in Eq. (11) and $C_{1,2}$ is an oscillating function of $g$ and $\delta$ given by

[TABLE]

The factor $2$ in (18) reflects the fact that two eigenmodes contribute to the signal. For a slab, one can compute the following expansion for the imaginary part of the first eigenvalue [7]:

[TABLE]

where the $\pm$ sign is merely a matter of convention (given that $\lambda_{2}=\lambda_{1}^{*}$). Note that the eigenmodes of the Bloch-Torrey operator in a slab are known explicitly and can be written in terms of Airy functions, see [7, 11, 12].

The conclusion of this computation is that the cross-term $C_{1,2}$ produces oscillations in the signal on top of the asymptotic decay (11). The definition (7) of $\beta_{n,m}$ implies that this cross-term is linked to the overlapping of the modes $v_{1}$ and $v_{2}$. In turn, this overlapping depends on the ratio between the width $\ell_{g}$ of the modes, and their spacing $L$: the smaller the ratio $\ell_{g}/L$, the smaller the overlapping, and thus, the smaller the oscillating term. In the limit of well-separated modes, $\beta_{1,2}=0$ and one has $C=2C_{1,1}$. The effect of overlapping of eigenmodes is illustrated in Fig. 1 where we plotted the transverse magnetization in a slab for two different gradient strengths: at low gradients the localization pockets overlap and at high gradients they are well-separated. Note that we plotted only the real part of the magnetization; the imaginary part is non-zero but it does not contribute to the signal since its integral over the slab is zero. Although we illustrated this overlapping effect on the simple example of a slab, the conclusion (and the previous computations) may be generalized to any geometry, where $L$ would denote the spacing between two localized eigenmodes.

2.4 The influence of $\Delta-\delta$

Until now we have considered a PGSE sequence with $\Delta=\delta$, i.e. no diffusion step between two gradient pulses. In this section we consider a more general PGSE sequence and investigate the influence of the diffusion step on the transverse magnetization. The main parameter of this study is the duration of the diffusion step, $\Delta-\delta$. Mathematically, the effect of this diffusion step is to multiply the transverse magnetization just after the first gradient pulse by the operator $\mathcal{D}=\exp((\Delta-\delta)D_{0}\nabla^{2})$. For simplicity we assume that the $180^{\circ}$ RF pulse occurs immediately after the first gradient pulse and not in the middle of the diffusion step (note that this assumption does not affect the final result). Then Eqs. (8), (9), and the consequent analysis remain applicable with the new definition for the coefficients $\beta_{n,m}$:

[TABLE]

The dependence on $\Delta-\delta$ is now hidden in the coefficients $\beta_{n,m}$. Note that if $\Delta=\delta$, then the operator $\mathcal{D}$ is the identity (since diffusion during zero time does not affect the magnetization) and we recover the previous definition (7) for $\beta_{n,m}$. If we consider the regime of long $\delta$ where only the first eigenmode contributes, Eq. (8) is reduced to Eq. (10) and thus the diffusion step results merely in a modification of the global amplitude of the magnetization through the new definition (21) of the coefficient $\beta_{1,1}$. In turn, the global amplitude of the signal in Eq. (11) is also affected by this coefficient and decreases with increasing $\Delta-\delta$.

The situation of two overlapping eigenmodes is more complex. The diffusion step also changes the values of the coefficients $\beta_{1,2}$, $\beta_{2,1}$, and hence the relative amplitudes of the oscillating term $C_{1,2}$ (see Eqs. (18) and (19)). Diffusion is expected to increase the width of the localized eigenmodes, thus increasing the overlapping between them and enhancing the oscillations in the signal. This effect should be stronger for already overlapping modes than for well separated ones, in other words, when the ratio $\ell_{g}/L$ is not very small, where $L$ denotes the distance between two localization pockets. The effect of the diffusion step on the localization regime in a one-dimensional setting was partly investigated in [11].

2.5 Unbounded domains

In the previous subsections, we assumed that diffusion takes place in a bounded domain $\Omega$, and this hypothesis ensured the validity of the decomposition of the transverse magnetization over the Bloch-Torrey eigenmodes basis (3). However, this assumption excludes some physically relevant cases, such as extracellular diffusion, which are described by unbounded domains. In this section, we discuss few known mathematical properties and the applicability of the previous results to such domains.

The most striking property of the Bloch-Torrey operator is that, in most unbounded domains, its spectrum is discrete [21], even though the spectrum of the Laplacian operator $\nabla^{2}$ is continuous. In other words, the $i\gamma gx$ term is responsible for the discreteness of the spectrum. However, there is no rigorous mathematical proof that the eigenmodes form a complete basis in such domains. We conjecture that this property still holds except for trivial cases (such as free space, which leads to an empty spectrum). In [20], a construction of approximate eigenmodes and eigenvalues was performed and Eqs. (14) and (20) are valid for unbounded domains as well.

Despite its practical importance, few theoretical studies addressed the validity of the Gaussian phase approximation in unbounded domains [10]. In the following sections, we investigate experimentally the emergence of the localization regime in such a domain (array of rods).

3 Material and Methods

Experiments were performed using hyperpolarized xenon-129 gas ($\gamma\approx 74\cdot 10^{6}~{}\rm{s^{-1}T^{-1}}$) continuously flowing through phantoms containing different diffusion barrier geometries. Utilizing gas diffusion compared to water diffusion entails a several orders of magnitude larger diffusion coefficient which allows probing structures on the millimeter scale, which can easily be constructed with 3D-printers. Due to the weak signal of thermally polarized xenon gas, hyperpolarized gas with a considerably higher NMR signal was employed. Hyperpolarization was achieved by Rb/Xe-129 spin-exchange optical pumping (SEOP) [27, 28, 29]. For technical reasons [27], a gas mixture (Air Liquide Deutschland GmbH, Düsseldorf, Germany) composed of xenon (0.95 Vol %) and nitrogen (8.75 Vol%) and helium-4 (rest) was used. The free diffusion coefficient of xenon in this gas mixture was measured to be $D_{0}=(3.7\pm 0.2)\cdot 10^{-5}~{}\rm{m^{2}s^{-1}}$ [30], which is one order of magnitude larger than for pure xenon gas [31]. The gas was transferred at a small constant flow (approximately $150~{}\rm{mL/min}$) to the phantom positioned in an in-house built xenon coil in the isocenter of the magnet of a $1.5~{}\rm{T}$ clinical MR scanner (Magnetom Symphony, A Tim System, Siemens Healthcare, Erlangen, Germany) with a maximal employed gradient amplitude of $32~{}\rm{mT/m}$. The experimental set-up and hyperpolarization process are detailed in [30, 32].

The phantoms used are illustrated in Fig. 2 and described in Table 1. Phantom $1$ contains parallel plates separated by a distance of $L=3~{}\rm{mm}$, built by the in-house workshop from PMMA. For the phantoms $2$ and $3$, two blocks containing cylindrical tubes (with two different diameters $2R$) in a hexagonal arrangement were 3D-printed. Here, the gas diffuses inside the cylinders. Since all cylinders are identical and isolated from each other, this setting is equivalent to a single cylinder of diameter $2R$. Phantoms $4$ and $5$ consist of cylindrical solid rods on a square grid attached to a base plate and a roof plate with holes for gas in and out flow, so that the gas diffuses outside the cylinders. The geometry is defined by the diameter of the rods ($2R$) and by the rod center-to-center distance ($L$).

Phantoms $2$-$5$ were printed with the PolyJet technology (Objet30 Pro, VeroClear as printing material, Stratasys, Ltd., Eden Prairie, MN, USA) and then inserted in a casing with a base area of $70.1~{}\rm{mm}\times 42.3~{}\rm{mm}$. Consequently, for phantoms $4$ and $5$, there are border regions between the array of rods and the enclosing walls in which diffusion also takes place. For all phantoms, surface relaxivity and permeation can be ignored for high-gradient experiments.

Phantom $1$ was positioned with the gas flow directed in the horizontal direction perpendicular to the main magnetic field and the gradient vector was pointing in the vertical direction in a sagittal slice of $50~{}\rm{mm}$ thickness orthogonal to the plates and to the gas flow direction, see Fig. 2. Phantoms $2$-$5$ were positioned with the gas flow directed parallel to the main magnetic field and the gradients were applied in the transversal plane of the scanner orthogonal to the gas flow direction (slice thickness $45~{}\rm{mm}$). For phantoms $4$ and $5$, measurements were taken with the gradient vector pointing along the left-right direction or the diagonal direction. In all cases, the gas flow did not influence the diffusion-weighted NMR signal.

A pulsed-gradient spin-echo (PGSE) sequence was applied as depicted in Fig. 3. The durations $\delta$ of the trapezoidal gradient pulses were set to $6~{}\rm{ms}$ and include flat top time plus the ramp-up time of $\epsilon=0.32~{}\rm{ms}$. The gradient separation time was $\Delta=9.34~{}\rm{ms}$. The diffusion-weighted signal was sampled by gradually increasing the gradient amplitude in 32 steps from 0 to $g_{\rm{max}}=32~{}\rm{mT/m}$, recording the spin echo signal and acquiring up to $15$ averages. The time between two consecutive $90^{\circ}$ excitation pulses, i.e. the time between two measurements, was set to $18~{}\rm{s}$ to restore the polarization in the phantom via gas flow. To account for fluctuations in the polarization level the recorded spin echo was averaged and normalized to an additional signal pre-readout directly after the $90^{\circ}$ excitation pulse. To obtain the diffusion-induced signal attenuation, all points were normalized to the point acquired without diffusion weighting (i.e., at $g=0$).

The numerical computation of the signal in a slab and in a cylinder is effectively reduced to that in an interval and a disk, respectively. For these simple shapes, the most efficient and accurate computation of the signal is realized with the matrix formalism [33, 34, 35, 3], in which the Bloch-Torrey equation is projected onto the basis of explicitly known Laplacian eigenfunctions to represent the signal via matrix products and exponentials (see [36, 37] for details). The matrix formalism was also used to compute the transverse magnetization in these two domains (see similar computations in [11]). In turn, the numerical computation of the signals for phantoms $4$ and $5$ (arrays of rods) was performed differently. While the matrix formalism could in principle be applied, the need for a numerical computation of Laplacian eigenfunctions in such structures makes this approach less efficient. Thus, we performed Monte Carlo simulations including the borders around the rod arrays with $2.5\cdot 10^{7}$ random walkers and $1.6\cdot 10^{5}$ steps per random walk trajectory. In order to compute the eigenmodes of the Bloch-Torrey operator in the rods geometry, we used the PDE solver from Matlab (The MathWorks, Natick, MA USA) in a square array of $3\times 3$ rods with Dirichlet boundary conditions on the outer boundary and we kept only the eigenmodes that were localized on the central rod. Since the distance from the central rod to the outer boundary is about $L$, which was much larger than $\ell_{g}$ and $\ell_{\delta}$ in our simulations, the effect of the outer boundary is negligible, so that the computed eigenmodes are very close to the ones for the infinitely periodic array of rods.

4 Results

We present the experimental and numerical results for diffusion in three different geometries: inside slabs, inside cylinders, and outside arrays of rods. The localization regime in the slab geometry was already investigated experimentally by Hürlimann et al [9], whereas only a few theoretical studies were devoted to the cylinder geometry [8, 10, 20]. The signal for an array of rods was previously computed numerically in [10] using a finite-element method [38].

For our particular gradient sequence and the parameters of the xenon gas mixture, one can compute $\ell_{\delta}=0.5~{}\rm{mm}$ and $\ell_{g}$ decreases from $0.8~{}\rm{mm}$ to $0.25~{}\rm{mm}$ for $g$ ranging from $1~{}\rm{mT/m}$ to $32~{}\rm{mT/m}$. So, by increasing the gradient, $\ell_{g}$ crosses $\ell_{\delta}$ and the localization regime emerges.

4.1 Slab geometry

We choose the axes such that the slab is orthogonal to the $z$-axis. Note that this convention is different from the one that we adopted in Sec. 2, where the gradient was directed along the x-axis. One can then decompose the diffusive motion independently along the three axes $x,y,z$ and get

[TABLE]

where $\theta$ is the angle between the gradient and the $z$-axis and $E_{\rm slab}(L,g)$ is the signal from an interval of length $L$. Here, the first factor is the signal attenuation due to diffusion in the lateral plane $xy$, which is almost free as outer boundaries are separated by distances that greatly exceed both $\ell_{g}$ and $\ell_{\delta}$ (about $40~{}\rm{mm}$). One gets the slowest decay by ensuring that the gradient is orthogonal to the slab, i.e. $\theta=0$, which was chosen in the experiments and numerical simulations. This is expected since in this situation the boundary restricts diffusion along the gradient direction the most.

The signals are presented in Fig. 4. At low gradients, the Gaussian phase approximation is valid with an effective diffusion coefficient $D$. Note, however, the deviation from the free diffusion signal $E=e^{-bD_{0}}$ due to restriction by the slab. A short-time analysis [39, 40, 8, 41] yields an apparent diffusion coefficient

[TABLE]

where $S/V$ is the surface-to-volume ratio of the confining domain ($S/V=2/L$ for a slab), $\eta\approx 0.9$ is a numerical prefactor that depends on the sequence, and $T=\Delta+\delta$ is the duration of the gradient sequence.

This is a short-time approximation in the sense that $\eta\sqrt{D_{0}T}/L$ should be small enough. In addition, this formula relies on the Gaussian phase approximation which requires small $b$-values. For our parameters, we get $D\approx 0.66D_{0}$ and the agreement between $e^{-bD}$ and the signal at low gradients is good.

At higher gradients, the signal deviates significantly from the Gaussian regime and follows the asymptotic decay (see Eq. (17) and Refs. [3, 7, 9, 8, 10])

[TABLE]

where the prefactor $C_{1,1}$ can be computed exactly [9, 11] and scales as $\ell_{g}/L$. One can interpret this prefactor as the fraction of spins inside the two layers of thickness $\sim\ell_{g}$ where the signal is localized. As we discussed in Sec. 2.4, $C_{1,1}$ decreases with increasing diffusion step duration $\Delta-\delta$. Note that the localization regime emerges at gradients as small as $10~{}\rm{mT/m}$. Since the width of the slab is much greater than $\ell_{g}$, one can treat the localization layers on both sides of the slab as independent from each other (see Fig. 1 and the discussion in Sec. 2.3).

We observe a remarkable agreement between the experimental data, exact solution via the matrix formalism, and the asymptotic relation (24). Note that the latter contains no fitting parameter (the prefactor $C_{1,1}$ was found by computing the eigenmodes numerically). Systematic minor deviations of the experimental points may be caused by weak misalignment of the gradient direction (i.e. $\theta$ slightly different from [math] in Eq. (22)) or weak surface relaxivity. Note that we performed all computations with a square gradient profile instead of a trapezoidal one (in other words, with $\epsilon=0$, see Fig. 3) and we checked that not accounting for the trapezoidal profile had a negligible influence on the computed signal due to the very short ramp-up time ($\epsilon=0.34~{}\rm{ms}$).

4.2 Diffusion inside a cylinder

If the diameter of the cylinder is much larger than the gradient length $\ell_{g}$ and the diffusion length $\ell_{\delta}$, then the cylinder geometry is locally similar to the slab geometry except for the smaller “localization pocket” where the magnetization does not vanish. An argument similar to that of Eq. (22) implies that the magnetization is localized where the angle between the gradient and the normal vector to the boundary is close to [math] (this condition becomes more and more restrictive at higher $b$-values according to Eqs. (22) and (24)). The transverse magnetization inside a cylinder, obtained with a matrix formalism computation, supports this argument (Fig. 5). In fact, as the gradient increases, the magnetization gradually transforms from a flat uniform profile to the one that is localized around two opposite points on the boundary of the cylinder and displays two independent pockets at sufficiently high gradients.

In Fig. 6 we show the signal for diffusion inside a cylinder of diameter $3.8~{}\rm{mm}$. Although the signal decays faster than in the slab, one observes a similar stretched-exponential behavior. Using Eqs. (14) and (17), one gets the asymptotic decay for the cylinder:

[TABLE]

where $C$ is given by Eq. (18). Here, $R$ is large enough so that there is no overlapping between the first two eigenmodes, and $C=2C_{1,1}$. The prefactor $C_{1,1}$ can be computed numerically from the eigenmodes and scales approximately as $(\ell_{g}/R)^{3/2}$. One observes the perfect agreement between experiment, matrix formalism computation, and asymptotic formula at high gradients (without any fitting parameter).Note also that Eq. (24) with only the leading term is not accurate here so that the correction terms in the exponential are indeed important.

For a cylinder of a smaller diameter ($2R=2~{}\rm{mm}$, see Fig. 7), the signal shows some oscillations that are reminiscent of diffusion-diffraction patterns [13, 42, 43, 44, 34]. This is the case discussed in Sec. 2.3 where two localization pockets overlap because $\ell_{g}/(2R)$ is not small enough. The signal is still given by Eq. (25) but one cannot neglect the cross-term $C_{1,2}$ in the expression (18) of $C$. The oscillations are then described by $C_{1,2}$ and appear on top of the asymptotic stretched-exponential decay. These oscillations shown in Fig. 7 are very well reproduced by the asymptotic formulas (19) and (20) with $L=2R$, where the coefficients $\mu_{1}$, $\mu_{2}$, and $\beta_{1,2}$ were computed numerically from the eigenmodes. We stress that these coefficients generally depend on $g$ and that $\beta_{1,2}$ additionally depends on the diffusion step duration $\Delta-\delta$. This overlapping phenomenon is supported by Fig. 8 which illustrates that the magnetization inside the cylinder is not well localized even at the highest gradient available.

4.3 Diffusion outside an array of rods

The geometry of phantoms $4$ and $5$ is defined by $L$, the center-to-center spacing between rods, and $2R$, the diameter of the rods. We consider three different cases: (a) phantom $4$ ($L=4~{}\rm{mm}$ and $2R=3.2~{}\rm{mm}$) with the gradient vector in the diagonal direction; (b) phantom $5$ ($L=3.4~{}\rm{mm}$ and $2R=3.2~{}\rm{mm}$) and gradient vector in the diagonal direction; (c) phantom $5$ and gradient vector in the horizontal direction.

The main difference between these three cases is the spacing $e_{\rm pock}$ between two neighboring rods along the gradient direction, i.e. the spacing between two neighboring localization pockets (see Fig. 9): $e_{\rm pock}=e_{\rm d}=\sqrt{2}L-2R=2.5~{}\rm{mm}$ in (a); $e_{\rm pock}=e_{\rm d}=1.6~{}\rm{mm}$ in (b); and $e_{\rm pock}=e_{\rm h}=0.2~{}\rm{mm}$ in (c). Figure 10 shows the signal for these three cases, ordered by descending $e_{\rm pock}$. Note that here the signal is formed by the magnetization localized near the rods and by the magnetization localized near the borders of the casing in which the phantom is enclosed. We did not plot the low-$b$, short time approximation $E=e^{-bD}$ here because the surface-to-volume ratio of the structure is too large so that the approximate formula (23) for $D$ is not valid.

First of all, one can note an excellent agreement between experimental data and Monte Carlo simulations. The high-$g$ asymptotic behavior of the eigenvalues and the signal in a rods geometry is similar to the one for cylinders in Eq. (25) except for a sign change due to the opposite curvature [20]:

[TABLE]

where $C$ is given by Eq. (18) and may be computed numerically from the eigenmodes. This formula matches very well the signal at high gradients in case (a) (see Fig. 10 (a)).

In the previous subsection, we already saw the signal without oscillations (Fig. 6) due to well-localized states, as well as the signal with oscillations (Fig. 7) due to a partial overlap of two localization pockets when the size of the confining domain is not very large compared to $\ell_{\delta}$ and $\ell_{g}$. Here the same phenomenon occurs. With $e_{\rm pock}$ much larger than $\ell_{g}$ (case (a)), there is little overlapping between the localization pockets of neighboring rods. This ensures the localization of the eigenmodes and the small amplitude of the oscillations in the signal (see Sec. 2.3), i.e. the validity of $C=2C_{1,1}$. In the case where the ratio between $e_{\rm pock}$ and $\ell_{g}$ is smaller, i.e. localization pockets overlap more (case (b)), more pronounced oscillations on top of the overall decay (26) arise (see Fig. 10 (b)). The signal is still described by Eq. (26), and the oscillations are contained in the cross-term $C_{1,2}$ from the expression (18) of $C$ and may be computed from Eqs. (19) and (20) (with $L=e_{\rm pock}$). Systematic deviations between the exact signal and the asymptotic formulas may be attributed to the truncation of Eq. (14), neglecting higher-order modes in the expression of the signal, and not accounting for the borders in the experimental setup.

In turn, in case (c) $e_{\rm pock}$ is smaller than $\ell_{g}$ even at the highest gradient available. If one could approximate the small space between two neighboring rods as a slab, then the ratio $\ell_{g}/e_{\rm pock}$ would be too large for localization to be relevant (see Sec. 2.3). Our conjecture is that the residual signal at high gradients may be interpreted as a kind of motional narrowing inside the small gaps of width $e_{\rm pock}$ between the rods (see [45, 3]).

5 Discussion and Conclusion

We have observed and described the localization regime in three geometries: slab, cylinder, and array of circular obstacles (rods). The localization regime appears whenever the gradient length $\ell_{g}=(\gamma g/D_{0})^{-1/3}$ is much smaller than the diffusion length $\ell_{\delta}=(D_{0}\delta)^{1/2}$ and any relevant geometrical length scale $\ell_{s}$ of the medium along the gradient direction. Thus, it is universal at high gradients and non-narrow pulses. In this regime, the transverse magnetization is localized near the obstacles, boundaries, or membranes of the sample. For this reason, the signal is particularly sensitive to the microstructure of the medium. In particular, possible oscillations of the signal are caused by a partial overlap between localized magnetization pockets and thus contain information about mutual arrangements of obstacles.

Let us clarify the role of the conditions (i) $\ell_{g}\ll\ell_{\delta}$ and (ii) $\ell_{g}\ll\ell_{s}$. Condition (i) ensures that the eigenmode decomposition of the transverse magnetization may be truncated to its first terms. In turn, the signal decays exponentially with $\delta$ and its dependence on $g$ is essentially determined by the first eigenvalue of the Bloch-Torrey operator. Furthermore, condition (ii) is necessary for the localization of the eigenmodes of the Bloch-Torrey operator and the validity of the high-$g$ expansion of its eigenvalues. Thus, both conditions are required for the localization of the transverse magnetization and the stretched-exponential decay of the signal with $g$.

An extreme case where condition (i) is not satisfied would be narrow gradient pulses experiments ($g\to\infty$ and $\delta\to 0$). Although these experiments require high gradients, they do not achieve the localization regime. In fact, whereas the localization regime emerges when $(\gamma g)^{2}\delta^{3}D_{0}\gg 1$, narrow-pulse experiments correspond to $\gamma g\to\infty$ and $\delta\to 0$ such that $\gamma g\delta=q$ is a finite value. This leads to $(\gamma g)^{2}\delta^{3}D_{0}\to 0$. In other words, the signal attenuation is not produced by the encoding step but by the subsequent diffusion step with $g=0$. This is evident from the fact that the signal is left unchanged if one sets $\Delta=0$, i.e. no diffusion time between two short gradient pulses. On the other hand, condition (ii) is typically not satisfied in the motional narrowing regime ($g\to 0$ and $\delta\to\infty$). Then, the eigenmodes of the Bloch-Torrey operator are close to the eigenmodes of the Laplace operator. In particular, the first Laplacian eigenmode is constant if one assumes impermeable, non-relaxing boundaries. Therefore, the transverse magnetization at long times is uniform inside the sample.

Slab and cylinder are confined geometries that may model an intracellular space. It is well-known that such domains produce non-Gaussian signals, for example in the limit of narrow pulses (e.g. diffusion-diffraction patterns, see [13, 42, 43, 44, 34]). However, the signal from extracellular space is almost always assumed to be Gaussian, and non-Gaussian effects are attributed to multiple contributing pools. In other words, one assumes that the obstacles may be treated as an effective medium with an effective diffusivity $D$. The measurements performed on the array of rods show that the extracellular signal is not Gaussian at high gradients and that it simply results from the localization of the signal at the outer boundaries of the obstacles. Ignoring this effect may lead to false interpretations created by commonly used fitting models.

We stress that the localization regime in our setting starts to emerge with moderate gradients of about $10~{}\rm{mT/m}$, or at $bD_{0}$ about $4$. These gradients are easily achieved in most clinical scanners. Note that the localization regime emerges under the condition $\ell_{\delta}/\ell_{g}\gg 1$ which can be restated as $bD_{0}\propto\gamma^{2}D_{0}g^{2}\delta^{3}\gg 1$. In order to rescale the experimental conditions from xenon gas to water, we compute the ratio $(\gamma^{2}D_{0})_{\rm xenon}/(\gamma^{2}D_{0})_{\rm water}\approx 10^{3}$. This means that in order to have the same value of $\gamma^{2}D_{0}g^{2}\delta^{3}$, one has to increase $g$ and $\delta$ such that $g^{2}\delta^{3}$ is $10^{3}$ times larger with water than with xenon. For example, the experiments by Hürlimann et al were performed with $\delta=60~{}\rm{ms}$ which is approximately $10$ times longer than in our experiments, and gradients of comparable magnitude as ours (around $20~{}\rm{mT/m}$).

Theoretical aspects of diffusion NMR have been broadly explored at weak gradients by means of perturbation theory. While surface-to-volume ratio, overall curvature, permeability and surface relaxivity could be estimated at short times [39, 40, 41], the structural organization of barriers is accessible at long times [46, 47]. In biomedical applications, the apparent diffusion coefficient and kurtosis are employed as biomarkers to detect stroke, tumors, lesions, partial tissue destruction, etc. (see [18] and references therein). In spite of its considerable progress in medicine and material sciences, the comprehensive theory of diffusion NMR remains to be elaborated. The recent mathematical advances show that the gradient term presents a singular perturbation to the Laplace operator in unbounded domains, resulting in the discrete spectrum of the Bloch-Torrey operator [21]. But even in bounded domains, the presence of the gradient term may lead to branching points of the spectrum and thus a finite radius of convergence of the cumulant expansion [7]. These mathematical evidences urge for revising the conventional perturbation theory paradigm and developing non-perturbative approaches to diffusion NMR. Even so the localization regime is yet poorly understood and exploiting its potential advantages is still challenging in experiments (partly due to strongly attenuated signals), the high sensitivity of the signal to the microstructure at high gradients is a promising avenue for creating new experimental protocols. If former theoretical efforts were essentially focused on eliminating the dephasing effects and reducing the mathematical problem to pure diffusion, future developments have to aim at exploiting the advantages of high gradients and accentuating the differences from pure diffusion.