Nicolas Moutal
email: nicolas.moutal@polytechnique.edu

Kerstin Demberg

Denis Grebenkov
email: denis.grebenkov@polytechnique.edu

Tristan Anselm Kuder

Localization regime in diffusion NMR: theory and experiments

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.

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) [START_REF] Callaghan | Principles of Nuclear Magnetic Resonance Microscopy[END_REF][START_REF] Price | NMR studies of translational motion: Principles and applications[END_REF][START_REF] Grebenkov | NMR survey of reflected Brownian motion[END_REF]. One may probe finerscale details by increasing the magnetic field gradient g or the b-value (b ∝ 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 ms/µm 2 , for the reconstruction of neuronal fiber pathways in the human brain [START_REF] Wedeen | The geometric structure of the brain fiber pathways[END_REF]. 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 [START_REF] Grebenkov | NMR survey of reflected Brownian motion[END_REF][START_REF] Kiselev | Fundamentals of diffusion MRI physics[END_REF][START_REF] Frøhlich | Effect of impermeable boundaries on diffusion-attenuated MR signal[END_REF]. 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. [START_REF] Stoller | Transverse spin relaxation in inhomogeneous magnetic fields[END_REF], 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" [START_REF] De Swiet | Decay of nuclear magnetization by bounded diffusion in a constant field gradient[END_REF] was observed experimentally by Hürlimann et al. for water in slabs of width 160 µm at gradients as small as 20 mT/m and b-values of about 10 ms/µm 2 [START_REF] Hürlimann | Spin echoes in a constant gradient and in the presence of simple restriction[END_REF]. It features a very strong deviation from the Gaussian phase approximation signal, with ln(E) ∝ (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 ℓ g = (γg/D 0 ) -1/3 near the boundaries contributes to the measured signal (with γ being the gyromagnetic ratio of the diffusing nuclei). Here, ℓ 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 [START_REF] Grebenkov | Diffusion MRI/NMR at high gradients: Challenges and perspectives[END_REF]. The influence of the permeability of the boundaries was later studied in [START_REF] Grebenkov | Exploring diffusion across permeable barriers at high gradients. II. Localization regime[END_REF][START_REF] Grebenkov | The complex Airy operator on the line with a semipermeable barrier[END_REF] 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 ℓ 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 [START_REF] Hürlimann | Spin echoes in a constant gradient and in the presence of simple restriction[END_REF] 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 [START_REF] Callaghan | Diffraction-like effects in NMR diffusion studies of fluids in porous solids[END_REF][START_REF] Niendorf | Biexponential diffusion attenuation in various states of brain tissue -implications for diffusionweighting imaging[END_REF][START_REF] Jensen | Diffusion kurtosis imaging: the quantification of non-Gaussian water diffusion by means of magnetic resonance imaging[END_REF][START_REF] Jacob | 3D MRI of non-Gaussian 3 He gas diffusion in the rat lung[END_REF][START_REF] Laun | Determination of the defining boundary in nuclear magnetic resonance diffusion experiments[END_REF][START_REF] Kiselev | Fundamentals of diffusion MRI physics[END_REF][START_REF] Novikov | Quantifying brain microstructure with diffusion MRI: Theory and parameter estimation[END_REF], 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. [START_REF] Wedeen | The geometric structure of the brain fiber pathways[END_REF], we show that our theoretical analysis and numerical computations are in excellent agreement with the experimental data. The characteristic stretchedexponential 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.

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 [START_REF] Torrey | Bloch equations with diffusion terms[END_REF]:

∂M ∂t = (D 0 ∇ 2 -iγgx)M , (1)
where we chose the x-axis as the gradient direction. The differential operator on the right-hand side of the equation, D 0 ∇ 2iγ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 ∇ 2 . However, the iγgx term makes the whole operator non-Hermitian, which affects deeply its spectral properties [START_REF] Stoller | Transverse spin relaxation in inhomogeneous magnetic fields[END_REF][START_REF] Grebenkov | The complex Airy operator on the line with a semipermeable barrier[END_REF][START_REF] Grebenkov | On spectral properties of the Bloch-Torrey operator in two dimensions[END_REF][START_REF] Almog | Spectral semiclassical analysis of a complex Schrödinger operator in exterior domains[END_REF]: 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 [START_REF] Helffer | Spectral theory and its applications[END_REF][START_REF] Moiseyev | Non-Hermitian quantum mechanics[END_REF]. 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 [START_REF] Grebenkov | NMR survey of reflected Brownian motion[END_REF][START_REF] Stoller | Transverse spin relaxation in inhomogeneous magnetic fields[END_REF][START_REF] De Swiet | Decay of nuclear magnetization by bounded diffusion in a constant field gradient[END_REF][START_REF] Hürlimann | Spin echoes in a constant gradient and in the presence of simple restriction[END_REF][START_REF] Grebenkov | Exploring diffusion across permeable barriers at high gradients. II. Localization regime[END_REF][START_REF] Grebenkov | Diffusion MRI/NMR at high gradients: Challenges and perspectives[END_REF][START_REF] Herberthson | Dynamics of local magnetization in the eigenbasis of the Bloch-Torrey operator[END_REF]. 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 [START_REF] Ziener | Spin dephasing in a magnetic field[END_REF][START_REF] Ziener | Spin echoes: full numerical solution and breakdown of approximative solutions[END_REF] and references therein).

Eigenmode decomposition

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

In bounded domains, the term iγ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 [START_REF] Grebenkov | On spectral properties of the Bloch-Torrey operator in two dimensions[END_REF]). 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" [START_REF] Stoller | Transverse spin relaxation in inhomogeneous magnetic fields[END_REF] 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 (r) is replaced by its complex conjugate v * m (r)):

Ω v n (r)v m (r) dr =        0 if n m 1 if n = m . (2)
This allows one to project functions onto the {v n } n≥1 basis by using standard formulas. Let us assume that the initial condition is a uniform magnetization after the exciting 90 • radio frequency (RF) pulse, i.e. M(r, t = 0) = 1/V, where V is the volume of the domain Ω. Hence we decompose the transverse magnetization as:

M(r, t) = 1 √ V n≥1 µ n v n (r)e -λ n t , (3)
where µ n are complex coefficients given by

µ n = 1 √ V Ω v n (r) dr , (4)
and λ n are the corresponding eigenvalues:

-(D 0 ∇ 2 -iγgx)v n = λ n v n . (5)
We set the minus sign in front of the Bloch-Torrey operator to have eigenvalues with a positive real part, thus we sort λ n by increasing real part.

At the end of the first gradient pulse (i.e., at t = δ), a refocusing 180 • RF pulse is applied, which is equivalent to applying a complex conjugation to the magnetization M:

M(r, δ + 0) = 1 √ V n≥1 µ * n v * n (r)e -λ * n δ = 1 √ V n≥1 µ * n e -λ * n δ m≥1 β n,m v m (r) , (6)
where we have decomposed each v * n over the eigenmode basis (v m ) m≥1 , with complex coefficients β n,m :

β n,m = Ω v * n (r)v m (r) dr . (7)
Finally, after the second gradient pulse, the transverse magnetization is equal to

M(r, 2δ) = 1 √ V n,m≥1 µ * n β n,m v m (r)e -(λ * n +λ m )δ , (8)
and the resulting signal can be computed from

E = Ω M(r, 2δ) dr = n,m≥1 µ * n β n,m µ m e -(λ * n +λ m )δ . (9)
We stress again that the eigenvalues λ n and eigenmodes v n (r) depend on the gradient g, and so do the coefficients µ n and β n,m .

Localization regime and gradient length

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

M(r, 2δ) ≈ 1 √ V µ * 1 β 1,1 v 1 (r)e -2Re(λ 1 )δ , (10)
E ≈ C 1,1 e -2Re(λ 1 )δ , (11)
where the coefficient (1, 1)) are small enough. Note that the eigenvalues may come in complex conjugated pairs, in which case the above expressions [START_REF] Grebenkov | Diffusion MRI/NMR at high gradients: Challenges and perspectives[END_REF] and [START_REF] Grebenkov | Exploring diffusion across permeable barriers at high gradients. II. Localization regime[END_REF] would involve four terms (two diagonal terms and two cross terms) with similar exponential decays e -2Re(λ 1 )δ . We show in Sec. 2.3 that these additional terms may produce oscillations in the signal on top of the overall decay [START_REF] Grebenkov | Exploring diffusion across permeable barriers at high gradients. II. Localization regime[END_REF]. For clarity, we first focus on Eq. ( 11) in this subsection.

C 1,1 = |µ 1 | 2 β 1,
In this regime, the behavior of the magnetization, and thus the signal, is then dictated by the dependence of Re(λ 1 ) on the gradient g. In [START_REF] Stoller | Transverse spin relaxation in inhomogeneous magnetic fields[END_REF][START_REF] De Swiet | Decay of nuclear magnetization by bounded diffusion in a constant field gradient[END_REF][START_REF] Grebenkov | Exploring diffusion across permeable barriers at high gradients. II. Localization regime[END_REF][START_REF] Grebenkov | The complex Airy operator on the line with a semipermeable barrier[END_REF][START_REF] Grebenkov | On spectral properties of the Bloch-Torrey operator in two dimensions[END_REF], an asymptotic expansion for the λ n at high gradients was derived. In particular, the first term of the expansion of the first eigenvalue is universal:

2Re(λ 1 ) = |a ′ 1 | D 0 ℓ 2 g + O(ℓ -3/2 g ) , (g → ∞) (12)
where ℓ g = (γg/D 0 ) -1/3 is the gradient length and a ′ 1 ≈ -1.02 is the first zero of the derivative of the Airy function Ai(z). This formula is universal in the sense that it does not depend on the geometry of the domain. By combining Eqs. [START_REF] Grebenkov | Diffusion MRI/NMR at high gradients: Challenges and perspectives[END_REF] and [START_REF] Grebenkov | The complex Airy operator on the line with a semipermeable barrier[END_REF], one gets the stretchedexponential behavior, which is a hallmark of the localization regime, ln(E) ∝ (ℓ δ /ℓ g ) 2 ∝ (bD 0 ) 1/3 , (bD 0 ≫ 1) [START_REF] Callaghan | Diffraction-like effects in NMR diffusion studies of fluids in porous solids[END_REF] where ℓ δ = √ D 0 δ is the diffusion length during the gradient pulse, and b = 2γ 2 g 2 δ 3 /3 for this gradient profile.

In order to offer an interpretation for the gradient length ℓ 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 √ D 0 t, and thus have accumulated a random phase difference of the order of (γgt) √ D 0 t. The nuclei have practically uncorrelated phases when this difference is comparable with 2π, which yields a decorrelation time: t g = (γ 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 ℓ g is then the typical distance traveled by these particles during time t g : ℓ g = D 0 t g .

If the gradient length ℓ g is much smaller than the typical size of the medium ℓ s and the diffusion length ℓ δ , 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 ℓ g close to the boundary: this is the localization regime. In regard to Eq. [START_REF] Grebenkov | Diffusion MRI/NMR at high gradients: Challenges and perspectives[END_REF], the localization regime corresponds to the localization of the first eigenmodes of the Bloch-Torrey operator at high gradients over a width ℓ g [START_REF] Grebenkov | Exploring diffusion across permeable barriers at high gradients. II. Localization regime[END_REF][START_REF] Grebenkov | Diffusion MRI/NMR at high gradients: Challenges and perspectives[END_REF][START_REF] Grebenkov | On spectral properties of the Bloch-Torrey operator in two dimensions[END_REF]. This localization occurs where the gradient is orthogonal to the boundary, i.e. where the boundary prevents diffusion along the gradient direction the most.

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 [START_REF] Grebenkov | Exploring diffusion across permeable barriers at high gradients. II. Localization regime[END_REF][START_REF] Grebenkov | The complex Airy operator on the line with a semipermeable barrier[END_REF][START_REF] Grebenkov | On spectral properties of the Bloch-Torrey operator in two dimensions[END_REF]. 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 [START_REF] Grebenkov | On spectral properties of the Bloch-Torrey operator in two dimensions[END_REF] 2Re

(λ k,l ) = |a ′ k | D 0 ℓ 2 g + (2l -1) D 0 |R c | 1/2 ℓ 3/2 g - √ 3D 0 2|a ′ k |R c ℓ g + O(ℓ -1/2 g ) , (14)
where a ′ 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

2Re(λ k,l 1 ,l 2 ) = |a ′ k | D 0 ℓ 2 g + D 0 ℓ 3/2 g 2l 1 -1 |R c,1 | 1/2 + 2l 2 -1 |R c,2 | 1/2 - √ 3D 0 2|a ′ k |ℓ g 1 R c,1 + 1 R c,2 + O(ℓ -1/2 g ) , (15)
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 [START_REF] Niendorf | Biexponential diffusion attenuation in various states of brain tissue -implications for diffusionweighting imaging[END_REF] 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 /(|R c | 1/2 ℓ 3/2 g ) which means that although the localization regime emerges when ℓ δ ≫ ℓ g , the asymptotic decay ( 11) is established at a longer time:

ℓ δ ≫ |R c | 1/4 ℓ 3/4
g (assuming that ℓ g < |R c | in the localization regime, which is always the case for basic geometries). In other words, in the intermediate regime

ℓ g ≪ ℓ δ ≪ |R c | 1/4 ℓ 3/4
g , the transverse magnetization is localized at the boundary of the domain, but the truncated formulas [START_REF] Grebenkov | Diffusion MRI/NMR at high gradients: Challenges and perspectives[END_REF] and [START_REF] Grebenkov | Exploring diffusion across permeable barriers at high gradients. II. Localization regime[END_REF] 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 δ, 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 ℓ g is much smaller than L. In fact, as we stated in Sec. 2.2, ℓ 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 ℓ g < 0.38L. Due to the left-right symmetry (x → -x), the first two eigenmodes v 1 (x) and v 2 (x) satisfy the identity v 2 (x) = v * 1 (-x) and λ 2 = λ * 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 δ, the magnetization may be represented by a superposition of v 1 and v 2 :

M(x, 2δ) ≈ c 1,1 + c 2,1 e -2iIm(λ 1 )δ e -2Re(λ 1 )δ v 1 (x) + c 2,2 + c 1,2 e 2iIm(λ 1 )δ e -2Re(λ 1 )δ v 2 (x) , (16)
with c n,m = µ * n β n,m / √ V, and the signal is given by

E ≈ Ce -2Re(λ 1 )δ , (17)
with

C = 2(C 1,1 + Re(C 1,2 )) , (18)
where C 1,1 was previously introduced in Eq. ( 11) and C 1,2 is an oscillating function of g and δ given by

C 1,2 = µ * 1 β 1,2 µ 2 e 2iIm(λ 1 )δ . ( 19
)
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 [START_REF] Stoller | Transverse spin relaxation in inhomogeneous magnetic fields[END_REF]:

2Im(λ 1 ) = ±       γgL - √ 3|a ′ 1 | D 0 ℓ 2 g       + O(ℓ -3/2 g ) , (20)
where the ± sign is merely a matter of convention (given that λ 2 = λ * 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 [START_REF] Stoller | Transverse spin relaxation in inhomogeneous magnetic fields[END_REF][START_REF] Grebenkov | Exploring diffusion across permeable barriers at high gradients. II. Localization regime[END_REF][START_REF] Grebenkov | The complex Airy operator on the line with a semipermeable barrier[END_REF].

The conclusion of this computation is that the crossterm C 1,2 produces oscillations in the signal on top of the asymptotic decay [START_REF] Grebenkov | Exploring diffusion across permeable barriers at high gradients. II. Localization regime[END_REF]. The definition [START_REF] Stoller | Transverse spin relaxation in inhomogeneous magnetic fields[END_REF] of β 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 ℓ g of the modes, and their spacing L: the smaller the ratio ℓ g /L, the smaller the overlapping, and thus, the smaller the oscillating term. In the limit of well-separated modes, β 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.

The influence of ∆ -δ

Until now we have considered a PGSE sequence with ∆ = δ, 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, ∆ -δ. Mathematically, the effect of this diffusion step is to multiply the transverse magnetization just after the first gradient pulse by the operator D = exp (∆ -δ)D 0 ∇ 2 . For simplicity we assume that the 180 • 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), [START_REF] Hürlimann | Spin echoes in a constant gradient and in the presence of simple restriction[END_REF], and the consequent analysis remain applicable with the new definition for the coefficients β n,m :

β n,m = Ω (Dv * n )(r)v m (r) dr . (21)
The dependence on ∆ -δ is now hidden in the coefficients β n,m . Note that if ∆ = δ, then the operator D is the identity (since diffusion during zero time does not affect the magnetization) and we recover the previous definition [START_REF] Stoller | Transverse spin relaxation in inhomogeneous magnetic fields[END_REF] for β n,m . If we consider the regime of long δ 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 [START_REF] Almog | Spectral semiclassical analysis of a complex Schrödinger operator in exterior domains[END_REF] of the coefficient β 1,1 . In turn, the global amplitude of the signal in Eq. ( 11) is also affected by this coefficient and decreases with increasing ∆ -δ.

The situation of two overlapping eigenmodes is more complex. The diffusion step also changes the values of the coefficients β 1,2 , β 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 ℓ 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 [START_REF] Grebenkov | Exploring diffusion across permeable barriers at high gradients. II. Localization regime[END_REF].

Unbounded domains

In the previous subsections, we assumed that diffusion takes place in a bounded domain Ω, 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 [START_REF] Almog | Spectral semiclassical analysis of a complex Schrödinger operator in exterior domains[END_REF], even though the spectrum of the Laplacian operator ∇ 2 is continuous. In other words, the iγ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 [START_REF] Grebenkov | On spectral properties of the Bloch-Torrey operator in two dimensions[END_REF], 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 [START_REF] Grebenkov | Diffusion MRI/NMR at high gradients: Challenges and perspectives[END_REF]. In the following sections, we investigate experimentally the emergence of the localization regime in such a domain (array of rods). Fig. 1: The real part of the transverse magnetization M(x) in a slab for two different gradient strengths, in the long δ regime, as well as the real part of the first two eigenmodes v 1 (x) and v 2 (x) weighted by the coefficients in front of v 1 (x) and v 2 (x) in Eq. ( 16). We checked that Eq. ( 16), i.e. the superposition of v 1 (x) and v 2 (x), reproduces perfectly the exact magnetization M(x) in this regime. The computations were performed with a matrix formalism (described in Sec. 3). In both cases, one observes the localization of the magnetization near the endpoints of the interval. We chose a constant ratio ℓ δ /ℓ g = 2.5 for both figures, so that the amplitude of the magnetization is approximately the same (see Eq. ( 13)). (left) ℓ g /L ≈ 0.2: one can see some overlapping between the two eigenmodes. (right) ℓ g /L ≈ 0.1: there is almost no overlapping of the eigenmodes.

Material and Methods

Experiments

rier 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) [START_REF] Walker | Spin-exchange optical pumping of noble-gas nuclei[END_REF][START_REF] Oros | Hyperpolarized xenon in NMR and MRI[END_REF][START_REF] Shah | Measurement of rubidium and xenon absolute polarization at high temperatures as a means of improved production of hyperpolarized 129 Xe[END_REF]. For technical reasons [START_REF] Walker | Spin-exchange optical pumping of noble-gas nuclei[END_REF], a gas mixture (Air Liquide Deutschland GmbH, Dsseldorf, 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 ± 0.2) • 10 -5 m 2 s -1 [START_REF] Kuder | Diffusion pore imaging by hyperpolarized Xenon-129 nuclear magnetic resonance[END_REF], which is one order of magnitude larger than for pure xenon gas [START_REF] Acosta | Diffusion in binary gas mixtures studied by NMR of hyperpolarized gases and molecular dynamics simulations[END_REF]. The gas was transferred at a small constant flow (approximately 150 mL/min) to the phantom positioned in an in-house built xenon coil in the isocenter of the magnet of a 1.5 T clinical MR scanner (Magnetom Symphony, A Tim System, Siemens Healthcare, Erlangen, Germany) with a maximal employed gradient amplitude of 32 mT/m. The experimental set-up and hyperpolarization process are detailed in [START_REF] Kuder | Diffusion pore imaging by hyperpolarized Xenon-129 nuclear magnetic resonance[END_REF][START_REF] Demberg | Nuclear magnetic resonance diffusion pore imaging: Experimental phase detection by double diffusion encoding[END_REF]. 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 mm, built by the inhouse 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 centerto-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 mm × 42.3 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 highgradient 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 mm thick- ness 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 mm). For phantoms 4 and 5, measurements were taken with the gradient vector pointing along the leftright 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 δ of the trapezoidal gradient pulses were set to 6 ms and include flat top time plus the ramp-up time of ǫ = 0.32 ms. The gradient separation time was ∆ = 9.34 ms. The diffusion-weighted signal was sampled by gradually increasing the gradient amplitude in 32 steps from 0 to g max = 32 mT/m, recording the spin echo signal and acquiring up to 15 averages. The time between two consecutive 90 • excitation pulses, i.e. the time between two measurements, was set to 18 s to restore the polariza- tion 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 prereadout directly after the 90 • 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 [START_REF] Caprihan | A multiplenarrow-pulse approximation for restricted diffusion in a timevarying field gradient[END_REF][START_REF] Callaghan | A simple matrix formalism for spin echo analysis of restricted diffusion under generalized gradient waveforms[END_REF][START_REF] Barzykin | Exact solution of the Torrey-Bloch equation for a spin echo in restricted geometries[END_REF][START_REF] Grebenkov | NMR survey of reflected Brownian motion[END_REF], 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 [START_REF] Grebenkov | Laplacian eigenfunctions in NMR. I. A numerical tool[END_REF][START_REF] Grebenkov | Pulsed-gradient spin-echo monitoring of restricted diffusion in multilayered structures[END_REF] for details). The matrix formalism was also used to compute the transverse magnetization in these two domains (see similar computations in [START_REF] Grebenkov | Exploring diffusion across permeable barriers at high gradients. II. Localization regime[END_REF]). 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 • 10 7 random walkers and 1.6 • 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 × 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 ℓ g and ℓ δ 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.

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 [START_REF] Hürlimann | Spin echoes in a constant gradient and in the presence of simple restriction[END_REF], whereas only a few theoretical studies were devoted to the cylinder geometry [START_REF] De Swiet | Decay of nuclear magnetization by bounded diffusion in a constant field gradient[END_REF][START_REF] Grebenkov | Diffusion MRI/NMR at high gradients: Challenges and perspectives[END_REF][START_REF] Grebenkov | On spectral properties of the Bloch-Torrey operator in two dimensions[END_REF]. The signal for an array of rods was previously computed numerically in [START_REF] Grebenkov | Diffusion MRI/NMR at high gradients: Challenges and perspectives[END_REF] using a finite-element method [START_REF] Nguyen | A finite element method to solve the Bloch-Torrey equation applied to diffusion magnetic resonance imaging[END_REF].

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

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

E = exp -bD 0 sin 2 θ E slab (L, g cos θ) , ( 22
)
where θ is the angle between the gradient and the z-axis and E 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 ℓ g and ℓ δ (about 40 mm). One gets the slowest decay by ensuring that the gradient is orthogonal to the slab, i.e. θ = 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 [START_REF] Mitra | Diffusion propagator as a probe of the structure of porous media[END_REF][START_REF] Mitra | Short-time behavior of the diffusion coefficient as a geometrical probe of porous media[END_REF][START_REF] De Swiet | Decay of nuclear magnetization by bounded diffusion in a constant field gradient[END_REF][START_REF] Moutal | Probing surfaceto-volume ratio of an anisotropic medium by diffusion NMR with general gradient encoding[END_REF] yields an apparent diffusion coefficient

D ≈ D 0 1 -η 4 3 √ π S V D 0 T , (23)
where S /V is the surface-to-volume ratio of the confining domain (S /V = 2/L for a slab), η ≈ 0.9 is a numerical prefactor that depends on the sequence, and T = ∆+δ is the duration of the gradient sequence. This is a short-time approximation in the sense that η √ 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 ≈ 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. [START_REF] Grebenkov | NMR survey of reflected Brownian motion[END_REF][START_REF] Stoller | Transverse spin relaxation in inhomogeneous magnetic fields[END_REF][START_REF] Hürlimann | Spin echoes in a constant gradient and in the presence of simple restriction[END_REF][START_REF] De Swiet | Decay of nuclear magnetization by bounded diffusion in a constant field gradient[END_REF][START_REF] Grebenkov | Diffusion MRI/NMR at high gradients: Challenges and perspectives[END_REF]) The signal for free diffusion e -bD 0 is indicated by a dotted line, whereas the low-b, short-time approximation e -bD , where D is given by Eq. ( 23), is plotted as a dashed-dotted line. The high-gradient asymptotic regime (24) appears as a dashed line.

E slab ≈ 2C 1,1 exp -|a ′ 1 |ℓ 2 δ /ℓ 2 g , (24)
where the prefactor C 1,1 can be computed exactly [START_REF] Hürlimann | Spin echoes in a constant gradient and in the presence of simple restriction[END_REF][START_REF] Grebenkov | Exploring diffusion across permeable barriers at high gradients. II. Localization regime[END_REF] and scales as ℓ g /L. One can interpret this prefactor as the fraction of spins inside the two layers of thickness ∼ ℓ g where the signal is localized. As we discussed in Sec. 2.4, C 1,1 decreases with increasing diffusion step duration ∆ -δ. Note that the localization regime emerges at gradients as small as 10 mT/m. Since the width of the slab is much greater than ℓ 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 [START_REF] Herberthson | Dynamics of local magnetization in the eigenbasis of the Bloch-Torrey operator[END_REF]. 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. θ slightly different from 0 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 ǫ = 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 (ǫ = 0.34 ms).

Diffusion inside a cylinder

If the diameter of the cylinder is much larger than the gradient length ℓ g and the diffusion length ℓ δ , 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 nor- mal vector to the boundary is close to 0 (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 mm. Although the signal decays faster than in the slab, one observes a similar stretchedexponential behavior. Using Eqs. ( 14) and [START_REF] Laun | Determination of the defining boundary in nuclear magnetic resonance diffusion experiments[END_REF], one gets the asymptotic decay for the cylinder:

E cyl ≈ C exp        -|a ′ 1 | ℓ 2 δ ℓ 2 g - ℓ 2 δ R 1/2 ℓ 3/2 g + √ 3ℓ 2 δ 2|a ′ 1 |Rℓ g        , (25)
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 (ℓ 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 mm, see Fig. 7), the signal shows some oscillations that are reminiscent of diffusion-diffraction patterns [START_REF] Callaghan | Diffraction-like effects in NMR diffusion studies of fluids in porous solids[END_REF][START_REF] Callaghan | Diffusion in porous systems and the influence of pore morphology in pulsed gradient spin-echo nuclear magnetic resonance studies[END_REF][START_REF] Callaghan | Pulsed-gradient spin-echo NMR for planar, cylindrical, and spherical pores under conditions of wall relaxation[END_REF][START_REF] Linse | The validity of the short-gradientpulse approximation in NMR studies of restricted diffusion. Simulations of molecules diffusing between planes, in cylinders, and spheres[END_REF][START_REF] Callaghan | A simple matrix formalism for spin echo analysis of restricted diffusion under generalized gradient waveforms[END_REF]. This is the case discussed in Sec. 2.3 where two localization pockets overlap because ℓ g /(2R) is not small enough. The signal is still given by Eq. ( 25) but one The signal for free diffusion is indicated by a dotted line, whereas the short-time approximation e -bD , where D is given by Eq. ( 23), is plotted as a dashed-dotted line.

The high-gradient asymptotic regime (25) appears as a dashed line.

cannot neglect the cross-term C 1,2 in the expression [START_REF] Novikov | Quantifying brain microstructure with diffusion MRI: Theory and parameter estimation[END_REF] 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 [START_REF] Torrey | Bloch equations with diffusion terms[END_REF] and [START_REF] Grebenkov | On spectral properties of the Bloch-Torrey operator in two dimensions[END_REF] with L = 2R, where the coefficients µ 1 , µ 2 , and β 1,2 were computed numerically from the eigenmodes. We stress that these coefficients generally depend on g and that β 1,2 additionally depends on the diffusion step duration ∆-δ. 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.

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 mm and 2R = 3.2 mm) with the gradient vector in the diagonal direction; (b) phantom 5 (L = 3.4 mm and 2R = 3.2 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 pock between two neighboring rods along the gradient direction, i.e. the spacing between two neighboring localization pockets (see Fig. 9): e pock = e d = √ 2L -2R = 2.5 mm in (a); e pock = e d = 1.6 mm in (b); and e pock = e h = 0.2 mm in (c). Figure 10 shows the signal for these three cases, ordered by descending e 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 surfaceto-volume ratio of the structure is too large so that the approximate formula [START_REF] Moiseyev | Non-Hermitian quantum mechanics[END_REF] 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 [START_REF] Grebenkov | On spectral properties of the Bloch-Torrey operator in two dimensions[END_REF]:

E rods ≈ C exp        -|a ′ 1 | ℓ 2 δ ℓ 2 g - ℓ 2 δ R 1/2 ℓ 3/2 g - √ 3ℓ 2 δ 2|a ′ 1 |Rℓ g        , ( 26
)
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 ℓ δ and ℓ g . Here the same phenomenon occurs. With e pock much larger than ℓ 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 pock and ℓ 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 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 pock is smaller than ℓ g even at the highest gradient available. If one could approximate the small space between two neighboring rods as a slab, then the ratio ℓ g /e 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 pock between the rods (see [START_REF] Robertson | Spin-echo decay of spins diffusing in a bounded region[END_REF][START_REF] Grebenkov | NMR survey of reflected Brownian motion[END_REF]).

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 ap- pears whenever the gradient length ℓ g = (γg/D 0 ) -1/3 is much smaller than the diffusion length ℓ δ = (D 0 δ) 1/2 and any relevant geometrical length scale ℓ 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) ℓ g ≪ ℓ δ and (ii) ℓ g ≪ ℓ 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 δ 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 stretchedexponential decay of the signal with g.

An extreme case where condition (i) is not satisfied would be narrow gradient pulses experiments (g → ∞ and δ → 0). Although these experiments require high gradients, they do not achieve the localization regime. In fact, whereas the localization regime emerges when (γg) 2 δ 3 D 0 ≫ 1, narrow-pulse experiments correspond to γg → ∞ and δ → 0 such that γgδ = q is a finite value. This leads to (γg) 2 δ 3 D 0 → 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 ∆ = 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 → 0 and δ → ∞). 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 [START_REF] Callaghan | Diffraction-like effects in NMR diffusion studies of fluids in porous solids[END_REF][START_REF] Callaghan | Diffusion in porous systems and the influence of pore morphology in pulsed gradient spin-echo nuclear magnetic resonance studies[END_REF][START_REF] Callaghan | Pulsed-gradient spin-echo NMR for planar, cylindrical, and spherical pores under conditions of wall relaxation[END_REF][START_REF] Linse | The validity of the short-gradientpulse approximation in NMR studies of restricted diffusion. Simulations of molecules diffusing between planes, in cylinders, and spheres[END_REF][START_REF] Callaghan | A simple matrix formalism for spin echo analysis of restricted diffusion under generalized gradient waveforms[END_REF]). 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 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 ℓ δ /ℓ g ≫ 1 which can be restated as bD 0 ∝ γ 2 D 0 g 2 δ 3 ≫ 1. In order to rescale the experimental conditions from xenon gas to water, we compute the ratio (γ 2 D 0 ) xenon /(γ 2 D 0 ) water ≈ 10 3 . This means that in order to have the same value of γ 2 D 0 g 2 δ 3 , one has to increase g and δ such that g 2 δ 3 is 10 3 times larger with water than with xenon. For example, the experiments by Hürlimann et al were performed with δ = 60 ms which is approximately 10 times longer than in our experiments, and gradients of comparable magnitude as ours (around 20 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 [START_REF] Mitra | Diffusion propagator as a probe of the structure of porous media[END_REF][START_REF] Mitra | Short-time behavior of the diffusion coefficient as a geometrical probe of porous media[END_REF][START_REF] Moutal | Probing surfaceto-volume ratio of an anisotropic medium by diffusion NMR with general gradient encoding[END_REF], the structural organization of barriers is accessible at long times [START_REF] Novikov | Random walks with barriers[END_REF][START_REF] Novikov | Revealing mesoscopic structural universality with diffusion[END_REF]. In biomedical applications, the apparent diffusion coefficient and kurtosis are employed as biomarkers to detect stroke, tumors, lesions, partial tissue destruction, etc. (see [START_REF] Novikov | Quantifying brain microstructure with diffusion MRI: Theory and parameter estimation[END_REF] 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 [START_REF] Almog | Spectral semiclassical analysis of a complex Schrödinger operator in exterior domains[END_REF]. 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 [START_REF] Stoller | Transverse spin relaxation in inhomogeneous magnetic fields[END_REF]. 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 accen-tuating the differences from pure diffusion.

  were performed using hyperpolarized xenon-129 gas (γ ≈ 74 • 10 6 s -1 T -1 ) continuously flowing through phantoms containing different diffusion bar-

  rods 17 × 10 rods 2R = 3.2 mm; L = 4 mm b V = 1.48 mm; b H = 1.35 mm 5 array of rods 20 × 10 rods 2R = 3.2 mm; L = 3.4 mm b V = 0.78 mm; b H = 1.08 mm

Fig. 2 :

 2 Fig. 2: Phantoms with corresponding magnified schematic depiction showing relevant length scales. (a) Phantom 1: Parallel plates. (b) Phantoms 2 and 3: Cylindrical tubes. (c) Phantoms 4 and 5: Cylindrical rods on a square grid.

Fig. 3 :

 3 Fig. 3: PGSE sequence. The spin echo forms at time TE.

Fig. 4 :

 4 Fig.4: Signal attenuation for phantom 1 (a slab of width 3 mm). Experimental results are shown by full circles and matrix formalism computation by a solid line. The signal for free diffusion e -bD 0 is indicated by a dotted line, whereas the low-b, short-time approximation e -bD , where D is given by Eq. (23), is plotted as a dashed-dotted line. The high-gradient asymptotic regime (24) appears as a dashed line.

Fig. 5 :

 5 Fig. 5: Transverse magnetization computed by matrix formalism inside phantom 2 (a cylinder of diameter 3.8 mm) for four values of the gradient g: (a) 2 mT/m, (b) 5 mT/m, (c) 10 mT/m and (d) 32 mT/m. The direction of the gradient is indicated by an arrow.

Fig. 6 :

 6 Fig.6: Signal attenuation for diffusion inside phantom 2 (a cylinder of diameter 3.8 mm). Experimental results are shown by full circles and matrix formalism computation by a solid line. The signal for free diffusion is indicated by a dotted line, whereas the short-time approximation e -bD , where D is given by Eq. (23), is plotted as a dashed-dotted line. The high-gradient asymptotic regime (25) appears as a dashed line.

Fig. 7 :

 7 Fig.7: Signal attenuation for diffusion inside phantom 3 (a cylinder of diameter 2 mm). Experimental results are shown by full circles and matrix formalism computation by a solid line. The signal for free diffusion is indicated by a dotted line, whereas the short-time approximation e -bD , where D is given by Eq. (23), is plotted as a dashed-dotted line. The high-gradient asymptotic regime (25) appears as a dashed line.

Fig. 8 :

 8 Fig. 8: Transverse magnetization computed by matrix formalism inside phantom 3 (a cylinder of diameter 2 mm) for four values of the gradient g: (a) 2 mT/m, (b) 5 mT/m, (c) 10 mT/m and (d) 32 mT/m. The direction of the gradient is indicated by an arrow.

Fig. 9 :

 9 Fig. 9: Schematic representation of the rods showing the lengths e d and e h .

Fig. 10 :

 10 Fig. 10: Signal attenuation for diffusion in phantoms 4 and 5 (array of rods of diameter 3.2 mm and center-to-center spacing 4 mm and 3.4 mm, respectively). (a) phantom 4 with the gradient in the diagonal direction, spacing e pock = 2.5 mm; (b) phantom 5 and diagonal gradient direction, spacing e pock = 1.6 mm; (c) phantom 5 and horizontal gradient direction, spacing e pock = 0.2 mm. Experimental results are shown by full circles and Monte Carlo simulations by a solid line. The signal for free diffusion is indicated by a dotted line. The asymptotic regime (26) appears as a dashed line. It is not shown for (c) as this regime is not applicable here, see the text.

Table 1 :

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