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

Denis S. Grebenkov  Laboratoire de Physique de la Matière Condensée,
CNRS – Ecole Polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau,
France<EMAIL_ADDRESS>
###### 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.

##  1\. Introduction

Diffusion magnetic resonance imaging (dMRI) is a non-invasive technique with
multiple applications in medicine, neurosciences and material sciences [1, 2,
3, 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^{\circ}$ 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, 6, 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(\bm{x},t)$ of the nuclei in a confining domain
$\Omega\subset\mathbb{R}^{d}$:

(1) $\partial_{t}m(\bm{x},t)=D_{0}\Delta
m(\bm{x},t)-i\gamma(\bm{G}(t)\cdot\bm{x})m(\bm{x},t)\quad(\bm{x}\in\Omega),$

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 $\bm{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(\bm{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(\bm{x},t)\bigr{|}_{\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
$\bm{x}$ is too small to be measured, only its integral over the confining
domain (or a voxel) is accessible in experiments:

(2) $S=\int\limits_{\Omega}d\bm{x}\,m(\bm{x},t).$

This macroscopic signal that can be accessed as a function of the gradient
profile $\bm{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
$\bm{G}(t)=\bm{G}$. Denoting by $x$ the coordinate axis in the direction of
the gradient, one has $(\bm{G}\cdot\bm{x})=Gx$, where $G=|\bm{G}|$ is the
gradient amplitude, and $x$ is the projection of $\bm{x}$ onto the direction
of $\bm{G}$. Defining the Bloch-Torrey operator as

(3) $\mathcal{B}_{g}=-\Delta+igx\quad(g=\gamma G/D_{0}),$

one can formally solve the Bloch-Torrey equation as
$m(\bm{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, 13, 14, 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

(4)
$S=\int\limits_{\Omega}d\bm{x}\,\biggl{(}e^{-D_{0}\delta\mathcal{B}_{-g}}e^{-D_{0}\delta\mathcal{B}_{g}}\frac{1}{|\Omega|}\biggr{)},$

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,\ldots$) to represent the above signal as [17, 18, 19, 20]

(5)
$S=\sum\limits_{j,j^{\prime}=1}^{\infty}C_{j,j^{\prime}}^{(g)}\,e^{-D_{0}\delta(\lambda_{j}^{(-g)}+\lambda_{j^{\prime}}^{(g)})},$

where the coefficients

(6)
$C_{j,j^{\prime}}^{(g)}=\frac{1}{|\Omega|}\left(\int\limits_{\Omega}d\bm{x}\,v_{j}^{(-g)}(\bm{x})\right)\left(\int\limits_{\Omega}d\bm{x}\,v_{j}^{(-g)}(\bm{x})\,v_{j^{\prime}}^{(g)}(\bm{x})\right)\left(\int\limits_{\Omega}d\bm{x}\,v_{j^{\prime}}^{(g)}(\bm{x})\right)$

characterize the overlap between two eigenfunctions $v_{j}^{(-g)}$ and
$v_{j^{\prime}}^{(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\mathrm{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
$\mathcal{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 [22, 23, 6]. 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, 33, 34, 35, 36, 37, 38, 39, 40, 41,
42, 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 is biological samples was recently 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 $\mathcal{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 $\mathcal{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 $\mathcal{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 (A), the matrix elements for a sphere (B) and
for a capped cylinder (C), as well as a simple orthogonalization procedure for
eigenfunctions with degenerate eigenvalues (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\geq 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)}(\bm{x})$
satisfying

(7)
$\mathcal{B}_{g}\,v_{j}^{(g)}(\bm{x})=\lambda_{j}^{(g)}\,v_{j}^{(g)}(\bm{x})\quad(\bm{x}\in\Omega),\qquad\partial_{n}v_{j}^{(g)}(\bm{x})=0\quad(\bm{x}\in\partial\Omega).$

As the Bloch-Torrey operator is not Hermitian for $g>0$, the standard scalar
product in $L_{2}(\Omega)$,
$(u,v)=\int\nolimits_{\Omega}d\bm{x}\,u^{*}(\bm{x})v(\bm{x})$, is replaced by
a bilinear form $\langle
u,v\rangle=\int\nolimits_{\Omega}d\bm{x}\,u(\bm{x})v(\bm{x})$. In particular,
the eigenfunctions $\\{v_{j}^{(g)}\\}$ are in general not orthogonal to each
other, $(v_{j}^{(g)},v_{j^{\prime}}^{(g)})\neq 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

(8) $\bigl{(}\lambda_{j}^{(g)}-\lambda_{j^{\prime}}^{(g)}\bigr{)}\langle
v_{j}^{(g)},v_{j^{\prime}}^{(g)}\rangle=0,$

so that if the eigenvalues $\lambda_{j}^{(g)}$ and
$\lambda_{j^{\prime}}^{(g)}$ are not equal, then $\langle
v_{j}^{(g)},v_{j^{\prime}}^{(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\nolimits_{\Omega}d\bm{x}\,[v_{j}^{(g)}]^{2}=0$.
In general, however, this integral is not zero, and we normalize the
eigenfunctions to have

(9) $\langle v_{j}^{(g)},v_{j}^{(g)}\rangle=1.$

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 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 $\mathrm{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)}(\bm{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: $\bm{G}=G\bm{e}_{z}$. The
Bloch-Torrey operator can be written in spherical coordinates
$(r,\theta,\phi)$ as

(10)
$\mathcal{B}_{g}^{z}=-\biggl{(}\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}\biggr{)}+igr\xi,$

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:

(11)
$\hat{\mathcal{B}}_{g}=-\biggl{(}\partial_{r}^{2}+\frac{2}{r}\partial_{r}+\frac{1}{r^{2}}\partial_{\xi}(1-\xi^{2})\partial_{\xi}\biggr{)}+igr\xi.$

In fact, most former studies were focused in this operator and its matrix
representation on the basis of the Laplace operator [47, 12, 13, 6, 46, 14,
9]. 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 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 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^{\prime}_{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,\ldots$ being the order of $j_{n}(z)$,
$k=0,1,2,\ldots$ enumerating the zeros $\alpha_{nk}$, and
$m=-n,-n+1,\ldots,n$. The associated eigenvalues
$\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.

$\displaystyle\begin{array}[]{c| c| c c c| c c c c c}j&{\bf 1}&{\bf
2}&3&4&5&6&{\bf 7}&8&9\\\ \hline\cr
nkm&000&100&10(-1)&101&20(-1)&201&200&20(-2)&202\\\
\lambda_{j}^{(0)}&0&4.33&4.33&4.33&11.17&11.17&11.17&11.17&11.17\\\
\end{array}$ $\displaystyle\begin{array}[]{c| c| c c c c c c c}j&{\bf 10}&{\bf
11}&12&13&14&15&16&17\\\ \hline\cr
nkm&010&300&30(-2)&302&30(-1)&301&30(-3)&303\\\
\lambda_{j}^{(0)}&20.19&20.38&20.38&20.38&20.38&20.38&20.38&20.38\\\
\end{array}$ Table 1.  First 17 eigenvalues $\lambda_{nkm}=\lambda_{j}^{(0)}$
of the (negative) Laplace operator, $\mathcal{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$). Figure 1.  Real part of the
first 17 eigenvalues $\lambda_{j}^{(g)}$ of the Bloch-Torrey operator
$\mathcal{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{\mathcal{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
$\mathrm{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}^{\rm
s}\approx 5.622$, $g_{2}^{\rm s}\approx 12.1$, $g_{3}^{\rm s}\approx 20.1$,
and $g_{4}^{\rm s}\approx 23.84$. The eigenvalues and eigenfunctions were
constructed via the matrix formalism (see A), in which the matrices were
truncated at 333.

###  3.1. Eigenvalues

Figure 1 summarizes the spectral properties of the Bloch-Torrey operator
$\mathcal{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)},\ldots,\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 $\mathcal{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
$\mathcal{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{\mathcal{B}}_{g}$. In other words, the
difference between the spectra of the operators $\mathcal{B}_{g}^{z}$ and
$\hat{\mathcal{B}}_{g}$ is the presence of additional eigenvalues (shown by
dashed lines) in the former case.

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}^{\rm s}$, at which real eigenvalues merge to
become complex-conjugate pairs: $g_{1}^{\rm s}\approx 5.622$, $g_{2}^{\rm
s}\approx 12.1$, $g_{3}^{\rm s}\approx 20.1$, and $g_{4}^{\rm s}\approx
23.84$. Note that the branch points $g_{1}^{\rm s}$ and $g_{4}^{\rm s}$ are of
order 2 (i.e., two simple eigenvalues merge here), while the branch points
$g_{2}^{\rm s}$ and $g_{3}^{\rm 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{\mathcal{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}^{\rm s}$; the branches $\lambda_{3}^{(g)},\ldots\lambda_{6}^{(g)}$
corresponding to $\lambda_{10(-1)}$, $\lambda_{101}$, $\lambda_{20(-1)}$,
$\lambda_{201}$ (with $|m|=1$) merge at $g_{2}^{\rm 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 $\mathcal{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 $z$,
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}^{\rm 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). 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}^{\rm 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}^{\rm 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}^{\rm 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}^{\rm 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}^{\rm s}$,
the associated eigenfunctions exhibit a reflection symmetry:
$v_{2}^{(g)}(\bm{x})=[v_{1}^{(g)}(R_{z}\bm{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)},\ldots,v_{6}^{(g)}$ shown in Fig. 1 at $g=0$ and $g=15$).

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}^{\rm s}\approx 5.622$. Note that color range changes between
different panels.

###  3.3. Macroscopic signal

When the duration $\delta$ of the gradient pulses is sufficiently long, only
few eigenmodes with small $\mathrm{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):

(14) $S\approx
e^{-2D_{0}\delta\lambda_{1}^{(g)}}\biggl{[}C_{1,1}^{(g)}+2\mathrm{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)})}\biggr{]},$

with the coefficients $C_{j,j^{\prime}}^{(g)}$ given by Eq. (6); note that we
used the property $\mathcal{B}_{-g}=\mathcal{B}_{g}^{*}$ that implies
$C_{j^{\prime},j}^{(g)}=[C_{j,j^{\prime}}^{(g)}]^{*}$. One can distinguish two
scenarios according to whether the eigenvalues $\lambda_{1}^{(g)}$ and
$\lambda_{2}^{(g)}$ are real or complex.

(i) When $0<g<g_{1}^{\rm 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. (14) can be neglected, yielding the one-mode approximation for the signal,

(15) $S\approx S_{\rm one}=C_{1,1}^{(g)}\,e^{-2D_{0}\delta\lambda_{1}^{(g)}}.$

Note that this approximation is not valid when $g$ is close to the branch
point $g_{1}^{\rm s}$.

(ii) When $g>g_{1}^{\rm 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. (14) can be written as

(16) $S\approx S_{\rm
two}=2e^{-2D_{0}\delta\,\mathrm{Re}\\{\lambda_{1}^{(g)}\\}}\biggl{[}C_{1,1}^{(g)}+\mathrm{Re}\bigl{\\{}C_{1,2}^{(g)}\,e^{2iD_{0}\delta\,\mathrm{Im}\\{\lambda_{1}^{(g)}\\}}\bigr{\\}}\biggr{]}.$

Figure 3 shows the dependence of the coefficients $C_{1,1}^{(g)}$ and
$\mathrm{Re}\\{C_{1,2}^{(g)}\\}$ on $g$. When $g$ approaches the branch point
$g_{1}^{\rm 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
$\mathrm{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}^{\rm s}$.

Figure 3.  The coefficients $C_{1,1}^{(g)}$ and
$\mathrm{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}^{\rm s}\approx 5.622$, at which
both coefficients diverge.

The accuracy of the approximations (15, 16) 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. (16). 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, 26, 20], so that the last factor in Eq. (16)
is approximately $e^{2i\gamma GR\delta}$. One thus retrieves diffusion-
diffraction patterns [51, 52, 53, 54, 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 the next-order corrections
to $\mathrm{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]:

(17)
$\mathrm{Re}\\{\lambda_{1}^{(g)}\\}=\frac{|a^{\prime}_{1}|}{2\ell_{g}^{2}}+\frac{1}{\sqrt{R}\ell_{g}^{3/2}}-\frac{\sqrt{3}}{4|a^{\prime}_{1}|R\ell_{g}}+O(\ell_{g}^{-1/2})\qquad(g\to\infty),$

where $\ell_{g}=(\gamma G/D_{0})^{-1/3}=g^{-1/3}$, and
$a^{\prime}_{1}\approx-1.02$ is the first zero of the derivative of the Airy
function $\rm{Ai}(z)$. Substituting the leading order of Eq. (17) into Eq.
(16), one retrieves the stretched exponential decay of the signal: $\ln
S\propto G^{2/3}$.

Figure 4.  Signal as a function of the gradient pulse duration $\delta$ for a
sphere of radius $R=10~{}\mu$m, with $\gamma=2.675\cdot 10^{8}$ rad/T/s
(protons), $D_{0}=2.3\cdot 10^{-9}~{}\rm{m}^{2}/\rm{s}$ (water molecules),
$G=17$ mT/m (a) and $G=129$ 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 (15, 16) 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$ 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$ mT/m (or $g=15$, panel (b)).

##  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
$\bm{G}=G(\bm{e}_{x}\cos\eta+\bm{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

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

where

(19) $\mathcal{B}_{g}^{\rm
d}=-\biggl{(}\underbrace{\partial_{r}^{2}+\frac{1}{r}\partial_{r}+\frac{1}{r^{2}}\partial_{\theta}^{2}}_{=\Delta_{\rm
d}}\biggr{)}+igr\cos\theta,\qquad\mathcal{B}_{g}^{\rm
i}=-\partial_{z}^{2}+igz$

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

The matrix elements of the operator $\mathcal{B}_{g}^{(\eta)}$ are derived in
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,\ldots$ being the order of
the Bessel function $J_{n}(z)$ of the first kind, $k=0,1,2,\ldots$ being the
index of the zeros $\alpha_{nk}$ of $J^{\prime}_{n}(z)$, $l$ distinguishing
between $s_{1}(z)=\cos(z)$ and $s_{2}(z)=\sin(z)$, and $m=0,1,2,\ldots$
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.

$j$ | 1 | 2-3 | 4-5 | 6 | 7-8 | 9 | 10-11 | 12-13
---|---|---|---|---|---|---|---|---
$nklm$ | 0010 | 10(1-2)0 | 20(1-2)0 | 0011 | 10(1-2)1 | 0110 | 30(1-2)0 | 20(1-2)1
$\lambda_{j}^{(0)}$ | 0 | 3.39 | 9.33 | 9.87 | 13.26 | 14.68 | 17.65 | 19.20
Table 2.  First 13 eigenvalues $\lambda_{nklm}=\lambda_{j}^{(0)}$ of the
(negative) Laplace operator, $\mathcal{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)}$.

###  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
$\mathcal{B}_{g}^{(\pi/2)}=\mathcal{B}_{g}^{z}=\mathcal{B}_{0}^{\rm
d}+\mathcal{B}_{g}^{\rm i}$, or lies in the transverse $xy$ plane so that
$\mathcal{B}_{g}^{(0)}=\mathcal{B}_{g}^{x}=\mathcal{B}_{g}^{\rm
d}+\mathcal{B}_{0}^{\rm i}$.

Figure 5.  Real part of the first 13 eigenvalues of the Bloch-Torrey operators
$\mathcal{B}_{g}^{z}$ (a) and $\mathcal{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}^{\rm i}\approx 18.06$ (a) and $g_{1}^{\rm d}\approx 3.76$, $g_{2}^{\rm
d}\approx 9.39$, $g_{3}^{\rm 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$.

Figure 5 presents real parts of the first 13 eigenvalues of the Bloch-Torrey
operators $\mathcal{B}_{g}^{z}$ (panel (a)) and $\mathcal{B}_{g}^{x}$ (panel
(b)). Let us first inspect the spectrum of $\mathcal{B}_{g}^{z}$, which is the
sum of the (negative) Laplace operator $\mathcal{B}_{0}^{\rm d}=-\Delta_{\rm
d}$ in the disk and the Bloch-Torrey operator $\mathcal{B}_{g}^{\rm i}$ on the
interval $(-H/2,H/2)$. The spectrum of $\mathcal{B}_{g}^{\rm i}$ was
thoroughly investigated in [21, 18, 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 range
of $g$. Moreover, Stoller et al. studied the branch points of
$\mathcal{B}_{g}^{\rm i}$ and found an explicit formula [21], which reads in
our notations (see also discussion in [31]):

(20) $g_{k}^{\rm i}=\sqrt{3}\frac{27}{4}j_{k}^{2}\,,\qquad\textrm{where}\quad
J_{-2/3}(j_{k})=0\quad(k=1,2,\ldots).$

In particular, one gets $g_{1}^{\rm i}\approx 18.06$ and $g_{2}^{\rm i}\approx
229.35$. The position of the first branch point $g_{1}^{\rm 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_{\rm 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}^{\rm i},g_{2}^{\rm
i},\ldots$. Note that the “shifted” eigenvalues correspond to the
eigenfunctions $u_{nkl}^{\rm d}$ of $-\Delta_{\rm d}$ that are orthogonal to
$u_{001}^{\rm d}=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
$\mathcal{B}_{g}^{x}$, which is the sum of $\mathcal{B}_{g}^{\rm d}$ and the
second derivative $\mathcal{B}_{0}^{\rm i}=-\partial_{z}^{2}$ on the interval.
The spectrum of the former operator is simply replicated and shifted
vertically by the eigenvalues $\pi^{2}m^{2}/H^{2}$ ($m=1,2,\ldots$) 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}^{\rm d}\approx
3.76$, $g_{2}^{\rm d}\approx 9.39$, and $g_{3}^{\rm d}\approx 13.87$, which
are replicated along the vertical axis by adding $\pi^{2}m^{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. (18) was independent of $g$ and thus just shifted
vertically the spectrum of the other. For intermediate angles $\eta$, both
terms in Eq. (18) 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}^{\rm d}$ and $\mathcal{B}_{g\sin\eta}^{\rm i}$, and
see how their features change. Most importantly, even though the
eigenfunctions are still factored along the longitudinal and transverse
directions $\bm{e}_{z}$ and $\bm{e}_{x}$, the nonzero gradients along
$\bm{e}_{x}$ and $\bm{e}_{z}$ directions break the symmetries of these factors
so that all eigenfunctions may contribute to the signal.

Figure 6.  Real part of the first 13 eigenvalues of the Bloch-Torrey operator
$\mathcal{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}^{\rm d}\approx
3.76/\cos(\eta)$, $g_{2}^{\rm d}\approx 9.36/\cos(\eta)$, $g_{3}^{\rm
d}\approx 13.87/\cos(\eta)$, and $g_{1}^{\rm i}\approx 18.06/\sin(\eta)$
associated to the operators $\mathcal{B}_{g}^{\rm d}$ and
$\mathcal{B}_{g}^{\rm i}$. Note that each pair of twice degenerate eigenvalues
appears as a single branch on panel (b). Truncation order was $320$.

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}^{\rm
d}$ and $\mathcal{B}_{g\sin\eta}^{\rm 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 (Fig. 7), in which the angle $\eta$
is tuned to make equal the rescaled branch points $g_{1}^{\rm d}/\cos\eta$ and
$g_{1}^{\rm i}/\sin\eta$, i.e., $\eta=\tan^{-1}(g_{1}^{\rm d}/g_{1}^{\rm
i})\approx 1.3661\approx 78.3^{\circ}$. The position $g_{1}^{\rm c}=g_{1}^{\rm
d}/\cos\eta\approx 18.5$ of the first group of branch points is indicated by a
vertical line. One can see more branches that merge at a single branch point
than in previous two panels. 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}^{\rm 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 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.

Figure 7.  Real part of the first 16 eigenvalues $\lambda_{j}^{(g)}$ of the
Bloch-Torrey operator $\mathcal{B}_{g}^{z}$ 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
$\mathrm{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 A), in which the matrices were
truncated at 320.

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, for a smooth boundary, the localization occurs at specific boundary
points $\bm{x}_{b}$, at which the normal vector $\bm{n}(\bm{x}_{b})$ is
parallel to the gradient $\bm{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
$\mathcal{B}_{g}^{\rm d}$ and an eigenfunction of $\mathcal{B}_{g}^{\rm i}$
for the disk and the interval, respectively. At high enough $g$, both factors
are localized: the eigenfunction of $\mathcal{B}_{g}^{\rm d}$ is localized at
either of two opposite points of the disk, while the eigenfunction of
$\mathcal{B}_{g}^{\rm i}$ is localized at either of two endpoints of the
interval. 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.

Figure 8.  $xz$ projection of the real part of four eigenfunctions
$v_{j}^{(g)}$ of the Bloch-Torrey operator $\mathcal{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.

###  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

(21) $\mathcal{B}_{g}=\frac{1}{R^{2}}\bigl{(}-\bar{\Delta}_{\rm
d}+ig\cos\eta\,R^{3}\,\bar{r}\cos\theta\bigr{)}+\frac{1}{H^{2}}\bigl{(}-\partial^{2}_{\bar{z}}+ig\sin\eta\,H^{3}\,\bar{z}\bigr{)},$

where bar denotes rescaled quantities: $\bar{r}=r/R$, $\bar{z}=z/H$,
$\bar{\Delta}_{\rm d}=R^{2}\Delta_{\rm 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

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

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 (22). The structure of the eigenfunctions of
$\mathcal{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
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 $\mathcal{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 $\mathcal{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
$\mathcal{B}_{g}^{\rm d}$ and $\mathcal{B}_{g}^{\rm 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, 58, 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^{\prime}}^{(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 $\mathcal{B}_{g}^{x}$, $\mathcal{B}_{g}^{y}$ and
$\mathcal{B}_{g}^{z}$ for three orthogonal directions allows one to analyze
double-pulsed field-gradient experiments [60, 61, 62, 15, 63, 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.

## 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 $\mathcal{B}_{g}$. It is inspired from the matrix
formalism [12, 6, 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 the radius of the capped cylinder.

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

(23) $v_{j}^{(g)}(\bm{x})=\sum\limits_{k}X_{j,k}^{(g)}u_{k}(\bm{x}),$

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

$\lambda_{j}^{(g)}\sum\limits_{k}X_{j,k}^{(g)}u_{k}(\bm{x})=\lambda_{j}^{(g)}v_{j}^{(g)}(\bm{x})=(-\Delta+igx)v_{j}^{(g)}(\bm{x})=\sum\limits_{k}X_{j,k}^{(g)}(-\Delta+igx)u_{k}(\bm{x}).$

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

(24)
$R^{2}\lambda_{j}^{(g)}X_{j,j^{\prime}}^{(g)}=\sum\limits_{k}X_{j,k}^{(g)}(\Lambda_{k,j^{\prime}}+i\bar{g}B_{k,j^{\prime}}),$

where $\Lambda_{k,j^{\prime}}=\delta_{k,j^{\prime}}\lambda_{k}R^{2}$,
$B_{k,j^{\prime}}=\int\nolimits_{\Omega}d\bm{x}\,u_{k}(\bm{x})\,(x/R)\,u^{*}_{j^{\prime}}(\bm{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

(25) $\Lambda^{(g)}X=X(\Lambda+i\bar{g}B),$

where $\Lambda^{(g)}$ is the diagonal matrix of eigenvalues
$R^{2}\lambda_{n}^{(g)}$ of the Bloch-Torrey operator $\mathcal{B}_{g}$. As a
consequence, a numerical 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 (23) that implies

(26)
$\int\limits_{\Omega}d\bm{x}\,v_{j}^{(g)}(\bm{x})\,v_{j^{\prime}}^{(g)}(\bm{x})=[XWX^{\dagger}]_{j,j^{\prime}}=\delta_{j,j^{\prime}},$

where

(27)
$W_{k,k^{\prime}}=\int\limits_{\Omega}d\bm{x}\,u_{k}(\bm{x})\,u_{k^{\prime}}(\bm{x}).$

Since there is no complex-conjugate in Eq. (27), the matrix $W$ is not
necessarily identity (see, e.g., Eq. (46) for a sphere). We recall that the
integral in Eq. (26) may be zero for $j=j^{\prime}$ 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 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^{\prime}}^{(g)}$ from the
spectral expansion (5) of the signal can be written as

(28)
$C_{j,j^{\prime}}^{(g)}=\mu_{j}^{(-g)}\,\Gamma_{j,j^{\prime}}^{(g)}\,\mu_{j^{\prime}}^{(g)},$

where

(29)
$\mu_{j}^{(g)}=\frac{1}{\sqrt{|\Omega|}}\int\limits_{\Omega}d\bm{x}\,v_{j}^{(g)}(\bm{x})$

is the projection of the eigenfunction $v_{j}^{(g)}(\bm{x})$ onto a constant,
and

(30)
$\Gamma_{j,j^{\prime}}^{(g)}=\int\limits_{\Omega}d\bm{x}\,v_{j}^{(-g)}(\bm{x})\,v_{j^{\prime}}^{(g)}(\bm{x}).$

is the overlap between two eigenfunctions $v_{j}^{(-g)}(\bm{x})$ and
$v_{j^{\prime}}^{(g)}(\bm{x})$. The integrals in Eqs. (29, 30) can be computed
directly by using the representation (23):

(31)
$\mu_{j}=\frac{1}{\sqrt{|\Omega|}}\int\limits_{\Omega}d\bm{x}\sum\limits_{k}X_{j,k}u_{k}(\bm{x})=[XU]_{j},$

where

(32)
$U_{k}=\frac{1}{\sqrt{|\Omega|}}\int\limits_{\Omega}d\bm{x}\,u_{k}(\bm{x}).$

Since $u_{k}(\bm{x})$ are orthogonal to $u_{0}(\bm{x})=1/\sqrt{|\Omega|}$ for
Neumann boundary condition, one gets

(33) $\mu_{j}=X_{j,0}.$

Similarly,

(34)
$\Gamma_{j,j^{\prime}}=\sum\limits_{k}X_{j,k}^{*}\sum\limits_{k^{\prime}}X_{j^{\prime},k^{\prime}}\int\limits_{\Omega}d\bm{x}\,u_{k}^{*}(\bm{x})\,u_{k^{\prime}}(\bm{x})=[X^{*}X^{\dagger}]_{j,j^{\prime}}.$

##  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{\mathcal{B}}_{g}$, the matrix representation
was derived in [6, 46]:

(35)
$\Lambda_{nk,n^{\prime}k^{\prime}}=\delta_{n,n^{\prime}}\delta_{k,k^{\prime}}\frac{\alpha_{nk}^{2}}{R^{2}}\,,$

where $\alpha_{nk}$ are the positive zeros of $j^{\prime}_{n}(z)$ (with
$n=0,1,2,\ldots$), enumerated by $k=0,1,2,\ldots$, and

(36) $B_{nk,n^{\prime}k^{\prime}}=\delta_{n,n^{\prime}\pm
1}\frac{n+n^{\prime}+1}{(2n+1)(2n^{\prime}+1)}\beta_{nk}\beta_{n^{\prime}k^{\prime}}\frac{\alpha_{nk}^{2}+\alpha_{n^{\prime}k}^{2}-n(n^{\prime}+1)-n^{\prime}(n+1)+1}{(\alpha_{nk}^{2}-\alpha_{n^{\prime}k^{\prime}}^{2})^{2}}\,,$

with

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

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 $\mathcal{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

(38)
$\Lambda_{nkm,n^{\prime}k^{\prime}m^{\prime}}=\delta_{m,m^{\prime}}\Lambda_{nk,n^{\prime}k^{\prime}}=\delta_{m,m^{\prime}}\delta_{n,n^{\prime}}\delta_{k,k^{\prime}}\frac{\alpha_{nk}^{2}}{R^{2}}\,.$

In turn, the matrix $B$ representing the gradient term, depends on the
direction $\bm{e}_{G}$ of the gradient $\bm{G}=G\bm{e}_{G}$. Encoding this
direction in spherical coordinates by angles $\theta_{G}$ and $\phi_{G}$ as

(39)
$\bm{e}_{G}=\sin\theta_{G}\cos\phi_{G}\bm{e}_{x}+\sin\theta_{G}\sin\phi_{G}\bm{e}_{y}+\cos\theta_{G}\bm{e}_{z},$

one can represent this gradient by the matrix

$\displaystyle B_{nkm,n^{\prime}k^{\prime}m^{\prime}}$
$\displaystyle=\frac{1}{R}\int\limits_{\Omega}d\bm{x}\,[u_{nkm}(\bm{x})]^{*}(\bm{e}_{G}\cdot\bm{x})u_{n^{\prime}k^{\prime}m^{\prime}}(\bm{x})$
$\displaystyle=\biggl{[}\sin\theta_{G}\cos\phi_{G}B^{x}+\sin\theta_{G}\sin\phi_{G}B^{y}+\cos\theta_{G}B^{z}\biggr{]}_{nkm,n^{\prime}k^{\prime}m^{\prime}},$

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

(41)
$B_{nkm,n^{\prime}k^{\prime}m^{\prime}}^{z}=\delta_{m,m^{\prime}}B_{nk,n^{\prime}k^{\prime}}\sqrt{1-\frac{m^{2}}{(\max\\{n,n^{\prime}\\})^{2}}}\qquad(|m|\leq\min\\{n,n^{\prime}\\}).$

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

(42) $\displaystyle B_{nkm,(n+1)k^{\prime}m^{\prime}}^{x}$
$\displaystyle=\frac{B_{nk,(n+1)k^{\prime}}}{2}\left(\delta_{m^{\prime},m-1}\frac{\sqrt{(n-m+1)(n-m+2)}}{n+1}\right.$
$\displaystyle\left.-\delta_{m^{\prime},m+1}\frac{\sqrt{(n+m+1)(n+m+2)}}{n+1}\right),$
(43) $\displaystyle B_{nkm,(n-1)k^{\prime}m^{\prime}}^{x}$
$\displaystyle=-\frac{B_{nk,(n-1)k^{\prime}}}{2}\left(\delta_{m^{\prime},m-1}\frac{\sqrt{(n+m-1)(n+m)}}{n}\right.$
$\displaystyle\left.-\delta_{m^{\prime},m+1}\frac{\sqrt{(n-m-1)(n-m)}}{n}\right),$

and

(44) $\displaystyle B_{nkm,(n+1)k^{\prime}m^{\prime}}^{y}$
$\displaystyle=i\frac{B_{nk,(n+1)k^{\prime}}}{2}\left(\delta_{m^{\prime},m-1}\frac{\sqrt{(n-m+1)(n-m+2)}}{n+1}\right.$
$\displaystyle\left.+\delta_{m^{\prime},m+1}\frac{\sqrt{(n+m+1)(n+m+2)}}{n+1}\right),$
(45) $\displaystyle B_{nkm,(n-1)k^{\prime}m^{\prime}}^{y}$
$\displaystyle=-i\frac{B_{nk,(n-1)k^{\prime}}}{2}\left(\delta_{m^{\prime},m-1}\frac{\sqrt{(n+m-1)(n+m)}}{n}\right.$
$\displaystyle\left.+\delta_{m^{\prime},m+1}\frac{\sqrt{(n-m-1)(n-m)}}{n}\right),$

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

(46)
$W_{nkm,n^{\prime}k^{\prime}m^{\prime}}=(-1)^{m}\delta_{n,n^{\prime}}\delta_{k,k^{\prime}}\delta_{m,-m^{\prime}}.$

##  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=\\{\bm{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 gradients, it is therefore
enough to consider two reduced Bloch-Torrey operators:
$\hat{\mathcal{B}}^{z}_{g}=-\partial_{z}^{2}+igz$ on the interval
$(-H/2,H/2)$, and
$\hat{\mathcal{B}}^{xy}_{g}=-(\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 [46, 6]. The spectral properties of $\hat{\mathcal{B}}^{z}_{g}$
were thoroughly investigated in [21, 18, 24], while the spectrum of
$\hat{\mathcal{B}}^{xy}_{g}$ was discussed in [17, 31] (see references
therein). However, more sophisticated pulsed-gradient sequences with several
gradients directions require the knowledge of the whole Bloch-Torrey operator
$\mathcal{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

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

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

(48) $u_{nkl}^{\rm
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\left\\{\begin{array}[]{l
l}\cos(n\theta)&(l=1),\\\ \sin(n\theta)&(l=2),\\\ \end{array}\right.$

where $J_{n}(z)$ is the Bessel function of the first kind, $\alpha_{nk}$ are
the positive zeros of $J^{\prime}_{n}(z)$ enumerated by $k=0,1,2,\ldots$, and

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

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

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

In general, the eigenvalue $\lambda_{nklm}$ is twice degenerate for $n>0$ and
simple for $n=0$ (in this case, $u_{0k2}^{\rm 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:

(56) $\displaystyle\Lambda$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{ccccc}\Lambda_{\rm d}&0&0&0&\ldots\\\
0&\Lambda_{\rm d}+(\pi^{2}/H^{2})\mathbb{I}&0&0&\ldots\\\ 0&0&\Lambda_{\rm
d}+(4\pi^{2}/H^{2})\mathbb{I}&0&\ldots\\\ 0&0&0&\Lambda_{\rm
d}+(9\pi^{2}/H^{2})\mathbb{I}&\ldots\\\ \ldots&\ldots&\ldots&\ldots&\ldots\\\
\end{array}\right),$ (62) $\displaystyle B^{x,y}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{ccccc}B^{x,y}_{\rm d}&0&0&0&\ldots\\\
0&B^{x,y}_{\rm d}&0&0&\ldots\\\ 0&0&B^{x,y}_{\rm d}&0&\ldots\\\
0&0&0&B^{x,y}_{\rm d}&\ldots\\\ \ldots&\ldots&\ldots&\ldots&\ldots\\\
\end{array}\right)$ (68) $\displaystyle B^{z}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{ccccc}B^{\rm i}_{0,0}\mathbb{I}&B^{\rm
i}_{0,1}\mathbb{I}&B^{\rm i}_{0,2}\mathbb{I}&B^{\rm
i}_{0,3}\mathbb{I}&\ldots\\\ B^{\rm i}_{1,0}\mathbb{I}&B^{\rm
i}_{1,1}\mathbb{I}&B^{\rm i}_{1,2}\mathbb{I}&B^{\rm
i}_{1,3}\mathbb{I}&\ldots\\\ B^{\rm i}_{2,0}\mathbb{I}&B^{\rm
i}_{2,1}\mathbb{I}&B^{\rm i}_{2,2}\mathbb{I}&B^{\rm
i}_{2,3}\mathbb{I}&\ldots\\\ B^{\rm i}_{3,0}\mathbb{I}&B^{\rm
i}_{3,1}\mathbb{I}&B^{\rm i}_{3,2}\mathbb{I}&B^{\rm
i}_{3,3}\mathbb{I}&\ldots\\\ \ldots&\ldots&\ldots&\ldots&\ldots\\\
\end{array}\right),$

where $\mathbb{I}$ is the identity matrix, $B^{\rm i}_{m,m^{\prime}}$ are the
matrix elements for the interval:

(69) $B^{\rm
i}_{m,m^{\prime}}=((-1)^{m+m^{\prime}}-1)\sqrt{2-\delta_{m,0}}\sqrt{2-\delta_{m^{\prime},0}}\,\frac{m^{2}+m^{\prime
2}}{\pi^{2}(m^{2}-m^{\prime 2})^{2}}\quad(m\neq m^{\prime}),$

and $B^{\rm i}_{m,m}=0$, with $m,m^{\prime}=0,1,2,\ldots$. In turn,
$\Lambda_{\rm d}$ and $B_{\rm 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

(70) $[\Lambda_{\rm
d}]_{nkl,n^{\prime}k^{\prime}l^{\prime}}=\delta_{n,n^{\prime}}\delta_{k,k^{\prime}}\delta_{l,l^{\prime}}\frac{\alpha_{nk}^{2}}{R^{2}}$

and

(71) $[B_{\rm d}^{x}]_{nk1,n^{\prime}k^{\prime}1}=B^{\rm
d}_{nk,n^{\prime}k^{\prime}},\qquad[B_{\rm
d}^{x}]_{nk2,n^{\prime}k^{\prime}2}=B^{\rm
d}_{nk,n^{\prime}k^{\prime}}(1-\delta_{n+n^{\prime},1}),$ (72) $[B_{\rm
d}^{x}]_{nk1,n^{\prime}k^{\prime}2}=[B_{\rm
d}^{x}]_{nk2,n^{\prime}k^{\prime}1}=0,$

where

(73) $B^{\rm d}_{nk,n^{\prime}k^{\prime}}=\delta_{n,n^{\prime}\pm
1}(1+\delta_{n,0}+\delta_{n^{\prime},0})^{1/2}\beta_{nk}\beta_{n^{\prime}k^{\prime}}\frac{\alpha_{nk}^{2}+\alpha_{n^{\prime}k^{\prime}}^{2}-2nn^{\prime}}{(\alpha_{nk}^{2}-\alpha_{n^{\prime}k^{\prime}}^{2})^{2}}\,.$

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

(74) $[B_{\rm d}^{y}]_{nk1,n^{\prime}k^{\prime}1}=[B_{\rm
d}^{y}]_{nk2,n^{\prime}k^{\prime}2}=0,$ (75) $[B_{\rm
d}^{y}]_{nk1,(n+1)k^{\prime}2}=B^{\rm d}_{nk,(n+1)k^{\prime}},\quad[B_{\rm
d}^{y}]_{nk1,(n-1)k^{\prime}2}=-B^{\rm
d}_{nk,(n-1)k^{\prime}}(1-\delta_{n+n^{\prime},1}),$

and

(76) $[B_{\rm d}^{y}]_{nk2,(n-1)k^{\prime}1}=B^{\rm
d}_{nk,(n-1)k^{\prime}},\quad[B_{\rm d}^{y}]_{nk2,(n+1)k^{\prime}1}=-B^{\rm
d}_{nk,(n+1)k^{\prime}}(1-\delta_{n+n^{\prime},1}).$

##  Appendix D Orthogonalization of eigenfunctions with degenerate eigenvalues

The relation (8) ensures the orthogonality of eigenfunctions $v_{j}^{(g)}$ and
$v_{j^{\prime}}^{(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},\ldots,j_{n}$ such that
$\lambda_{j_{1}}^{(g)}=\lambda_{j_{2}}^{(g)}=\ldots=\lambda_{j_{n}}^{(g)}$,
the associated eigenfunctions $v_{j_{i}}^{(g)}$ form an eigenspace of
dimension $n$, in which any $n$ independent linear 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^{\prime}}$ be two eigenfunctions with the same
eigenvalue (we dropped here the superscript (g) for brevity), and their
(non)-orthogonality is characterized by the matrix

(77) $C=\left(\begin{array}[]{cc}\langle v_{j},v_{j}\rangle&\langle
v_{j},v_{j^{\prime}}\rangle\\\ \langle v_{j^{\prime}},v_{j}\rangle&\langle
v_{j^{\prime}},v_{j^{\prime}}\rangle\\\ \end{array}\right).$

We aim at constructing two linear combinations,

(78)
$\hat{v}_{j}=av_{j}+bv_{j^{\prime}},\qquad\hat{v}_{j^{\prime}}=cv_{j}+dv_{j^{\prime}},$

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

(79) $\left(\begin{array}[]{c}\hat{v}_{j}\\\ \hat{v}_{j^{\prime}}\\\
\end{array}\right)=\left(\begin{array}[]{c c}a&b\\\ c&d\\\
\end{array}\right)\left(\begin{array}[]{c}v_{j}\\\ v_{j^{\prime}}\\\
\end{array}\right),$

so that

(80)
$\left(\begin{array}[]{cc}\langle\hat{v}_{j},\hat{v}_{j}\rangle&\langle\hat{v}_{j},\hat{v}_{j^{\prime}}\rangle\\\
\langle\hat{v}_{j^{\prime}},\hat{v}_{j}\rangle&\langle\hat{v}_{j^{\prime}},\hat{v}_{j^{\prime}}\rangle\\\
\end{array}\right)=\left(\begin{array}[]{c c}a&b\\\ c&d\\\
\end{array}\right)C\left(\begin{array}[]{c c}a&c\\\ b&d\\\
\end{array}\right).$

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

(81) $\frac{1}{(ad-bc)^{2}}\left(\begin{array}[]{c c}b^{2}+d^{2}&-ab-cd\\\
-ab-cd&a^{2}+c^{2}\\\ \end{array}\right)=C.$

To proceed, we parameterize the unknown coefficients as:

(82) $a=A^{-1}\cos\alpha,\quad b=A^{-1}\sin\alpha,\quad
c=-B^{-1}\sin\alpha,\quad d=B^{-1}\cos\alpha.$

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. (81) yields three
equations:

(83) $\displaystyle A^{2}\cos^{2}\alpha+B^{2}\sin^{2}\alpha$ $\displaystyle=$
$\displaystyle C_{1,1},$ (84) $\displaystyle
A^{2}\sin^{2}\alpha+B^{2}\cos^{2}\alpha$ $\displaystyle=$ $\displaystyle
C_{2,2},$ (85) $\displaystyle(A^{2}-B^{2})\sin(2\alpha)$ $\displaystyle=$
$\displaystyle 2C_{1,2}.$

One can solve these equations as

(86) $\displaystyle\alpha$ $\displaystyle=$
$\displaystyle\frac{1}{2}\mathrm{atan}\left(\frac{2C_{1,2}}{C_{1,1}-C_{2,2}}\right),$
(87) $\displaystyle A^{2}$ $\displaystyle=$
$\displaystyle\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}.$

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:

(88) $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)}.$

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

In practice, once the matrix $X$ of coefficients in Eq. (23) is found by
solving the eigenvalue problem (25), one can evaluate the matrix
$XWX^{\dagger}$ that represents the orthogonality of the eigenfunctions
$v_{j}^{(g)}$. According to Eq. (26), 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 non-
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 nondiagonal elements and apply the above orthonormalization procedure
for each such pair. This procedure was applied for most spectral computations
in this work.

## References

  * [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