# Spectral properties of the Bloch–Torrey operator in three dimensions

Denis S Grebenkov

<figure><img src="image_1.png" /></figure>

► To cite this version:

Denis S Grebenkov. Spectral properties of the Bloch–Torrey operator in three dimensions. Journal of Physics A: Mathematical and Theoretical, 2024, 57 (12), pp.125201. 10.1088/1751-8121/ad2d6d. hal-04749296

HAL Id: hal-04749296

https://hal.science/hal-04749296v1

Submitted on 23 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
Spectral properties of the Bloch-Torrey operator in
three dimensions

Denis S. Grebenkov

E-mail: denis.grebenkov@polytechnique.edu

Laboratoire de Physique de la Matière Condensée,

CNRS – Ecole Polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau,
France

**Abstract.** We consider the Bloch-Torrey operator, $-\Delta + igx$, that governs the time evolution of the transverse magnetization in diffusion magnetic resonance imaging (dMRI). Using the matrix formalism, we compute numerically the eigenvalues and eigenfunctions of this non-Hermitian operator for two bounded three-dimensional domains: a sphere and a capped cylinder. We study the dependence of its eigenvalues and eigenfunctions on the parameter $g$ and on the shape of the domain (its eventual symmetries and anisotropy). In particular, we show how an eigenfunction drastically changes its shape when the associated eigenvalue crosses a branch (or exceptional) point in the spectrum. Potential implications of this behavior for dMRI are discussed.

PACS numbers: 02.50.-r, 05.40.-a, 02.70.Rr, 05.10.Gg

*Keywords:* Bloch-Torrey operator, non-Hermitian operator, branch point, diffusion-weighted NMR, localization, pulsed-gradient spin-echo, microstructure

Submitted to: *J. Phys. A: Math. Gen.*
# 1. Introduction

Diffusion magnetic resonance imaging (dMRI) is a non-invasive technique with multiple applications in medicine, neurosciences and material sciences [1–4]. In a typical setting, a static magnetic field $B_0$ is applied along the $z$ axis to create the local magnetization of the nuclei (e.g., protons). A radio-frequency (rf) 90° pulse allows one to turn the local magnetization into the transverse $xy$ plane, in which it starts to precess around the $z$ axis. If the static field $B_0$ is superimposed with a spatially inhomogeneous magnetic field, the Larmor frequency of each precessing nucleus depends on its spatial location, allowing one to encode random trajectories of these nuclei that are hindered by the environment and thus contain potentially exploitable information on its structural properties. Many theoretical and numerical approaches have been developed to study this fundamental problem (see reviews [5–7] and references therein). The most common microscopic description of this phenomenon relies on the Bloch-Torrey equation [8] that governs time evolution of the transverse magnetization $m(x, t)$ of the nuclei in a confining domain $\Omega \subset \mathbb{R}^d$:

$$ \partial_t m(x, t) = D_0 \Delta m(x, t) - i\gamma(\mathbf{G}(t) \cdot \mathbf{x})m(x, t) \quad (x \in \Omega), \quad (1) $$

where $D_0$ is the constant (self-)diffusion coefficient of the nuclei (e.g., water molecules), $\Delta$ is the Laplace operator, $\gamma$ is the gyromagnetic ratio of the nuclei, and $\mathbf{G}(t)$ is the gradient profile of the applied magnetic field, which is set and controlled by the experimental setup. The Bloch-Torrey equation is usually complemented by the uniform initial condition, $m(x, 0) = m_0 = 1/|\Omega|$, reflecting the homogeneous excitation of the nuclei at time $t = 0$ by the rf pulse in the volume $|\Omega|$ of the confining domain $\Omega$. The confining microstructure is incorporated via an appropriate boundary condition. A typical situation of an impenetrable inert surface $\partial\Omega$ is described by Neumann boundary condition, $\partial_n m(x, t)|_{\partial\Omega} = 0$, stating that the magnetization flux across the surface is zero, where $\partial_n$ is the normal derivative oriented outwards the confining domain $\Omega$. Surface relaxation due to magnetic impurities on the boundary or nuclear exchange across permeable membranes can also be described by modifying the boundary condition [9, 10]. In addition, $T_1$ and $T_2$ bulk relaxation mechanisms can be included into Eq. (1). Since the transverse magnetization in any point $x$ is too small to be measured, only its integral over the confining domain (or a voxel) is accessible in experiments:

$$ S = \int_{\Omega} d\mathbf{x} m(\mathbf{x}, t). \quad (2) $$

This macroscopic signal that can be accessed as a function of the gradient profile $\mathbf{G}(t)$, aggregates the microstructural features in a very sophisticated way through the boundary condition to the Bloch-Torrey equation (1). The imaginary unit $i$ in front of the last term of Eq. (1), which represents precession of the nuclei in the transverse plane, makes this classical diffusion-reaction problem challenging. In fact, the differential operator governing time evolution is not Hermitian that results in numerous unexpected
features such as the failure of perturbative approaches at high gradients, localization near specific points on the boundary, or branch points in the spectrum [5, 6, 11].

In order to understand the intricate relation between the microstructure and the signal, one can focus on piecewise constant gradient profiles and study the magnetization evolution during one constant gradient pulse, i.e., to set $G(t) = G$. Denoting by $x$ the coordinate axis in the direction of the gradient, one has $(G \cdot x) = Gx$, where $G = |G|$ is the gradient amplitude, and $x$ is the projection of $x$ onto the direction of $G$. Introducing the Bloch-Torrey operator as

$$ \mathcal{B}_g = -\Delta + igx \quad (g = \gamma G/D_0), \qquad (3) $$

one can formally solve the Bloch-Torrey equation as $m(x, t) = \exp(-D_0 \mathcal{B}_g t)m_0$. In other words, the effect of a constant gradient pulse is represented by the evolution operator $\exp(-D_0 \mathcal{B}_g t)$. One can also deal with more sophisticated gradient profiles by representing them as a sequence of constant gradient pulses and combining the corresponding evolution operators [12–15]. For instance, in a standard Stejskal-Tanner pulsed-gradient spin echo (PGSE) sequence with two rectangular gradient pulses of duration $\delta$ and opposite directions [16], the signal can be written as

$$ S = \int_{\Omega} d\boldsymbol{x} \left( e^{-D_0 \delta \mathcal{B}_g} e^{-D_0 \delta \mathcal{B}_g} \frac{1}{|\Omega|} \right), \qquad (4) $$

where $e^{-D_0 \delta \mathcal{B}_g}$ represents the evolution from the initial uniform magnetization $m_0 = 1/|\Omega|$ during the first gradient pulse, and $e^{-D_0 \delta \mathcal{B}_g}$ describes the evolution during the second gradient pulse with the opposite direction (for simplicity, we assumed here that the second pulse starts immediately after the first one). When the Bloch-Torrey operator has a discrete spectrum, one can use its eigenvalues $\lambda_j^{(g)}$ and eigenfunctions $v_j^{(g)}$ (enumerated by $j = 1, 2, \dots$) to represent the above signal as [17–20]

$$ S = \sum_{j,j'=1}^{\infty} C_{j,j'}^{(g)} e^{-D_0 \delta(\lambda_j^{(g)} + \lambda_{j'}^{(g)})}, \qquad (5) $$

where the coefficients

$$ C_{j,j'}^{(g)} = \frac{1}{|\Omega|} \left( \int_{\Omega} d\boldsymbol{x} \, v_j^{(-g)}(\boldsymbol{x}) \right) \left( \int_{\Omega} d\boldsymbol{x} \, v_j^{(-g)}(\boldsymbol{x}) v_{j'}^{(g)}(\boldsymbol{x}) \right) \left( \int_{\Omega} d\boldsymbol{x} \, v_{j'}^{(g)}(\boldsymbol{x}) \right) \qquad (6) $$

characterize the overlap between two eigenfunctions $v_j^{(-g)}$ and $v_{j'}^{(g)}$, and their projections onto a constant. As a consequence, the macroscopic signal $S$ and its dependence on the microstructure are fully determined by the spectral properties of the Bloch-Torrey operator. Moreover, when gradient pulses are long and/or strong enough such that $D_0 \delta \text{Re}\{\lambda_1^{(g)}\} \gg 1$, the above expansion can be truncated to few terms, yielding a practical approximation for the signal, as discussed below.

The seminal paper by Stoller, Happer and Dyson provided the first thorough analysis of the Bloch-Torrey operator in one dimension (for an interval and a half-line) [21]. In particular, they showed that the spectrum is discrete, while the eigenvalues
$\lambda_k^{(g)}$ of the Bloch-Torrey operator $B_g$ behave as $\lambda_k^{(g)} \propto g^{2/3} \propto G^{2/3}$ at large $G$, that results in the specific long-time decay of the signal, $\ln S \propto G^{2/3}t$, with unexpected “anomalous” dependence $G^{2/3}$ on the gradient. This behavior is drastically different from the common quadratic dependence, $\ln S \propto G^2$, that appears at small gradients in both slow-diffusion and motional-narrowing regimes [6, 22, 23]. The spectral analysis was later extended to different classes of confining domains, including an arbitrary array of permeable intervals [18, 24, 25], a disk and a sphere [17], bounded planar domains [19, 26], the exterior of compact domains [27, 28], and periodic domains [29, 30]. Most focus was on the large-$G$ asymptotic behavior of the eigenvalues and on the localization of eigenfunctions. Moreover, the whole structure of the spectrum, including the existence of branch points (also known as exceptional or diabolic points), was investigated [31]. The existence of branch points is a peculiar feature of non-Hermitian operators (see, e.g. [32–43] and references therein). The “anomalous” $G^{2/3}$-dependence of $\ln S$ was first confirmed experimentally by Hürlimann et al. for diffusion of water molecules between two parallel planes [44], and later for gas diffusion in cylindrical phantoms [20]. Experimental evidence for the localization regime in biological samples was reported [45].

In this paper, we extend the recent analysis from Ref. [31] that was focused on planar domains, into three dimensions. First, we uncover the behavior of eigenvalues and eigenfunctions of the Bloch-Torrey operator for a sphere. While the spherical confinement is one of the most archetypical models in this field, a systematic study of the spectral properties of $B_g$ in this setting is still missing. In particular, we analyze the dependence of eigenvalues on the gradient and reveal the existence of branch points in the spectrum of $B_g$ for this domain. We also discuss one-mode and two-modes approximations of the macroscopic signal. Second, we analyze the spectrum of the Bloch-Torrey operator for a capped cylinder that exhibits structural anisotropy. We show how the structure of the spectrum depends on the gradient direction, in particular, how the branch points can be tuned experimentally. The structure of the underlying eigenfunctions is discussed.

The paper is organized as follows. In Sec. 2, we recall some basic spectral properties of the Bloch-Torrey operator $B_g$. Section 3 presents the detailed analysis for the case of a sphere; in particular, we discuss the dependence of the eigenvalues on $g$, the branch points in the spectrum, and the drastic change of eigenfunctions at these points. In turn, Sec. 4 focuses on a capped cylinder that exhibits anisotropy and allows us to reveal its impact onto the spectrum. Section 5 concludes the paper by summarizing the main results and presenting their practical implications in diffusion MRI. Appendices contain technical discussions such as the description of the numerical procedure for constructing the spectrum of the Bloch-Torrey operator by using the matrix formalism (Appendix A), the matrix elements for a sphere (Appendix B) and for a capped cylinder (Appendix C), as well as a simple orthogonalization procedure for eigenfunctions with degenerate eigenvalues (Appendix D).
## 2. Summary of basic spectral properties

For a given bounded domain $\Omega$ with a smooth boundary $\partial\Omega$, we are interested in the spectral properties of the Bloch-Torrey operator $\mathcal{B}_g$ defined in Eq. (3). As the parameter $g = \gamma G/D_0$ is determined by the amplitude of the gradient used in diffusion MRI, we mainly focus on positive values $g \ge 0$.

In this section, we remind basic spectral properties of the non-Hermitian Bloch-Torrey operator for $g > 0$ (see further discussion in [31] and references therein). As $igx$ is a bounded perturbation of the (unbounded) Laplace operator, the spectrum is discrete, i.e., there is an infinite sequence of eigenvalues $\lambda_j^{(g)}$ and eigenfunctions $v_j^{(g)}(x)$ satisfying

$$ \mathcal{B}_g v_j^{(g)}(x) = \lambda_j^{(g)} v_j^{(g)}(x) \quad (x \in \Omega), \quad \partial_n v_j^{(g)}(x) = 0 \quad (x \in \partial\Omega). \quad (7) $$

The eigenfunctions are in general complex-valued.

Since the Bloch-Torrey operator is not Hermitian for $g > 0$, the standard scalar product in $L_2(\Omega)$, $(u, v) = \int_{\Omega} d\mathbf{x} u^*(\mathbf{x})v(\mathbf{x})$, is replaced by a bilinear form $\langle u, v \rangle = \int_{\Omega} d\mathbf{x} u(\mathbf{x})v(\mathbf{x})$. In particular, the eigenfunctions $\{v_j^{(g)}\}$ are in general not orthogonal to each other, $(v_j^{(g)}, v_{j'}^{(g)}) \ne 0$, as it would be for Hermitian operators (e.g., for $\mathcal{B}_0$). In turn, one can easily show by the Green's formula that

$$ (\lambda_j^{(g)} - \lambda_{j'}^{(g)}) \langle v_j^{(g)}, v_{j'}^{(g)} \rangle = 0, \quad (8) $$

so that if the eigenvalues $\lambda_j^{(g)}$ and $\lambda_{j'}^{(g)}$ are not equal, then $\langle v_j^{(g)}, v_{j'}^{(g)} \rangle = 0$. It is worth stressing that $\langle v, v \rangle$ is not a norm of $v$; in particular, there exist special values of $g$ (so-called branch points, see below), at which $\langle v_j^{(g)}, v_j^{(g)} \rangle = \int_{\Omega} d\mathbf{x} [v_j^{(g)}]^2 = 0$. In general, however, this integral is not zero, and we normalize the eigenfunctions to have

$$ \langle v_j^{(g)}, v_j^{(g)} \rangle = 1. \quad (9) $$

This condition fixes the normalization up to a factor $\pm 1$.

The eigenvalues are in general complex-valued, with positive real parts that accumulate at $+\infty$. It is therefore convenient to order the eigenvalues according to their increasing real parts. However, we will adopt a different ordering procedure. In fact, the eigenvalues $\lambda_j^{(g)}$ can be understood as different branches in the complex plane $\mathbb{C}$ of a multi-valued function $\lambda(g)$ defined implicitly as the solution of the transcendental equation $\det(\mathcal{B}_g - \lambda(g)\mathcal{I}) = 0$ for any fixed $g$, where $\mathcal{I}$ is the identity operator (see [31] for more details). This formal definition resembles the practical procedure for computing the eigenvalues when the Bloch-Torrey operator $\mathcal{B}_g$ is represented by an infinite-dimensional matrix, which is then truncated and diagonalized numerically (see Appendix A). The eigenvalue branches $\lambda_j^{(g)}$ can merge and split at branch points but, apart from these points, they are smooth functions of $g$. We use this property to order the eigenvalues $\lambda_j^{(g)}$ according to the increasing order of Laplacian eigenvalues $\lambda_j^{(0)}$. In other words, one first orders the eigenvalues at $g = 0$ and then preserves their order by continuity of branches as $g$ increases. At each branch point, the order of merged eigenvalues is lost but they can be re-ordered in any convenient way. This ordering procedure does not
ensure an increasing order of Re{$\lambda_j^{(g)}$} for any g but it facilitates the visualization and interpretation of the spectrum. Most importantly, the associated eigenfunctions $v_j^{(g)}(x)$ also change smoothly with g and preserve their symmetries, except for branch points (see below).

### 3. Bloch-Torrey operator for a sphere

We consider restricted diffusion inside a sphere of radius R with reflecting boundary and apply the gradient along the z axis: $G = Ge_z$. The Bloch-Torrey operator can be written in spherical coordinates $(r, \theta, \phi)$ as

$$ \mathcal{B}_g^z = -\left(\partial_r^2 + \frac{2}{r}\partial_r + \frac{1}{r^2}\partial_\xi(1-\xi^2)\partial_\xi + \frac{1}{r^2(1-\xi^2)}\partial_\phi^2\right) + igr\xi, \quad (10) $$

where $\xi = \cos\theta$. Since the gradient operator does not depend on the azimuthal angle $\phi$ and the initial transverse magnetization is uniform, the considered problem is axisymmetric with respect to the z axis. In other words, the Bloch-Torrey equation does not change the uniformity with respect to $\phi$, i.e., the transverse magnetization remains independent of $\phi$. For this reason, one often considers the reduced Bloch-Torrey operator without the azimuthal part:

$$ \hat{\mathcal{B}}_g = -\left(\partial_r^2 + \frac{2}{r}\partial_r + \frac{1}{r^2}\partial_\xi(1-\xi^2)\partial_\xi\right) + igr\xi. \quad (11) $$

In fact, most former studies were focused in this operator and its matrix representation on the basis of the Laplace operator [6, 9, 12–14, 46, 47]. In turn, the full operator $\mathcal{B}_g^z$ in Eq. (10) is needed to deal with gradient pulses in different directions or with inhomogeneous initial magnetization. The related extension of the matrix formalism was introduced in [15]. In Appendix B, we recall the matrix elements for constructing the eigenvalues and eigenfunctions of both operators. Similarly, one can introduce the Bloch-Torrey operators $\mathcal{B}_g^x$ and $\mathcal{B}_g^y$ when the gradient is applied along x and y coordinates, respectively. Even though these two operators have different matrix representations (see Appendix B), the rotational invariance of the sphere ensures that the spectra of the three operators $\mathcal{B}_g^x$, $\mathcal{B}_g^y$ and $\mathcal{B}_g^z$ are identical. In turn, their eigenfunctions can be matched by an appropriate rotation of spherical coordinates (i.e., by choosing the spherical coordinates with the z axis aligned with the desired gradient direction). For this reason, we focus on the operator $\mathcal{B}_g^z$ in the following and compare its spectral properties to those of the reduced operator $\hat{\mathcal{B}}_g$.

When there is no gradient ($g = 0$), the eigenbasis of the (negative) Laplace operator $\mathcal{B}_0^z = -\Delta$ is fairly well known; in particular, the separation of variables yields the Laplacian eigenfunctions $u_{nkm}(r, \theta, \phi) \propto j_n(\alpha_{nk}r/R)P_n^m(\cos\theta)e^{im\phi}$, where $j_n(z)$ is the spherical Bessel function of the first kind, $P_n^m(z)$ is the associated Legendre polynomial, and $\alpha_{nk}$ are the positive zeros of the derivative $j'_n(z)$ ensuring the Neumann boundary condition. Here each Laplacian eigenfunction is parameterized by a triple index $nkm$ that reflects its symmetries, with $n = 0, 1, 2, ...$ being the order of $j_n(z)$, $k = 0, 1, 2, ...$ enumerating the zeros $\alpha_{nk}$, and $m = -n, -n+1, ..., n$. The associated eigenvalues
<table><caption>**Table 1.** First 17 eigenvalues $\lambda_{nkm} = \lambda_j^{(0)}$ of the (negative) Laplace operator, $B_0 = -\Delta$, in the unit sphere with reflecting boundary. The position $j$ of the eigenvalue $\lambda_j^{(0)}$ in the ordered sequence is used to enumerate the branch $\lambda_j^{(g)}$ for $g \neq 0$. Bold font highlights the indices of eigenfunctions that are axisymmetric (with $m = 0$).</caption><thead><tr><th>j</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th></tr></thead><tbody><tr><td>nkm</td><td>000</td><td>100</td><td>10(-1)</td><td>101</td><td>20(-1)</td><td>201</td><td>200</td><td>20(-2)</td><td>202</td></tr><tr><td>λ<sub>j</sub><sup>(0)</sup></td><td>0</td><td>4.33</td><td>4.33</td><td>4.33</td><td>11.17</td><td>11.17</td><td>11.17</td><td>11.17</td><td>11.17</td></tr><tr><td>j</td><td>10</td><td>11</td><td>12</td><td>13</td><td>14</td><td>15</td><td>16</td><td>17</td><td></td></tr><tr><td>nkm</td><td>010</td><td>300</td><td>30(-2)</td><td>302</td><td>30(-1)</td><td>301</td><td>30(-3)</td><td>303</td><td></td></tr><tr><td>λ<sub>j</sub><sup>(0)</sup></td><td>20.19</td><td>20.38</td><td>20.38</td><td>20.38</td><td>20.38</td><td>20.38</td><td>20.38</td><td>20.38</td><td></td></tr></tbody></table>

$\lambda_{nkm} = \alpha_{nk}^2/R^2$ do not depend on $m$ and are thus $(2n + 1)$ times degenerate. Writing these eigenvalues in an increasing order (see Table 1), we use the position $j$ of each eigenvalue in the sequence to enumerate the branches $\lambda_j^{(g)}$. Some ambiguities in the eigenvalue ordering procedure caused by the degeneracy of the Laplacian eigenvalues can be fixed manually.

### 3.1. Eigenvalues

Figure 1 summarizes the spectral properties of the Bloch-Torrey operator $B_g^z$ for the unit sphere ($R = 1$). One sees the real part of first 17 eigenvalues $\lambda_j^{(g)}$ as functions of $g$. At $g = 0$, one retrieves the Laplacian eigenvalues $\lambda_{nkm}$ with their degeneracies. For instance, three branches $\lambda_2^{(g)}$, $\lambda_3^{(g)}$ and $\lambda_4^{(g)}$ start from 4.33 at $g = 0$ but two of them coincide for all $g$, $\lambda_3^{(g)} \equiv \lambda_4^{(g)}$, resulting in a single upper curve. Similarly, five branches $\lambda_5^{(g)}, \dots, \lambda_9^{(g)}$ start from 11.17 at $g = 0$, but $\lambda_5^{(g)} \equiv \lambda_6^{(g)}$ for all $g$ result in a single lower curve, and $\lambda_8^{(g)} \equiv \lambda_9^{(g)}$ result in a single upper curve. These preserved degeneracies are related to the fact that the Bloch-Torrey operator $B_g^z$ does not affect the azimuthal angle. In fact, as $P_n^{-m}(x) = (-1)^m P_n^m(x)$, two Laplacian eigenfunctions corresponding to $+m$ and $-m$ exhibit the identical dependence on $r$ and $\theta$ and therefore remain indistinguishable even in the presence of the applied gradient along $z$ coordinate. As a consequence, the dependence $e^{im\phi}$ of the Laplacian eigenfunctions on the angle $\phi$ is preserved for the eigenfunctions of the Bloch-Torrey operator $B_g^z$. As the integral of the eigenfunctions containing the factor $e^{im\phi}$ with $m \neq 0$ over the sphere $\Omega$ vanishes, they do not contribute to the macroscopic signal. The related eigenvalues are shown by dashed lines. In turn, the eigenvalues shown by solid lines correspond to the eigenfunctions $v_j^{(g)}$ that inherited their independence of $\phi$ from the Laplacian eigenfunctions $u_{nk0}$ and thus do contribute to the signal. Expectedly, these eigenvalues could be directly obtained by diagonalizing the reduced Bloch-Torrey operator $\hat{B}_g$. In other words, the difference between the spectra of the operators $B_g^z$ and $\hat{B}_g$ is the presence of additional eigenvalues (shown by dashed lines) in the former case.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Real part of the first 17 eigenvalues $\lambda_j^{(g)}$ of the Bloch-Torrey operator $B_g^z$ for the unit sphere ($R=1$). Dashed lines indicate the eigenvalues that do not contribute to the macroscopic signal and thus do not appear in the spectrum of the reduced Bloch-Torrey operator $\hat{B}_g$. Colored snapshots show the $xz$ projection of the real part of the corresponding eigenfunction, evaluated at $g=0$ (on the left) and at $g=15$ (near vertical dashed line). Color indicates changes of $\text{Re}\{v_j^{(g)}\}$ from -1.5 (dark blue) to 1.5 (dark red), with the colorbar shown at right bottom, being the same for all snapshots. The values of $|m|$ determining the dependence $e^{im\phi}$ on the azimuthal angle $\phi$ are shown on the left. Four branch points are seen: $g_1^s \approx 5.622$, $g_2^s \approx 12.1$, $g_3^s \approx 20.1$, and $g_4^s \approx 23.84$. The eigenvalues and eigenfunctions were constructed via the matrix formalism (see Appendix A), in which the matrices were truncated at 333.</figcaption></figure>
The rotation invariance of the sphere implies the PT symmetry of the Bloch-Torrey operator [48, 49]. As a consequence, its eigenvalues are either real, or form complex-conjugate pairs (see [20, 31] for further discussions). This general property is confirmed on Fig. 1. Moreover, one can observe four branch (or exceptional) points $g_i^s$, at which real eigenvalues merge to become complex-conjugate pairs: $g_1^s \approx 5.622$, $g_2^s \approx 12.1$, $g_3^s \approx 20.1$, and $g_4^s \approx 23.84$. Note that the branch points $g_1^s$ and $g_4^s$ are of order 2 (i.e., two simple eigenvalues merge here), while the branch points $g_2^s$ and $g_3^s$ are of order 4 (two pairs of twice degenerate eigenvalues merge). To our knowledge, this is the first observation of a branch point of order 4 for the Bloch-Torrey operator (the previous studies [21, 31] revealed only branch points of order 2). As discussed earlier, twice degenerate eigenvalues correspond to the eigenfunctions that do not contribute to the signal. In particular, the reduced Bloch-Torrey operator $\hat{B}_g$ seems to possess only branch points of order 2. We also note that the preserved dependence of eigenfunctions on $\phi$ via $e^{im\phi}$ implies a simple branching rule: only the branches containing at $g=0$ the Laplacian eigenvalues $\lambda_{nkm}$ with the same $|m|$ can merge. For instance, the branches $\lambda_1^{(g)}$ and $\lambda_2^{(g)}$ corresponding to $\lambda_{000}$ and $\lambda_{100}$ (with $m=0$) merge at $g_1^s$; the branches $\lambda_3^{(g)}, \dots, \lambda_6^{(g)}$ corresponding to $\lambda_{10(-1)}, \lambda_{101}, \lambda_{20(-1)}, \lambda_{201}$ (with $|m|=1$) merge at $g_2^s$, and so on. Note that our numerical study did not reveal branch points of other orders except 2 and 4. We expect that their existence is unlikely but a mathematical proof of this statement remains an open problem. We also stress that branch points should be distinguished from “crossing” points, at which two (or more) eigenvalues cross, without changing their properties. For instance, the pair of real eigenvalues $\lambda_5^{(g)} \equiv \lambda_6^{(g)}$ crosses a single real eigenvalue $\lambda_7^{(g)}$ at $g \approx 9.3$. Three corresponding eigenfunctions form an orthogonal basis of the subspace of dimension 3. In contrast, one (or more) eigenfunction disappears at the branch point (see further discussion in [31]).

## 3.2. Eigenfunctions

Figure 1 also presents the $xz$ projections of the real part of the first 17 eigenfunctions of the Bloch-Torrey operator $B_g^z$ at $g=0$ and $g=15$. These snapshots help to visualize how the geometric structure of each eigenfunction changes with $g$. As discussed earlier, the Laplacian eigenfunctions $u_{nkm}$ and $u_{nk(-m)}$ exhibit the same dependence on $\phi$ and thus keep this property in the presence of the gradient along the $z$ axis, as confirmed by snapshots at $g=15$. Note that the $xz$ projection of some eigenfunctions is close to 0 (green color); in fact, such an eigenfunction should be orthogonal to its pair and thus exhibit most variations in other projections. One also sees how the symmetries of the first six eigenfunctions change after the branch point.

Let us inspect this change in more detail. Figure 2 illustrates the drastic change in the shape of the eigenfunctions $v_1^{(g)}$ and $v_2^{(g)}$ when $g$ crosses the branch point $g_1^s \approx 5.622$. We first consider the eigenfunction $v_1^{(g)}$ (bottom row). The uniform property of $v_1^{(0)}$ is immediately broken for any $g > 0$, as confirmed by the second panel showing $v_1^{(g)}$ at $g=1$. A similar geometric pattern of $v_1^{(g)}$ was observed for even small $g$ (not shown).
<figure><img src="image_4.png" /><figcaption>**Figure 2.** xz projection of the real part of the eigenfunctions $v_1^{(g)}$ (bottom row) and $v_2^{(g)}$ (top row) of the Bloch-Torrey operator $\mathcal{B}_g^z$ for the unit sphere for several values of g. The associated eigenvalue is indicated on the top of each plot. The branch point is $g_1^s \approx 5.622$. Note that color range changes between different panels.</figcaption></figure>

It is worth noting, however, that $v_1^{(1)}$ varies from 0.489 to 0.493 and thus remains very close to a constant (as $v_1^{(0)}$). As g increases up to $g_1^s$, the shape of $v_1^{(g)}$ remains visually unchanged but its variations grow rapidly. This is the consequence of the normalization by $\langle v_1^{(g)}, v_1^{(g)} \rangle^{-1/2}$. In fact, as g approaches the branch point $g_1^s$, $\langle v_1^{(g)}, v_1^{(g)} \rangle$ vanishes and thus the normalization factor diverges, as discussed in [31]. At $g = 5.63 > g_1^s$, the shape of the eigenfunction $v_1^{(g)}$ has drastically changed and started to exhibit variations along the z axis, as imposed by the applied gradient. Further increase of g does not change this symmetry but enhances the localization of the eigenfunction on the South pole.

A similar behavior is observed for the second eigenfunction $v_2^{(g)}$: its shape, inherited from the Laplacian eigenfunction $u_{100}$, is preserved for $g < g_1^s$ and then drastically changes to another shape exhibiting variations along z axis. Moreover, as the eigenvalues $\lambda_1^{(g)}$ and $\lambda_2^{(g)}$ form a complex conjugate pair for $g > g_1^s$, the associated eigenfunctions exhibit the reflection symmetry: $v_2^{(g)}(x) = [v_1^{(g)}(R_z x)]^*$, where $R_z$ is the reflection with respect to the xy plane (i.e., z is replaced by $-z$). Finally, a similar behavior is observed (but not shown here) for other eigenfunctions that drastically change their shapes at the branch point of their eigenvalues (e.g., compare the eigenfunctions $v_3^{(g)}, \dots, v_6^{(g)}$ shown in Fig. 1 at $g = 0$ and $g = 15$).

### 3.3. Macroscopic signal

When the duration $\delta$ of the gradient pulses is sufficiently long, only few eigenmodes with small $\text{Re}\{\lambda_j^{(g)}\}$ do contribute to the signal. The structure of the spectrum shown in Fig. 1 suggests to keep only the first two eigenmodes in the spectral expansion (5):

$$S \approx e^{-2D_0\delta\lambda_1^{(g)}} [C_{1,1}^{(g)} + 2\text{Re}\{C_{1,2}^{(g)}\}e^{-D_0\delta(\lambda_2^{(g)} - \lambda_1^{(g)})} + C_{2,2}^{(g)}e^{-2D_0\delta(\lambda_2^{(g)} - \lambda_1^{(g)})}], \quad (12)$$

with the coefficients $C_{j,j'}^{(g)}$ given by Eq. (6); note that we used the property $\mathcal{B}_{-g} = \mathcal{B}_g^*$ that implies $C_{j',j}^{(g)} = [C_{j,j'}^{(g)}]^*$. One can distinguish two scenarios according to whether the eigenvalues $\lambda_1^{(g)}$ and $\lambda_2^{(g)}$ are real or complex.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The coefficients $C_{1,1}^{(g)}$ and $\text{Re}\{C_{1,2}^{(g)}\}$ from Eq. (6) characterizing the contributions of the first two eigenfunctions to the signal for the unit sphere. Vertical dashed line indicates the branch point $g_1^s \approx 5.622$, at which both coefficients diverge.</figcaption></figure>

(i) When $0 < g < g_1^s$, the eigenvalues $\lambda_1^{(g)}$ and $\lambda_2^{(g)}$ are real and simple. If $D_0\delta(\lambda_2^{(g)} - \lambda_1^{(g)}) \gg 1$, the last two terms in Eq. (12) can be neglected, yielding the one-mode approximation for the signal,

$$S \approx S_{\text{one}} = C_{1,1}^{(g)} e^{-2D_0\delta\lambda_1^{(g)}}. \qquad (13)$$

Note that this approximation is not valid when $g$ is close to the branch point $g_1^s$.

(ii) When $g > g_1^s$, the eigenvalue $\lambda_1^{(g)}$ is complex and paired with $\lambda_2^{(g)} = [\lambda_1^{(g)}]^*$. As a consequence, one has $C_{2,2}^{(g)} = C_{1,1}^{(g)}$ so that Eq. (12) can be written as

$$S \approx S_{\text{two}} = 2e^{-2D_0\delta \text{Re}\{\lambda_1^{(g)}\}} [C_{1,1}^{(g)} + \text{Re}\{C_{1,2}^{(g)} e^{2iD_0\delta \text{Im}\{\lambda_1^{(g)}\}}\}]. \qquad (14)$$

Figure 3 shows the dependence of the coefficients $C_{1,1}^{(g)}$ and $\text{Re}\{C_{1,2}^{(g)}\}$ on $g$. When $g$ approaches the branch point $g_1^s$, the normalization of the involved eigenfunctions $v_1^{(g)}$ and $v_2^{(g)}$ diverges, resulting in the divergence of these coefficients: $C_{1,1}^{(g)} \to +\infty$ and $\text{Re}\{C_{1,2}^{(g)}\} \to -\infty$. However, as discussed in [31], these diverging contributions to the signal compensate each other and thus imply no resonant behavior of the signal near $g_1^s$. In other words, the signal changes smoothly with $g$ even at the branch point.

The accuracy of the approximations (13, 14) is illustrated on Fig. 4. Expectedly, both approximations fail at very small $\delta$ when many eigenfunctions are needed in the spectral expansion (5) to get the signal. In turn, both approximations become accurate at larger $\delta$. On panel (b), one can also notice oscillations due to the second term in Eq. (14). Their period is controlled by the imaginary part of $\lambda_1^{(g)}$. At high gradients, the leading term of the large-$g$ asymptotic expansion of $\lambda_1^{(g)}$ is $gR$ [17, 20, 26], so that the last factor in Eq. (14) is approximately $e^{2i\gamma GR\delta}$. One thus retrieves a diffusion-diffraction pattern [51–55], which is more common for short gradient pulses (see also a comparison between the localization regime and narrow-pulse approximation in [11]). Note also that
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Signal as a function of the gradient pulse duration $\delta$ for a sphere of radius $R = 10~\mu\text{m}$, with $\gamma = 2.675 \cdot 10^8~\text{rad/T/s}$ (protons), $D_0 = 2.3 \cdot 10^{-9}~\text{m}^2/\text{s}$ (water molecules), $G = 17~\text{mT/m}$ (a) and $G = 129~\text{mT/m}$ (b). Circles present the exact spectral expansion (5), truncated to 333 terms; dashed line indicates the exact matrix representation (4) with the same truncation; solid lines show one-mode and two-modes approximations (13, 14) for panels (a) and (b), respectively. Their parameters are: $R^2\lambda_1^{(g)} \approx 0.188$ and $C_{1,1} \approx 1.14$ for $G = 17~\text{mT/m}$ (or $g = 2$, panel (a)); and $R^2\lambda_1^{(g)} \approx 4.67 + 6.68i$, $C_{1,1} \approx 1.12$ and $C_{1,2} \approx -0.46 + 0.18i$ for $G = 129~\text{mT/m}$ (or $g = 15$, panel (b)).</figcaption></figure>

the next-order corrections to $\text{Im}(\lambda_1^{(g)})$ can significantly alter this behavior. In turn, the real part of the first eigenvalue behaves at large $g$ as [17, 20]:

$$ \text{Re}\{\lambda_1^{(g)}\} = \frac{|a'_1|}{2\ell_g^2} + \frac{1}{\sqrt{R\ell_g^{3/2}}} - \frac{\sqrt{3}}{4|a'_1|R\ell_g} + O(\ell_g^{-1/2}) \quad (g \to \infty), \quad (15) $$

where $\ell_g = (\gamma G/D_0)^{-1/3} = g^{-1/3}$, and $a'_1 \approx -1.02$ is the first zero of the derivative of the Airy function $\text{Ai}(z)$. Substituting the leading order of Eq. (15) into Eq. (14), one retrieves the stretched exponential decay of the signal: $\ln S \propto G^{2/3}$.

## 4. Bloch-Torrey operator for capped cylinders

In this section, we consider restricted diffusion in a capped cylinder of radius $R$ and height $H$: $\Omega = \{(x, y, z) \in \mathbb{R}^3 : x^2 + y^2 < R^2, -H/2 < z < H/2\}$. Breaking rotational invariance, this shape allows us to investigate how the domain anisotropy can affect the spectrum of the Bloch-Torrey operator for a gradient in an arbitrary direction. Since the capped cylinder is axisymmetric, there is no difference between $x$ and $y$ directions so that one can focus on gradients in the $xz$ plane, for instance, by setting $G = G(e_x \cos \eta + e_z \sin \eta)$, where $\eta$ is the angle with respect to the horizontal axis in the $xz$ plane. The corresponding Bloch-Torrey operator, denoted as $\mathcal{B}_g^{(\eta)}$, reads in the cylindrical coordinates $(r, \theta, z)$ as

$$ \mathcal{B}_g^{(\eta)} = -\Delta + ig(x \cos \eta + z \sin \eta) = \mathcal{B}_{g \cos \eta}^d + \mathcal{B}_{g \sin \eta}^i, \quad (16) $$

where

$$ \mathcal{B}_g^d = -\underbrace{\left(\partial_r^2 + \frac{1}{r}\partial_r + \frac{1}{r^2}\partial_\theta^2\right)}_{=\Delta_d} + igr \cos \theta, \quad \mathcal{B}_g^i = -\partial_z^2 + igz \quad (17) $$
<table><caption>**Table 2.** First 13 eigenvalues $\lambda_{nklm} = \lambda_j^{(0)}$ of the (negative) Laplace operator, $B_0 = -\Delta$, in the capped cylinder with $R = H = 1$. The position $j$ of the eigenvalue $\lambda_j^{(0)}$ in the ordered sequence is used to enumerate the branch $\lambda_j^{(g)}$ for $g \neq 0$. Twice degenerate eigenvalues (corresponding to $l = 1$ and $l = 2$) are shown together, e.g., $\lambda_2^{(0)} = \lambda_3^{(0)}$.</caption><thead><tr><th>j</th><th>1</th><th>2-3</th><th>4-5</th><th>6</th><th>7-8</th><th>9</th><th>10-11</th><th>12-13</th></tr></thead><tbody><tr><td>nklm</td><td>0010</td><td>10(1-2)0</td><td>20(1-2)0</td><td>0011</td><td>10(1-2)1</td><td>0110</td><td>30(1-2)0</td><td>20(1-2)1</td></tr><tr><td>λ<sub>j</sub><sup>(0)</sup></td><td>0</td><td>3.39</td><td>9.33</td><td>9.87</td><td>13.26</td><td>14.68</td><td>17.65</td><td>19.20</td></tr></tbody></table>

are the Bloch-Torrey operators in the disk and in the interval, respectively. As these two operators act on different variables, the eigenfunctions of $B_g^{(\eta)}$ are factored, while its eigenvalues are obtained as all possible sums of the eigenvalues of $B_g^i$ and $B_g^d$. The operator $B_g^i$, also known as the (complex) Airy operator [50], was thoroughly studied in [18, 21, 24], whereas $B_g^d$ was analyzed in [17, 20, 26, 31]. We aim at understanding how their spectral properties are superimposed in the case of a capped cylinder.

The matrix elements of the operator $B_g^{(\eta)}$ are derived in Appendix C. In particular, the Laplacian eigenfunctions,

$$u_{nklm}(r, \theta, z) \propto J_n(\alpha_{nk}r/R)s_l(n\theta) \cos(\pi m(z + H/2)/H),$$

are enumerated by multi-index $nkjm$, with $n = 0, 1, 2, \dots$ being the order of the Bessel function $J_n(z)$ of the first kind, $k = 0, 1, 2, \dots$ being the index of the zeros $\alpha_{nk}$ of $J'_n(z)$, $l$ distinguishing between $s_1(z) = \cos(z)$ and $s_2(z) = \sin(z)$, and $m = 0, 1, 2, \dots$ characterizing oscillations along $z$ axis. The eigenvalues $\lambda_{nklm} = \alpha_{nk}^2/R^2 + \pi^2 m^2/H^2$ do not depend on $l$ and are in general either simple (for $n = 0$) or twice degenerate (for $n > 0$), but higher degeneracies are possible. As previously, we use the ordered sequence of these eigenvalues to enumerate the eigenvalue branches $\lambda_j^{(g)}$ of the Bloch-Torrey operators (see Table 2). Throughout this section, we fix $R = H = 1$ and then explore the anisotropy by changing the gradient direction (angle $\eta$), as explained below.

### 4.1. Parallel and perpendicular directions of the gradient

We start with two simple cases when the gradient is either aligned with the cylinder axis and thus $B_g^{(\pi/2)} = B_g^z = B_0^d + B_g^i$, or lies in the transverse $xy$ plane so that $B_g^{(0)} = B_g^x = B_g^d + B_0^i$.

Figure 5 presents real parts of the first 13 eigenvalues of the Bloch-Torrey operators $B_g^z$ (panel (a)) and $B_g^x$ (panel (b)). Let us first inspect the spectrum of $B_g^z$, which is the sum of the (negative) Laplace operator $B_0^d = -\Delta_d$ in the disk and the Bloch-Torrey operator $B_g^i$ on the interval $(-H/2, H/2)$. The spectrum of $B_g^i$ was thoroughly investigated in [18, 21, 24]; for instance, Fig. 3 from [21] shows the real and imaginary parts of several eigenvalues. In particular, the branch of the first two eigenvalues $\lambda_1^{(g)}$ and $\lambda_2^{(g)}$ can be retrieved in Fig. 5(a), which zooms out Fig. 3 from [21] to a smaller
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Real part of the first 13 eigenvalues of the Bloch-Torrey operators $B_g^z$ (a) and $B_g^x$ (b) for a capped cylinder with $R = H = 1$. Dashed lines indicate the eigenvalues that do not contribute to the signal. Vertical lines indicate the positions of branch points: $g_1^i \approx 18.06$ (a) and $g_1^d \approx 3.76$, $g_2^d \approx 9.39$, $g_3^d \approx 13.87$ (b). Note that each pair of twice degenerate eigenvalues appears as a single branch on panel (a). Truncation order was 320.</figcaption></figure>

range of $g$. Moreover, Stoller et al. studied the branch points of $B_g^i$ and found an explicit formula [21], which reads in our notations (see also discussion in [31]):

$$g_k^i = \sqrt{3} \frac{27}{4} j_k^2, \quad \text{where} \quad J_{-2/3}(j_k) = 0 \quad (k = 1, 2, \dots). \qquad (18)$$

In particular, one gets $g_1^i \approx 18.06$ and $g_2^i \approx 229.35$. The position of the first branch point $g_1^i$, as indicated by the vertical line, is in excellent agreement with this prediction. The major difference between Fig. 5(a) and Fig. 3 from [21] is that the pair of eigenvalues $\lambda_1^{(g)}$ and $\lambda_2^{(g)}$ is replicated and shifted vertically by adding the eigenvalues of $-\Delta_d$, e.g., 3.3900, 9.3284, 14.6820, etc. As a consequence, there are infinitely many pairs of eigenvalues that branch at each value $g_1^i, g_2^i, \dots$. Note that the “shifted” eigenvalues correspond to the eigenfunctions $u_{nkl}^d$ of $-\Delta_d$ that are orthogonal to $u_{001}^d = \text{const}$, so that the resulting eigenfunctions $v_j^{(g)}$ do not contribute to the signal.

Let us now look at the spectrum of the Bloch-Torrey operator $B_g^x$, which is the
sum of $\mathcal{B}_g^d$ and the second derivative $\mathcal{B}_0^i = -\partial_z^2$ on the interval. The spectrum of the former operator is simply replicated and shifted vertically by the eigenvalues $\pi^2m^2/H^2$ ($m=1,2,\dots$) of $-\partial_z^2$. These shifted eigenvalues are shown by dashed lines in Fig. 5(b) because the associated eigenfunctions do not contribute to the signal due to the presence of the factor $\cos(\pi m(z+H/2)/H)$, whose integral over the interval $(-H/2, H/2)$ vanishes for any $m > 0$. In turn, the eigenvalues shown by solid lines correspond to $m = 0$ and repeat the spectrum of the Bloch-Torrey operator in the disk (compare with Fig. 9 from [31]). In particular, one observes here three branch points at $g_1^d \approx 3.76$, $g_2^d \approx 9.39$, and $g_3^d \approx 13.87$, which are replicated along the vertical axis by adding $\pi^2m^2/H^2$.

## 4.2. Changing gradient direction

The structure of the spectra for both considered operators $\mathcal{B}_g^z$ and $\mathcal{B}_g^x$ was rather simple because one of two terms in Eq. (16) was independent of $g$ and thus just shifted vertically the spectrum of the other. For intermediate angles $\eta$, both terms in Eq. (16) depend on $g$, and the angle $\eta$ controls rescaling of each spectrum through the factors $g \cos \eta$ and $g \sin \eta$. Changing $\eta$, one can “tune” continuously the spectra of $\mathcal{B}_{g \cos \eta}^d$ and $\mathcal{B}_{g \sin \eta}^i$, and see how their features change. Most importantly, even though the eigenfunctions are still factored along the longitudinal and transverse directions $e_z$ and $e_x$, the nonzero gradients along $e_x$ and $e_z$ directions break the symmetries of these factors so that all eigenfunctions may contribute to the signal.

Figure 6 shows the real parts of the first eigenvalues of the Bloch-Torrey operator $\mathcal{B}_g^{(\eta)}$ for two angles: $\eta = \pi/4$ and $\eta = \pi/3$. One sees how the two spectra of $\mathcal{B}_{g \cos \eta}^d$ and $\mathcal{B}_{g \sin \eta}^i$ are superimposed and tuned by the angle. In particular, as $g$ is multiplied by $\cos \eta$ or $\sin \eta$, the branch points of the spectra are rescaled by $1/\cos \eta$ and $1/\sin \eta$, respectively. Changing the gradient angle, one can continuously shift the positions of branch points and thus re-organize the branch structure of the spectrum.

Let us now focus on the particular setting, shown in Fig. 7, in which the angle $\eta$ is tuned to make equal the rescaled branch points $g_1^d / \cos \eta$ and $g_1^i / \sin \eta$, i.e., $\eta = \tan^{-1}(g_1^d/g_1^i) \approx 1.3661 \approx 78.3^\circ$. The position $g_1^c = g_1^d / \cos \eta \approx 18.5$ of the first group of branch points is indicated by a vertical line. As compared to Fig. 6, there are more branches that merge at a single branch point. For instance, there are four eigenvalues $\lambda_1^{(g)}, \lambda_2^{(g)}, \lambda_6^{(g)}, \lambda_7^{(g)}$ that merge at the first branch point, two eigenvalues $\lambda_3^{(g)}$ and $\lambda_8^{(g)}$ that merge at the second one, and four eigenvalues $\lambda_4^{(g)}, \lambda_5^{(g)}, \lambda_{12}^{(g)}, \lambda_{13}^{(g)}$ that merge at the third one. Curiously, these last four eigenvalues form two *distinct* branches for $g > g_1^c$: one pair $\lambda_4^{(g)}, \lambda_{12}^{(g)}$ and the other pair $\lambda_5^{(g)}, \lambda_{13}^{(g)}$ (one can notice a small deviation between them at $g = 30$, which is further increased at larger $g$). It turns out that the eigenvalues in each group keep the same value of the index $l$ that distinguished $s_1(n\theta) = \cos(n\theta)$ and $s_2(n\theta) = \sin(n\theta)$ in the eigenfunctions of the Laplace operator. According to Table 2, one sees that the eigenvalues with indices $j = 1, 2, 6, 7$ have $l = 1$, those with $j = 3, 8$ have $l = 2$, those with $j = 4, 12$ have $l = 1$, and those with $j = 5, 13$ have $l = 2$, etc. This is expected because the gradient is applied in the $xz$ plane and thus preserves the
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Real part of the first 13 eigenvalues of the Bloch-Torrey operator $B_g^{(\eta)}$ for a capped cylinder with $R = H = 1$, and the gradient applied in the $xz$ plane at the angle $\eta$ with respect to the $x$-axis: $\eta = \pi/4$ (**a**) and $\eta = \pi/3$ (**b**). Vertical dashed lines show the positions of rescaled branch points $g_1^d \approx 3.76/\cos(\eta)$, $g_2^d \approx 9.36/\cos(\eta)$, $g_3^d \approx 13.87/\cos(\eta)$, and $g_1^i \approx 18.06/\sin(\eta)$ associated to the operators $B_g^d$ and $B_g^i$. Note that each pair of twice degenerate eigenvalues appears as a single branch on panel (b). Truncation order was 320.</figcaption></figure>

distinction between $l = 1$ and $l = 2$ in the angular dependence (in the similar way as the gradient applied along $z$ axis preserved the dependence on the angle $\phi$ for a sphere, as discussed above). The insets of Fig. 7 illustrate the real parts of the $xz$ projections of the associated eigenfunctions $v_j^{(g)}$ at $g = 0$ and $g = 20$. One sees how the symmetries of eigenfunctions drastically change at the branch point.

This is further illustrated on Fig. 8, showing the real part of the $xz$ projection of four eigenfunctions $v_1^{(g)}$, $v_2^{(g)}$, $v_4^{(g)}$ and $v_6^{(g)}$ at different $g$. As earlier for the case of a sphere, the constant eigenfunction $v_1^{(0)}$ is rapidly destroyed by the gradient; in turn, the symmetries of other three eigenfunctions $v_2^{(0)}$, $v_4^{(0)}$ and $v_6^{(0)}$ are still visible (though slightly perturbed) at $g = 15$, which is below the branch point 18.5. When $g$ exceeds the branch point, the symmetries change, and the eigenfunctions start to be more and more localized. At first thought, the observed localization pattern is puzzling. In fact,
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Real part of the first 16 eigenvalues $\lambda_j^{(g)}$ of the Bloch-Torrey operator $B_g^{(\eta)}$ for the capped cylinder with $R = H = 1$, and the gradient applied in the $xz$ plane at the angle $\eta = \tan^{-1}(18.06/3.76) \approx 1.3661 \approx 78.3^\circ$ with respect to $x$ axis. Colored snapshots show the $xz$ projection of the real part of the corresponding eigenfunction, evaluated at $g = 0$ (on the left) and at $g = 20$ (near vertical dashed line). Color indicates changes of $\text{Re}\{v_j^{(g)}\}$ from -1.5 (dark blue) to 1.5 (dark red), with the colorbar shown at right bottom, being the same for all snapshots. The branch point is located at 18.5. The eigenvalues and eigenfunctions were constructed via the matrix formalism (see Appendix A), in which the matrices were truncated at 320.</figcaption></figure>
for a smooth boundary, the localization occurs at specific boundary points $x_b$, at which the normal vector $n(x_b)$ is parallel to the gradient $G$ [26]. In other words, the gradient direction determines the location of the localized eigenfunctions on the boundary. As the boundary of the capped cylinder is not smooth, the asymptotic analysis from [26] is not applicable. Moreover, Fig. 8 shows that the above selection rule is actually not valid here. In fact, the gradient is directed at the angle $\eta \approx 78.3^\circ$ in the $xz$ plane with respect to the $x$ axis. One might thus expect localization at left bottom and right top corners. However, the eigenfunctions $v_1^{(g)}$ and $v_2^{(g)}$ at $g = 100$ are localized at the other corners, namely, the left up and the right bottom corners, respectively. Moreover, $v_4^{(g)}$ is localized in the middle of the upper edge. For a capped cylinder, this behavior can be explained by the factored structure of eigenfunctions. As discussed earlier, every eigenfunction $v_j^{(g)}$ is the product of an eigenfunction of $B_g^d$ and an eigenfunction of $B_g^i$ for the disk and the interval, respectively. At high enough $g$, both factors are localized: the eigenfunction of $B_g^d$ is localized at either of two opposite points of the disk along the $x$ axis, while the eigenfunction of $B_g^i$ is localized at either of two endpoints of the interval $(-H/2, H/2)$ along the $z$ axis. The eigenfunction $v_j^{(g)}$ can thus be localized in or near any corner of the $xz$ projection. Further investigation of the localization in domains with nonsmooth boundaries presents an open mathematical problem.

### 4.3. Exploring the anisotropy

The above illustrations of the spectral properties of the Bloch-Torrey operator were realized for a particular capped cylinder with $H = R = 1$. While it is easy to replicate the above results for any capped cylinder, the overall structure of the spectrum does not change significantly. In fact, one can rewrite the Bloch-Torrey operator in the capped cylinder of radius $R$ and height $H$ as

$$B_g = \frac{1}{R^2}(-\bar{\Delta}_d + ig \cos \eta R^3 \bar{r} \cos \theta) + \frac{1}{H^2}(-\partial_{\bar{z}}^2 + ig \sin \eta H^3 \bar{z}), \quad (19)$$

where bar denotes rescaled quantities: $\bar{r} = r/R$, $\bar{z} = z/H$, $\bar{\Delta}_d = R^2 \Delta_d$. Setting $\bar{g} = g\sqrt{R^6 \cos^2 \eta + H^6 \sin^2 \eta}$ and introducing the angle $\bar{\eta}$ such that $\tan(\bar{\eta}) = (H^3/R^3) \tan(\eta)$, one has

$$B_g = \frac{1}{R^2} \bar{B}_{\bar{g} \cos \bar{\eta}}^d + \frac{1}{H^2} \bar{B}_{\bar{g} \sin \bar{\eta}}^i, \quad (20)$$

where the first term is the Bloch-Torrey operator in the unit disk (divided by $R^2$) with the dimensionless gradient $\bar{g} \cos \bar{\eta}$, whereas the second term is the Bloch-Torrey operator in the unit interval (divided by $H^2$) with the dimensionless gradient $\bar{g} \sin \bar{\eta}$. In the previous subsection, we used $R = H = 1$ and thus considered the sum of these two basic operators. In general, the prescribed lengths $R$ and $H$ re-define the gradient amplitude (from $g$ to $\bar{g}$) and the gradient angle (from $\eta$ to $\bar{\eta}$), as well as weighting factors $1/R^2$ and $1/H^2$ in the linear combination (20). The structure of the eigenfunctions of $B_g$ is thus affected by anisotropy of the capped cylinder only through $\bar{g}$ and $\bar{\eta}$; in turn, the spectrum is also controlled by the weights $1/R^2$ and $1/H^2$ that can reshape eigenvalue branches
<figure><img src="image_10.png" /><figcaption>**Figure 8.** $xz$ projection of the real part of four eigenfunctions $v_j^{(g)}$ of the Bloch-Torrey operator $B_g^{(\eta)}$ for a capped cylinder with $R = H = 1$. Different columns correspond to $j \in \{1, 2, 4, 6\}$, while different rows correspond to $g \in \{0, 15, 20, 100\}$. The gradient is applied in the $xz$ plane with the angle $\eta = \tan^{-1}(18.06/3.76) \approx 1.3661 \approx 78.3^\circ$ with respect to the $x$-axis. The associated eigenvalue is indicated on the top of each plot. The branch point is located at 18.5. Truncation order was 320.</figcaption></figure>

and shift the branch points. Developing experimental protocols that are sensitive to the shape of eigenfunctions will potentially allow to exploit this property in order to probe microscopic anisotropy of porous media at high gradients.

**5. Conclusion**

In this paper, we investigated the spectral properties of the Bloch-Torrey operator $B_g = -\Delta + igx$ in two three-dimensional domains: a sphere and a capped cylinder. These shapes are typical models in diffusion MRI, representing, e.g., a soma and a neuron in the brain tissue. While the general asymptotic behavior of eigenvalues and eigenfunctions was known in the limits of small and large $g$, there is no spectral result for the intermediate range of $g$, which is the most relevant for applications. In particular, we studied the structure of the spectrum, the dependence of eigenvalue branches on $g$, the existence of branch points, and drastic symmetry changes of eigenfunctions at branch points. Despite the geometric simplicity of the considered domains, we had to rely on the matrix formalism to construct the eigenvalues and eigenfunctions of $B_g$ numerically. We illustrated how different eigenvalue branches merge at branch points and how the symmetries of eigenfunctions, inherited from the Laplace operator, are destroyed at
these points. For a capped cylinder, we also showed the effect of anisotropy on the spectrum, in particular, how rotating the gradient direction allows one to rescale the spectra of the operators $B_g^d$ and $B_g^i$ in the orthogonal directions and thus to tune the branch points. The localization of eigenfunctions was shown to occur near the points with $z = \pm H/2$ and $r = R$, at which the cylindrical wall joins the top and bottom caps. As the boundary of a capped cylinder is not smooth at these points, the asymptotic behavior established in [26] is not applicable, and further analysis of localization in domains with nonsmooth boundaries is needed.

The present work lays the theoretical ground for various applications in diffusion MRI. As shown earlier, the macroscopic signal is getting more sensitive to the microstructure at high gradients [18, 56]. As a consequence, high-gradient diffusion MRI is a promising research direction with potential applications in material sciences, neurosciences and medicine [57–59]. An intuitive explanation of this enhanced sensitivity is that the localized eigenfunctions may probe selected boundary regions and thus access refined information on the microstructure. For instance, the coefficients $C_{j,j'}^{(g)}$ from Eq. (6) of the spectral expansion (5) are sensitive to the overlap between two eigenfunctions and thus may, potentially, probe a sort of spatial correlations between different boundary regions, at which these eigenfunctions are localized. Moreover, the access to the Bloch-Torrey operators $B_g^x, B_g^y$ and $B_g^z$ for three orthogonal directions allows one to analyze double-pulsed field-gradient experiments [15, 60–64] in terms of spectral expansions similar to Eq. (5). In particular, as the coefficients of such expansions involve different eigenfunctions of these non-commuting Bloch-Torrey operators, one may potentially reveal additional information on the microstructure such as its local anisotropy or curvature. More generally, an experimental exploration of eigenfunctions symmetry changes at branch points presents a very interesting but challenging task, and the developed spectral approach may pave a way towards new imaging modalities at high gradients.

## Data availability statement

No new data were created or analysed in this study.

## Acknowledgments

The author acknowledges the Alexander von Humboldt Foundation for support within a Bessel Prize award.

## Appendix A. Computation of spectral properties

In this Appendix, we extend the description of the numerical procedure from [31] that we use for computing the eigenvalues and eigenfunctions of the Bloch-Torrey operator $B_g$. It is inspired from the matrix formalism [6, 12, 46], in which the magnetization
is decomposed onto the complete basis of Laplacian eigenfunctions $u_k$ with Neumann boundary condition, which are known explicitly for simple domains (e.g., a disk and a sphere). Throughout this Appendix, we do not discuss mathematical aspects of the problem such the convergence of spectral representations. Our goal here is to provide a practical recipe for numerical computations. In order to deal with dimensionless quantities, we will rescale lengths by the “size” $R$ of the confining domain $\Omega$, e.g., by its (half-)diameter (this choice does not matter in practice). For the examples considered in the paper, $R$ is the radius of the sphere or of the capped cylinder.

We search an eigenfunction $v_j^{(g)}$ of the Bloch-Torrey operator as

$$v_j^{(g)}(\boldsymbol{x}) = \sum_k X_{j,k}^{(g)} u_k(\boldsymbol{x}), \qquad (A.1)$$

with unknown coefficients $X_{j,k}^{(g)}$. Substituting Eq. (A.1) into the eigenvalue problem (7), one gets

$$\lambda_j^{(g)} \sum_k X_{j,k}^{(g)} u_k(\boldsymbol{x}) = \lambda_j^{(g)} v_j^{(g)}(\boldsymbol{x}) = (-\Delta + igx) v_j^{(g)}(\boldsymbol{x}) = \sum_k X_{j,k}^{(g)} (-\Delta + igx) u_k(\boldsymbol{x}).$$

Multiplying this equation by $u_{j'}^*(\boldsymbol{x})R^2$, integrating over $\Omega$, and using orthogonality of Laplacian eigenfunctions, we get for any $j'$

$$R^2 \lambda_j^{(g)} X_{j,j'}^{(g)} = \sum_k X_{j,k}^{(g)} (\Lambda_{k,j'} + i\bar{g}B_{k,j'}), \qquad (A.2)$$

where $\Lambda_{k,j'} = \delta_{k,j'}\lambda_k R^2$, $B_{k,j'} = \int_{\Omega} d\boldsymbol{x} \, u_k(\boldsymbol{x}) (x/R) u_{j'}^*(\boldsymbol{x})$, $\lambda_k$ are the eigenvalues of the (negative) Laplace operator $-\Delta$, and we introduced the dimensionless parameter $\bar{g} = gR^3 = R^3\gamma G/D_0$. The multiplication by $R^2$ ensured that both matrices $\Lambda$ and $B$ are dimensionless. In a matrix form, one has

$$\Lambda^{(g)} X = X(\Lambda + i\bar{g}B), \qquad (A.3)$$

where $\Lambda^{(g)}$ is the diagonal matrix of eigenvalues $R^2\lambda_n^{(g)}$ of the Bloch-Torrey operator $B_g$. As a consequence, the diagonalization of the matrix $\Lambda + i\bar{g}B$ yields $\Lambda^{(g)}$ and $\tilde{X}$, where $\tilde{X}$ is the matrix whose columns contain left eigenvectors, from which $X$ is obtained by complex-conjugate transpose: $X = \tilde{X}^{\dagger,*}$, where $\dagger$ denotes transpose without complex conjugation, i.e., $[X^{\dagger}]_{k,j} = X_{j,k}$. For instance, one could use the matlab commands

[V,LambdaG,Xtilde] = eig(Lambda + 1i*gbar*B); X = Xtilde';

to get $\Lambda^{(g)}$ and $X$.

To ensure the normalization (9) of eigenfunctions, one can use the representation (A.1) that implies

$$\int_{\Omega} d\boldsymbol{x} \, v_j^{(g)}(\boldsymbol{x}) v_{j'}^{(g)}(\boldsymbol{x}) = [XWX^{\dagger}]_{j,j'} = \delta_{j,j'}, \qquad (A.4)$$

where

$$W_{k,k'} = \int_{\Omega} d\boldsymbol{x} \, u_k(\boldsymbol{x}) u_{k'}(\boldsymbol{x}). \qquad (A.5)$$
The eigenfunctions {$u_k$} of the Laplace operator with Neumann boundary condition can be chosen to be real-valued, in which case $W$ is the identity matrix. However, it may also be convenient to employ complex-valued Laplacian eigenfunctions. In this more general setting, the scalar product in $L_2(\Omega)$ includes complex-conjugation; as a consequence, even though two eigenfunctions $u_k$ and $u_{k'}$ are orthogonal to each other, $(u_k, u_{k'})_{L_2(\Omega)} = \int_{\Omega} u_k^* u_{k'} = 0$, their integral in Eq. (A.5) may not be zero, and the matrix $W$ is not necessarily identity. For instance, we used this convention in Sec. 3 by choosing $u_{nkm} \propto e^{im\phi}$ that yields some nonzero off-diagonal elements of the matrix $W$ in Eq. (B.12) for a sphere. We stress that this purely conventional issue does not affect any spectral property of the Bloch-Torrey operator.

We also recall that the integral in Eq. (A.4) may be zero for $j = j'$ at specific values of $g$. Moreover, if an eigenvalue $\lambda_j^{(g)}$ is degenerate, the associated eigenfunctions form an eigenspace so that the coefficients $X_{j,k}^{(g)}$ are not defined uniquely but up to a rotation in that eigenspace (see Appendix D for a simple orthogonalization procedure). We stress that this ambiguity does not affect the resulting macroscopic signal but may render the interpretation of eigenfunctions more sophisticated.

According to Eq. (6), the coefficients $C_{j,j'}^{(g)}$ from the spectral expansion (5) of the signal can be written as

$$C_{j,j'}^{(g)} = \mu_j^{(-g)} \Gamma_{j,j'}^{(g)} \mu_{j'}^{(g)}, \qquad (A.6)$$

where

$$\mu_j^{(g)} = \frac{1}{\sqrt{|\Omega|}} \int_{\Omega} d\boldsymbol{x} \, v_j^{(g)}(\boldsymbol{x}) \qquad (A.7)$$

is the projection of the eigenfunction $v_j^{(g)}(\boldsymbol{x})$ onto a constant, and

$$\Gamma_{j,j'}^{(g)} = \int_{\Omega} d\boldsymbol{x} \, v_j^{(-g)}(\boldsymbol{x}) v_{j'}^{(g)}(\boldsymbol{x}) \qquad (A.8)$$

is the overlap between two eigenfunctions $v_j^{(-g)}(\boldsymbol{x})$ and $v_{j'}^{(g)}(\boldsymbol{x})$. The integrals in Eqs. (A.7, A.8) can be computed directly by using the representation (A.1):

$$\mu_j = \frac{1}{\sqrt{|\Omega|}} \int_{\Omega} d\boldsymbol{x} \sum_k X_{j,k} u_k(\boldsymbol{x}) = [XU]_j, \qquad (A.9)$$

where

$$U_k = \frac{1}{\sqrt{|\Omega|}} \int_{\Omega} d\boldsymbol{x} \, u_k(\boldsymbol{x}). \qquad (A.10)$$

Since $u_k(\boldsymbol{x})$ are orthogonal to $u_0(\boldsymbol{x}) = 1/\sqrt{|\Omega|}$ for Neumann boundary condition, one gets

$$\mu_j = X_{j,0}. \qquad (A.11)$$

Similarly,

$$\Gamma_{j,j'} = \sum_k X_{j,k}^* \sum_{k'} X_{j',k'} \int_{\Omega} d\boldsymbol{x} \, u_k^*(\boldsymbol{x}) u_{k'}(\boldsymbol{x}) = [X^* X^\dagger]_{j,j'}. \qquad (A.12)$$
In practice, one can only construct a finite-dimensional approximation of the infinite-dimensional matrix $\Lambda + i\bar{g}B$ by using a large but finite number $N$ of eigenmodes of the Laplace operator (see Appendix B and Appendix C for details). The numerical diagonalization of the truncated matrix of size $N \times N$ yields $N$ eigenvalues, which are *expected* to converge to the eigenvalues of the Bloch-Torrey operator as $N$ goes to infinity. This conjecture was supported by numerical evidence: when computing a given number of eigenvalues (and eigenfunctions) by diagonalizing truncated matrices with larger and larger $N$, we observed that they rapidly become almost independent of $N$. Another indirect evidence for this convergence comes from the fact that the eigenvalues obtained from truncated matrices obey the large-$g$ asymptotic behavior derived for the eigenvalues of the Bloch-Torrey operator [20, 26]. At the same time, a rigorous proof of the convergence is still missing. In fact, for many non-Hermitian matrices, the eigenvalues are known to be very sensitive to perturbations (such as truncation) so that the convergence may fail. For example, the spectra of banded Toeplitz matrices of increasing sizes do not converge to the spectrum of their limiting operator acting on an appropriate infinite-dimensional space [65–67]. A systematic study of the convergence presents thus an interesting perspective.

## Appendix B. Matrix elements for a sphere

We summarize the matrix elements needed for computing the eigenfunctions of the Bloch-Torrey operator in a sphere of radius $R$ with reflecting boundary. For the reduced operator $\hat{B}_g$, the matrix representation was derived in [6, 46]:

$$ \Lambda_{nk,n'k'} = \delta_{n,n'} \delta_{k,k'} \frac{\alpha_{nk}^2}{R^2}, \qquad (B.1) $$

where $\alpha_{nk}$ are the positive zeros of $j'_n(z)$ (with $n = 0, 1, 2, ...$), enumerated by $k = 0, 1, 2, ...$, and

$$ B_{nk,n'k'} = \delta_{n,n' \pm 1} \frac{n+n'+1}{(2n+1)(2n'+1)} \beta_{nk} \beta_{n'k'} \frac{\alpha_{nk}^2 + \alpha_{n'k}^2 - n(n'+1) - n'(n+1) + 1}{(\alpha_{nk}^2 - \alpha_{n'k'}^2)^2}, \quad (B.2) $$

with

$$ \beta_{nk} = \left( \frac{(2n+1)\alpha_{nk}^2}{\alpha_{nk}^2 - n(n+1)} \right)^{1/2}, \quad \beta_{00} = \sqrt{3/2}. \qquad (B.3) $$

Here we use the double index $nk$ to enumerate the elements of the matrices $\Lambda$ and $B$.

The matrix elements of the full Bloch-Torrey operator $B_g$ were obtained in [15] that we reproduce below for completeness. As discussed in Sec. 3, the Laplacian eigenfunctions $u_{nkm}$ are now enumerated by a triple index $nkm$, with $m$ ranging from $-n$ to $n$, while the eigenvalues $\lambda_{nkm} = \alpha_{nk}^2/R^2$ do not depend on $m$ and thus $(2n+1)$ times degenerate. As a consequence, the matrix $\Lambda$ takes a block-diagonal form, with the elements

$$ \Lambda_{nkm,n'k'm'} = \delta_{m,m'} \Lambda_{nk,n'k'} = \delta_{m,m'} \delta_{n,n'} \delta_{k,k'} \frac{\alpha_{nk}^2}{R^2}. \qquad (B.4) $$
In turn, the matrix $B$ representing the gradient term, depends on the direction $e_G$ of the gradient $G = Ge_G$. Encoding this direction in spherical coordinates by angles $\theta_G$ and $\phi_G$ as

$$e_G = \sin \theta_G \cos \phi_G e_x + \sin \theta_G \sin \phi_G e_y + \cos \theta_G e_z, \quad (B.5)$$

one can represent this gradient by the matrix

$$B_{nkm,n'k'm'} = \frac{1}{R} \int_{\Omega} d\mathbf{x} [u_{nkm}(\mathbf{x})]^* (e_G \cdot \mathbf{x}) u_{n'k'm'}(\mathbf{x}) \quad (B.6) \\ = \left[ \sin \theta_G \cos \phi_G B^x + \sin \theta_G \sin \phi_G B^y + \cos \theta_G B^z \right]_{nkm,n'k'm'}$$

with three matrices $B^x$, $B^y$ and $B^z$, representing respectively the operators of multiplication by $x$, $y$, and $z$ in the Laplacian eigenbasis.

For the gradient along $z$ axis, one gets

$$B_{nkm,n'k'm'}^z = \delta_{m,m'} B_{nk,n'k'} \sqrt{1 - \frac{m^2}{(\max\{n, n'\})^2}} \quad (|m| \le \min\{n, n'\}). \quad (B.7)$$

Expectedly, one retrieves the matrix elements $B_{nk,n'k'}$ when $m = m' = 0$. For two other components, the only nonzero elements are

$$B_{nkm,(n+1)k'm'}^x = \frac{B_{nk,(n+1)k'}}{2} \left( \delta_{m',m-1} \frac{\sqrt{(n-m+1)(n-m+2)}}{n+1} - \delta_{m',m+1} \frac{\sqrt{(n+m+1)(n+m+2)}}{n+1} \right), \quad (B.8)$$

$$B_{nkm,(n-1)k'm'}^x = -\frac{B_{nk,(n-1)k'}}{2} \left( \delta_{m',m-1} \frac{\sqrt{(n+m-1)(n+m)}}{n} - \delta_{m',m+1} \frac{\sqrt{(n-m-1)(n-m)}}{n} \right), \quad (B.9)$$

and

$$B_{nkm,(n+1)k'm'}^y = i \frac{B_{nk,(n+1)k'}}{2} \left( \delta_{m',m-1} \frac{\sqrt{(n-m+1)(n-m+2)}}{n+1} + \delta_{m',m+1} \frac{\sqrt{(n+m+1)(n+m+2)}}{n+1} \right), \quad (B.10)$$

$$B_{nkm,(n-1)k'm'}^y = -i \frac{B_{nk,(n-1)k'}}{2} \left( \delta_{m',m-1} \frac{\sqrt{(n+m-1)(n+m)}}{n} + \delta_{m',m+1} \frac{\sqrt{(n-m-1)(n-m)}}{n} \right). \quad (B.11)$$

We also compute the elements of the matrix $W$ defined by Eq. (A.5):

$$W_{nkm,n'k'm'} = (-1)^m \delta_{n,n'} \delta_{k,k'} \delta_{m,-m'}. \quad (B.12)$$

We stress that this matrix is not the identity.
## Appendix C. Matrix elements for a capped cylinder

In order to study the effect of anisotropy, we consider diffusion in a capped cylinder of radius $R$ and height $H$: $\Omega = \{x = (x, y, z) \in \mathbb{R}^3 : x^2 + y^2 < R^2, -H/2 < z < H/2\}$. Since the lateral diffusion along the $z$ axis is independent from the transverse diffusion in the $xy$ plane, one usually considers separately the gradient encoding in these orthogonal directions. For a standard pulsed-gradient spin-echo sequence with two opposite gradient pulses, it is therefore enough to consider two reduced Bloch-Torrey operators: $\hat{B}_g^z = -\partial_z^2 + igz$ on the interval $(-H/2, H/2)$, and $\hat{B}_g^{xy} = -(\partial_r^2 + r^{-1}\partial_r) + igr \cos \theta$ for a disk of radius $R$. The matrix elements for both operators were given explicitly in [6, 46]. The spectral properties of $\hat{B}_g^z$ were thoroughly investigated in [18, 21, 24], while the spectrum of $\hat{B}_g^{xy}$ was discussed in [17, 31] (see references therein). However, more sophisticated pulsed-gradient sequences with several gradient directions require the knowledge of the whole Bloch-Torrey operator $B_g$. We summarize the matrix elements needed for constructing the spectrum of this operator.

The separation of variables allows one to get the eigenbasis of the Laplacian operator explicitly in cylindrical coordinates $(r, \theta, z)$ as

$$u_{nklm}(r, \theta, z) = u_{nkl}^{d}(r, \theta) \frac{\sqrt{2 - \delta_{m,0}}}{\sqrt{H}} \cos(\pi m(z + H/2)/H), \quad (C.1)$$

where $u_{nkl}^d(r, \theta)$ are the Laplacian eigenfunctions for a disk of radius $R$:

$$u_{nkl}^{d}(r, \theta) = \frac{\sqrt{2 - \delta_{n,0}}}{\sqrt{\pi R}} \frac{\beta_{nk}}{J_{n}(\alpha_{nk})} J_{n}(\alpha_{nk}r/R) \times \begin{cases} \cos(n\theta) & (l=1), \\ \sin(n\theta) & (l=2), \end{cases} \quad (C.2)$$

where $J_n(z)$ is the Bessel function of the first kind, $\alpha_{nk}$ are the positive zeros of $J'_n(z)$ enumerated by $k = 0, 1, 2, \dots$, and

$$\beta_{nk} = \frac{\alpha_{nk}}{\sqrt{\alpha_{nk}^2 - n^2}}, \quad \beta_{00} = 1. \quad (C.3)$$

One sees that the Laplacian eigenfunctions $u_{nklm}$ are enumerated by the multi-index $nklm$, with $m = 0, 1, 2, \dots$. The associated eigenvalue is simply

$$\lambda_{nklm} = \frac{\alpha_{nk}^2}{R^2} + \frac{\pi^2 m^2}{H^2}. \quad (C.4)$$

In general, the eigenvalue $\lambda_{nklm}$ is twice degenerate for $n > 0$ and simple for $n = 0$ (in this case, $u_{0k2}^d(r, \theta) \equiv 0$ is not an eigenfunction and thus excluded). However, one can get higher-order degeneracy for specific values of the aspect ratio $H/R$.

The structure of the Laplacian eigenfunctions allows one to construct explicitly the matrices $\Lambda$ and $B^i$ ($i = x, y, z$) representing the Laplace operator and the gradient along three coordinate axes. From the practical point of view, it is convenient to construct these matrices by reproducing their block structure:

$$\Lambda = \begin{pmatrix} \Lambda_d & 0 & 0 & 0 & \dots \\ 0 & \Lambda_d + (\pi^2/H^2)\mathbb{I} & 0 & 0 & \dots \\ 0 & 0 & \Lambda_d + (4\pi^2/H^2)\mathbb{I} & 0 & \dots \\ 0 & 0 & 0 & \Lambda_d + (9\pi^2/H^2)\mathbb{I} & \dots \\ \dots & \dots & \dots & \dots & \dots \end{pmatrix}, \quad (C.5)$$
$$B^{x,y} = \begin{pmatrix} B_d^{x,y} & 0 & 0 & 0 & \dots \\ 0 & B_d^{x,y} & 0 & 0 & \dots \\ 0 & 0 & B_d^{x,y} & 0 & \dots \\ 0 & 0 & 0 & B_d^{x,y} & \dots \\ \dots & \dots & \dots & \dots & \dots \end{pmatrix}, \qquad (C.6)$$

$$B^z = \begin{pmatrix} B_{0,0}^i \mathbb{I} & B_{0,1}^i \mathbb{I} & B_{0,2}^i \mathbb{I} & B_{0,3}^i \mathbb{I} & \dots \\ B_{1,0}^i \mathbb{I} & B_{1,1}^i \mathbb{I} & B_{1,2}^i \mathbb{I} & B_{1,3}^i \mathbb{I} & \dots \\ B_{2,0}^i \mathbb{I} & B_{2,1}^i \mathbb{I} & B_{2,2}^i \mathbb{I} & B_{2,3}^i \mathbb{I} & \dots \\ B_{3,0}^i \mathbb{I} & B_{3,1}^i \mathbb{I} & B_{3,2}^i \mathbb{I} & B_{3,3}^i \mathbb{I} & \dots \\ \dots & \dots & \dots & \dots & \dots \end{pmatrix}, \qquad (C.7)$$

where $\mathbb{I}$ is the identity matrix, $B_{m,m'}^i$ are the matrix elements for the interval:

$$B_{m,m'}^i = ((-1)^{m+m'} - 1) \sqrt{2-\delta_{m,0}} \sqrt{2-\delta_{m',0}} \frac{m^2+m'^2}{\pi^2(m^2-m'^2)^2} \quad (m \neq m'), \qquad (C.8)$$

and $B_{m,m}^i = 0$, with $m, m' = 0, 1, 2, \dots$. In turn, $\Lambda_d$ and $B_d^{x,y}$ are the matrices representing the Laplace operator and the gradient for the disk. As for the case of a sphere, there matrices were first derived explicitly in [6, 46] for the reduced Bloch-Torrey operator and then extended in [15]. We re-derive the extended expressions in a slightly different form. Skipping straightforward computations, we get

$$[\Lambda_d]_{nkl,n'k'l'} = \delta_{n,n'} \delta_{k,k'} \delta_{l,l'} \frac{\alpha_{nk}^2}{R^2} \qquad (C.9)$$

and

$$[B_d^x]_{nk1,n'k'1} = B_{nk,n'k'}^d, \quad [B_d^x]_{nk2,n'k'2} = B_{nk,n'k'}^d(1 - \delta_{n+n',1}), \qquad (C.10)$$

$$[B_d^x]_{nk1,n'k'2} = [B_d^x]_{nk2,n'k'1} = 0, \qquad (C.11)$$

where

$$B_{nk,n'k'}^d = \delta_{n,n' \pm 1} (1 + \delta_{n,0} + \delta_{n',0})^{1/2} \beta_{nk} \beta_{n'k'} \frac{\alpha_{nk}^2 + \alpha_{n'k'}^2 - 2nn'}{(\alpha_{nk}^2 - \alpha_{n'k'}^2)^2} \qquad (C.12)$$

is the matrix $B^d$ for the reduced Bloch-Torrey operator derived in [6, 46]. Similarly, one has

$$[B_d^y]_{nk1,n'k'1} = [B_d^y]_{nk2,n'k'2} = 0, \qquad (C.13)$$

$$[B_d^y]_{nk1,(n+1)k'2} = B_{nk,(n+1)k'}^d, \quad [B_d^y]_{nk1,(n-1)k'2} = -B_{nk,(n-1)k'}^d(1 - \delta_{n+n',1}), \qquad (C.14)$$

and

$$[B_d^y]_{nk2,(n-1)k'1} = B_{nk,(n-1)k'}^d, \quad [B_d^y]_{nk2,(n+1)k'1} = -B_{nk,(n+1)k'}^d(1 - \delta_{n+n',1}). \qquad (C.15)$$
## Appendix D. Orthogonalization of eigenfunctions with degenerate eigenvalues

The relation (8) ensures the orthogonality of eigenfunctions $v_j^{(g)}$ and $v_{j'}^{(g)}$ with respect to the bilinear form $\langle \cdot, \cdot \rangle$ if two associated eigenvalues are distinct. In turn, if an eigenvalue is $n$ times degenerate, i.e., there are distinct indices $j_1, \dots, j_n$ such that $\lambda_{j_1}^{(g)} = \lambda_{j_2}^{(g)} = \dots = \lambda_{j_n}^{(g)}$, the associated eigenfunctions $v_{j_i}^{(g)}$ form an eigenspace of dimension $n$, in which any $n$ linearly independent combinations of $v_{j_i}^{(g)}$ can be chosen as eigenfunctions. Even though this ambiguity does not affect the computation of the macroscopic signal via the matrix formalism, the coefficients in the spectral expansion (5) can be sensitive to this choice. Moreover, a proper graphical representation of each eigenfunction and its visual interpretation require to choose the linear combinations that respect the orthogonality. In this Appendix, we briefly describe a straightforward orthogonalization procedure for the case $n = 2$. This procedure was sufficient for the cases of a sphere and a capped cylinder. We also focus on the generic setting when $g$ is not a branch point (indeed, as one eigenfunction disappears at the branch point, the analysis of this particular situation is more subtle, see discussion in [31]).

Let $v_j$ and $v_{j'}$ be two eigenfunctions with the same eigenvalue (we dropped here the superscript (g) for brevity), and their (non)-orthogonality is characterized by the matrix

$$C = \begin{pmatrix} \langle v_j, v_j \rangle & \langle v_j, v_{j'} \rangle \\ \langle v_{j'}, v_j \rangle & \langle v_{j'}, v_{j'} \rangle \end{pmatrix}. \qquad (D.1)$$

We aim at constructing two linear combinations,

$$\hat{v}_j = av_j + bv_{j'}, \qquad \hat{v}_{j'} = cv_j + dv_{j'}, \qquad (D.2)$$

whose unknown coefficients $a, b, c, d$ are chosen to ensure the orthonormality of these combinations: $\langle \hat{v}_j, \hat{v}_{j'} \rangle = \delta_{j,j'}$. In a matrix form, we have

$$\begin{pmatrix} \hat{v}_j \\ \hat{v}_{j'} \end{pmatrix} = \begin{pmatrix} a & b \\ c & d \end{pmatrix} \begin{pmatrix} v_j \\ v_{j'} \end{pmatrix}, \qquad (D.3)$$

so that

$$\begin{pmatrix} \langle \hat{v}_j, \hat{v}_j \rangle & \langle \hat{v}_j, \hat{v}_{j'} \rangle \\ \langle \hat{v}_{j'}, \hat{v}_j \rangle & \langle \hat{v}_{j'}, \hat{v}_{j'} \rangle \end{pmatrix} = \begin{pmatrix} a & b \\ c & d \end{pmatrix} C \begin{pmatrix} a & c \\ b & d \end{pmatrix}. \qquad (D.4)$$

By equating the left-hand side to the identity matrix, one can multiply this equation on the left by $(a \ b)^{-1}$ and on the right by $(a \ c)^{-1}$, to get equations on the unknown coefficients $a, b, c, d$:

$$\frac{1}{(ad - bc)^2} \begin{pmatrix} b^2 + d^2 & -ab - cd \\ -ab - cd & a^2 + c^2 \end{pmatrix} = C. \qquad (D.5)$$

To proceed, we parameterize the unknown coefficients as:

$$a = A^{-1} \cos \alpha, \quad b = A^{-1} \sin \alpha, \quad c = -B^{-1} \sin \alpha, \quad d = B^{-1} \cos \alpha. \qquad (D.6)$$
In the Hermitian setting, $\alpha$ could be interpreted as a rotation angle, while A and B as rescaling factors. In our case, this is a formal representation; in particular, all three parameters $\alpha$, A, B can take complex values. Substituting these expressions into Eq. (D.5) yields three equations:

$$A^2 \cos^2 \alpha + B^2 \sin^2 \alpha = C_{1,1}, \qquad (D.7)$$

$$A^2 \sin^2 \alpha + B^2 \cos^2 \alpha = C_{2,2}, \qquad (D.8)$$

$$(A^2 - B^2) \sin(2\alpha) = 2C_{1,2}. \qquad (D.9)$$

One can solve these equations as

$$\alpha = \frac{1}{2}\text{atan}\left(\frac{2C_{1,2}}{C_{1,1} - C_{2,2}}\right), \qquad (D.10)$$

$$A^2 = \frac{C_{1,1} \cos^2 \alpha - C_{2,2} \sin^2 \alpha}{\cos^2 \alpha - \sin^2 \alpha}, \qquad B^2 = \frac{C_{2,2} \cos^2 \alpha - C_{1,1} \sin^2 \alpha}{\cos^2 \alpha - \sin^2 \alpha}. \qquad (D.11)$$

When $C_{1,1}$ is close to $C_{2,2}$, $\alpha$ is close to $\pm\pi/4$, so that the above expressions may be numerically unstable. In this case, it is more convenient to use another representation:

$$A^2 = \frac{C_{1,1} + C_{2,2}}{2} + \frac{C_{1,2}}{\sin(2\alpha)}, \qquad B^2 = \frac{C_{1,1} + C_{2,2}}{2} - \frac{C_{1,2}}{\sin(2\alpha)}. \qquad (D.12)$$

In this way, we have explicit expressions for the linear transformation from a non-orthogonal pair of eigenfunctions $v_j$ and $v_{j'}$ to an orthonormal pair of eigenfunctions $\hat{v}_j$ and $\hat{v}_{j'}$.

In practice, once the matrix X of coefficients in Eq. (A.1) is found by solving the eigenvalue problem (A.3), one can evaluate the matrix $XWX^\dagger$ that represents the orthogonality of the eigenfunctions $v_j^{(g)}$. According to Eq. (A.4), this matrix should be equal to the identity matrix. As discussed above, this is ensured by Eq. (8) for any pair of eigenfunctions with distinct eigenvalues. As a consequence, nonzero off-diagonal elements of the matrix $XWX^\dagger$ are only possible for pairs of eigenfunctions with the same eigenvalue. One can therefore search for such nonzero off-diagonal elements and apply the above orthonormalization procedure for each such pair. This procedure was applied for most spectral computations in this work.

[1] Callaghan PT 1993 *Principles of Nuclear Magnetic Resonance Microscopy* (Clarendon Press, Oxford)

[2] Price WS 2009 *NMR Studies of Translational Motion: Principles and Applications* (Cambridge University Press, Cambridge)

[3] Le Bihan D and Johansen-Berg H 2012 *Diffusion MRI* at 25: Exploring brain tissue structure and function *NeuroImage* **61** 324-341

[4] Novikov D, Fieremans E, Jespersen S, and Kiselev VG 2018 *Quantifying brain microstructure with diffusion MRI: Theory and parameter estimation* *NMR Biomed* e3998

[5] Axelrod S and Sen PN 2001 *Nuclear magnetic resonance spin echoes for restricted diffusion in an inhomogeneous field: Methods and asymptotic regimes* *J. Chem. Phys.* **114**, 6878

[6] Grebenkov DS 2007 *NMR Survey of Reflected Brownian Motion* *Rev. Mod. Phys.* **79** 1077-1137

[7] Kiselev VG 2017 *Fundamentals of diffusion MRI physics* *NMR Biomed.* **30** e3602

[8] Torrey HC 1956 *Bloch Equations with Diffusion Terms* *Phys. Rev.* **104** 563
[9] Grebenkov DS 2010 Pulsed-gradient spin-echo monitoring of restricted diffusion in multilayered structures *J. Magn. Reson.* **205** 181-195

[10] Nguyen DV, Li J-R, Grebenkov DS, and Le Bihan D 2014 A finite elements method to solve the Bloch-Torrey equation applied to diffusion magnetic resonance imaging *J. Comput. Phys.* **263** 283-302

[11] Moutal N and Grebenkov DS 2020 The localization regime in a nutshell *J. Magn. Reson.* **320** 106836

[12] Barzykin AV 1998 Exact Solution of the Torrey-Bloch Equation for a Spin Echo in Restricted Geometries *Phys. Rev. B* **58** 14171

[13] Barzykin AV 1999 Theory of Spin Echo in Restricted Geometries under a Step-wise Gradient Pulse Sequence *J. Magn. Reson.* **139** 342

[14] Grebenkov DS 2008 Laplacian Eigenfunctions in NMR I. A Numerical Tool *Conc. Magn. Reson.* **32A** 277-301

[15] Özarslan E, Shemesh N, and Basser PJ 2009 A general framework to quantify the effect of restricted diffusion on the NMR signal with applications to double pulsed field gradient NMR experiments *J. Chem. Phys.* **130** 104702

[16] Tanner JE and Stejskal EO 1968 Restricted Self-Diffusion of Protons in Colloidal Systems by the Pulsed-Gradient, Spin-Echo Method *J. Chem. Phys.* **49** 1768-1777

[17] de Swiet TM, and Sen PN 1994 Decay of nuclear magnetization by bounded diffusion in a constant field gradient *J. Chem. Phys.* **100** 5597

[18] Grebenkov DS 2014 Exploring diffusion across permeable barriers at high gradients. II. Localization regime *J. Magn. Reson.* **248** 164-176

[19] Herberthson M, Özarslan E, Knutsson H, Westin C-F 2017 Dynamics of local magnetization in the eigenbasis of the Bloch-Torrey operator *J. Chem. Phys.* **146** 124201

[20] Moutal N, Demberg K, Grebenkov DS, and Kuder TA 2019 Localization regime in diffusion NMR: theory and experiments *J. Magn. Reson.* **305** 162-174

[21] Stoller SD, Happer W, and Dyson FJ 1991 Transverse spin relaxation in inhomogeneous magnetic fields *Phys. Rev. A* **44** 7459

[22] Robertson B 1966 Spin-echo decay of spins diffusion in a bounded region *Phys. Rev.* **151** 273

[23] Neuman CH 1974 Spin echo of spins diffusion in a bounded medium *J. Chem. Phys.* **60** 4508

[24] Grebenkov DS, Helffer B, and Henry R 2017 The complex Airy operator on the line with a semi-permeable barrier *SIAM J. Math. Anal.* **49** 1844-1894

[25] Moutal N and Grebenkov DS 2019 Diffusion across semi-permeable barriers: spectral properties, efficient computation, and applications *J. Scient. Comp.* **81** 1630-1654

[26] Grebenkov DS and Helffer B 2018 On the spectral properties of the Bloch-Torrey operator in two dimensions *SIAM J. Math. Anal.* **50** 622-676

[27] Almog Y, Grebenkov DS, and Helffer B 2018 Spectral semi-classical analysis of a complex Schrödinger operator in exterior domains *J. Math. Phys.* **59** 041501

[28] Almog Y, Grebenkov DS, and Helffer B 2019 On a Schrödinger operator with a purely imaginary potential in the semiclassical limit *Comm. Part. Diff. Eq.* **44** 1542-1604

[29] Moutal N, Moutal A, and Grebenkov DS 2020 Diffusion NMR in periodic media: efficient computation and spectral properties *J. Phys. A: Math. Theor.* **53** 325201

[30] Grebenkov DS, Helffer B, and Moutal N 2021 On the spectral properties of the Bloch-Torrey equation in infinite periodically perforated domains, Chapter 10 in "Partial Differential Equations, Spectral Theory, and Mathematical Physics: The Ari Laptev Anniversary Volume", Eds. P. Exner, R. L. Frank, F. Gesztesy, H. Holden and T. Weidl, EMS Series of Congress Reports Vol. 18 (EMS Press, Berlin) pp. 177-195.

[31] Moutal N and Grebenkov DS 2022 Spectral branch points of the Bloch-Torrey operator *J. Phys. A: Math. Theor.* **55** 455201

[32] Berry MV 2004 Physics of Nonhermitian Degeneracies *Czech. J. Phys.* **54** 1039

[33] Heiss WD 2004 Exceptional points of non-Hermitian operators *J. Phys. A: Math. Gen.* **37** 2455-
2464

[34] Seyranian AP, Kirillov ON and Mailybaev AA 2005 Coupling of eigenvalues of complex matrices at diabolic and exceptional points *J. Phys. A: Math. Gen.* **38** 1723-1740

[35] Kirillov ON, Mailybaev AA and Seyranian AP 2005 Unfolding of eigenvalue surfaces near a diabolic point due to a complex perturbation *J. Phys. A: Math. Gen.* **38** 5531-5546

[36] Rubinstein J, Sternberg P, and Ma Q 2007 Bifurcation Diagram and Pattern Formation of Phase Slip Centers in Superconducting Wires Driven with Electric Currents *Phys. Rev. Lett.* **99** 167003

[37] Cartarius H, Main J, and Wunner G 2007 Exceptional Points in Atomic Spectra *Phys. Rev. Lett.* **99** 173003

[38] Cejnar P, Heinze S, and Macek M 2007 Coulomb Analogy for Non-Hermitian Degeneracies near Quantum Phase Transitions *Phys. Rev. Lett.* **99** 100601

[39] Klaiman S, Günther U, and Moiseyev N 2008 Visualization of Branch Points in PT-Symmetric Waveguides *Phys. Rev. Lett.* **101** 080402

[40] Chang C-H, Wang S-M, and Hong T-M 2009 Origin of branch points in the spectrum of PT-symmetric periodic potentials *Phys. Rev. A* **80** 042105

[41] Ceci S, Döring M, Hanhart C, Krewald S, Meissner U-G, and Svarc A 2011 Relevance of complex branch points for partial wave analysis *Phys. Rev. C* **84** 015205

[42] Shapiro B and Zarembo K 2017 On level crossing in random matrix pencils. I. Random perturbation of a fixed matrix, *J. Phys. A: Math. Theor.* **50** 045201

[43] Grøsfjeld T, Shapiro B, and Zarembo K 2019 On level crossing in random matrix pencils. II. Random perturbation of a random matrix *J. Phys. A: Math. Theor.* **52** 214001

[44] Hürlimann MD, Helmer KG, de Swiet TM, Sen PN, and Sotak CH 1995 Spin Echoes in a Constant Gradient and in the Presence of Simple Restriction *J. Magn. Reson. A* **113** 260

[45] Williamson NH, Ravin R, Benjamini D, Merkle H, Falgairolle M, O'Donovan MJ, Blivis D, Ide D, Cai TX, Ghorashi NS, Bai R, Basser PJ 2019 Magnetic resonance measurements of cellular and sub-cellular membrane structures in live and fixed neural tissue *eLife* **8** e51101

[46] Grebenkov DS 2008 Analytical solution for restricted diffusion in circular and spherical layers under inhomogeneous magnetic fields *J. Chem. Phys.* **128** 134702

[47] Callaghan PT 1997 A Simple Matrix Formalism for Spin Echo Analysis of Restricted Diffusion under Generalized Gradient Waveforms *J. Magn. Reson.* **129** 74

[48] Moiseyev N 2011 *Non-Hermitian quantum mechanics* (Cambridge University Press)

[49] El-Ganainy R, Makris KG, Khajavikhan M, Musslimani ZH, Rotter S, and Christodoulides DN 2018 Non-Hermitian physics and PT symmetry *Nat. Phys.* **14** 11-19

[50] Helffer B 2013 *Spectral theory and its applications* (Cambridge University Press)

[51] Callaghan PT, Coy A, MacGowan D, Packer KJ, and Zelaya FO 1991 Diffraction-Like Effects in NMR Diffusion Studies of Fluids in Porous Solids *Nature* **351** 467-469

[52] Cotts RM 1991 Diffusion and diffraction *Nature* **351** 443-444

[53] Sen PN, Hürlimann MD, and de Swiet TM 1995 Debye-Porod law of diffraction for diffusion in porous media *Phys. Rev. B* **51** 601

[54] Gibbs SJ 1997 Observations of diffusive diffraction in a cylindrical pore by PFG NMR *J. Magn. Reson.* **124** 223-226

[55] Özarslan E and Basser PJ 2007 MR diffusion-"diffraction" phenomenon in multi-pulse-field-gradient experiments *J. Magn. Reson.* **188** 285-294

[56] Grebenkov DS 2018 Diffusion MRI/NMR at high gradients: challenges and perspectives *Micro. Meso. Mater* **269** 79-82

[57] Wedeen VJ, Rosene DL, Wang R, Dai G, Mortazavi F, Hagmann P, Kaas JH, and Tseng W-YI 2012 The Geometric Structure of the Brain Fiber Pathways *Science* **335** 1628-1634

[58] Huang SY et al. 2021 Connectome 2.0: Developing the next-generation ultra-high gradient strength human MRI scanner for bridging studies of the micro-, meso- and macro-connectome *NeuroImage* **243** 118530

[59] Williamson NH, Witherspoon VJ, Cai TX, Ravin R, Horkay F, and Basser PJ 2023 Low-field,
high-gradient NMR shows diffusion contrast consistent with localization or motional averaging
of water near surfaces *Magn. Reson. Lett.* **3** 90-107

[60] Callaghan PT and Komlosh ME 2002 Locally anisotropic motion in a macroscopically isotropic system: displacement correlation measured using double pulsed gradient spin-echo NMR *Magn. Reson. Chem.* **40** S15-S19

[61] Komlosh ME, Horkay F, Freidlin RZ, Nevo U, Assaf Y, and Basser PJ 2007 Detection of microscopic anisotropy in gray matter and in a novel tissue phantom using double Pulsed Gradient Spin Echo MR *J. Magn. Reson.* **189** 38-45

[62] Özarslan E and Basser PJ 2008 Microscopic anisotropy revealed by NMR double pulsed field gradient experiments with arbitrary timing parameters *J. Chem. Phys.* **128** 154511

[63] Özarslan E 2009 Compartment shape anisotropy (CSA) revealed by double pulsed field gradient MR *J. Magn. Reson.* **199** 56-67

[64] Jespersen SN, Lundell H, Sønderby CK, and Dyrby TB 2013 Orientationally invariant metrics of apparent compartment eccentricity from double pulsed field gradient diffusion experiments *NMR Biomed.* **26** 1647-1662

[65] Reichel L and Trefethen LN 1992 Eigenvalues and pseudo-eigenvalues of Toeplitz matrices *Linear Algebra Appl.* **162-164** 153-185

[66] Trefethen LN 1997 Pseudospectra of linear operators *SIAM Rev.* **39** 383-406

[67] Trefethen LN and Embree M 2005 *Spectra and Pseudospectra: The Behavior of Nonnormal Matrices and Operators* (Princeton University Press, New Jersey)