Multidimensional encoding of restricted and anisotropic diffusion by double rotation of the q vector

Diffusion NMR and MRI methods building on the classic pulsed gradient spin-echo sequence are sensitive to many aspects of translational motion, including time and frequency dependence (“restriction”), anisotropy, and flow, leading to ambiguities when interpreting experimental data from complex heterogeneous materials such as living biological tissues. While the oscillating gradient technique specifically targets frequency dependence and permits control of the sensitivity to flow, tensor-valued encoding enables investigations of anisotropy in orientationally disordered materials. Here, we propose a simple scheme derived from the “double-rotation” technique in solid-state NMR to generate a family of modulated gradient waveforms allowing for comprehensive exploration of the 2D frequency–anisotropy space and convenient investigation of both restricted and anisotropic diffusion with a single multidimensional acquisition protocol, thereby combining the desirable characteristics of the oscillating gradient and tensor-valued encoding techniques. The method is demonstrated by measuring multicomponent isotropic Gaussian diffusion in simple liquids, anisotropic Gaussian diffusion in a polydomain lyotropic liquid crystal, and restricted diffusion in a yeast cell sediment.

1 Introduction

Magnetic field gradients applied during the dephasing and rephasing periods of a spin-echo sequence (Hahn, 1950) render the NMR signal sensitive to various aspects of translational motion including bulk diffusivity (Douglass and McCall, 1958), flow (Carr and Purcell, 1954), time and frequency dependence (“restriction”) (Woessner, 1963), anisotropy (Boss and Stejskal, 1965), and exchange (Kärger, 1969). Although the conventional and ubiquitous pulsed gradient spin-echo sequence by Stejskal and Tanner (1965) may give information about all of these aspects, more elaborate gradient modulations (Tanner, 1979; Cory et al., 1990; Callaghan and Manz, 1994; Mori and van Zijl, 1995) are required to unambiguously assign a certain mechanism to the experimental observations (Topgaard, 2017; Lundell and Lasič, 2020). Diffusion MRI methods incorporating such advanced diffusion encoding schemes have recently been shown to have potential for clinical research applications (Reymbaut et al., 2020); some notable examples are oscillating gradients to estimate cell sizes (Xu et al., 2021) and tensor-valued encoding to characterize cell shapes (Daimiel Naranjo et al., 2021) in breast tumors.

The sensitivity of the MRI signal to the various types of motion can be quantified with the tensor-valued encoding spectrum b(ω) (Topgaard, 2019a; Lundell and Lasič, 2020), the trace of which equals the dephasing power spectrum (Stepišnik, 1981) – relevant for isotropic restricted diffusion – and whose integral over ω equals the conventional b matrix (Basser et al., 1994) giving information about diffusion anisotropy. While most studies focus on either the frequency-dependent (Aggarwal, 2020) or tensorial (Reymbaut, 2020) aspects of the encoding, Lundell et al. (2019) suggested joining them into a common multidimensional framework. The approach was demonstrated with gradient waveforms deriving from the magic-angle spinning (MAS) technique in solid-state NMR spectroscopy (Andrew et al., 1959; Eriksson et al., 2013; Topgaard, 2013); however, these methods offer only limited access to the frequency and anisotropy dimensions.

Expanding on the results of Lundell et al. (2019), we take inspiration from the “double-rotation” (DOR) technique in solid-state NMR (Samoson et al., 1998) and derive a family of gradient waveforms for comprehensive exploration of, in particular, the frequency–anisotropy dimensions of b(ω), as quantified by the centroid frequency ω_cent (Arbabi et al., 2020) and encoding anisotropy b_Δ (Eriksson et al., 2015), in addition to the b value and b vector (Θ, Φ) of conventional diffusion tensor imaging (Kingsley, 2006). While ω_cent is key for characterizing restricted diffusion (Stepišnik and Callaghan, 2000), the variable b_Δ enables quantification of anisotropy in orientationally disordered materials (Eriksson et al., 2015) and estimation of nonparametric diffusion tensor distributions (de Almeida Martins and Topgaard, 2016; Topgaard, 2019b). The ability of the new gradient waveforms to give access to the complete 2D ω_cent–b_Δ plane is demonstrated by microimaging measurements on previously studied phantoms with well-defined restriction and anisotropy properties, namely water (Mills, 1973) and concentrated salt solution (Wadsö et al., 2009) with isotropic Gaussian diffusion, a lamellar liquid crystal giving anisotropic Gaussian diffusion (Topgaard, 2016), and a yeast cell sediment exhibiting isotropic restricted diffusion (Malmborg et al., 2006).

2 Theoretical background

To set the stage for later sections dealing with the proposed gradient waveforms to investigate both frequency-dependent and tensorial aspects of translational motion, we include a brief summary of the relevant theory here; greater detail can be found in the textbooks by Price (2009) and Callaghan (2011) as well as in the comprehensive review by Lundell and Lasič (2020). Readers already familiar with the background material may proceed directly to the design of gradient waveforms in Sect. 3 after noting Eqs. (32) and (34) with the definitions of the main variables ω_cent and b_Δ reporting on the sensitivity to restriction and anisotropy.

2.1 Encoding of translational motion by magnetic field gradients

Figure 1 illustrates the effects of a general gradient waveform g(t) on the NMR signal from an ensemble of spins undergoing restricted diffusion and flow within an infinite cylinder. As shown using the 2D and 3D plots of g(t) in Fig. 1a, both the magnitude and direction of the gradient vector are changing smoothly with time. Simultaneously, the spins spread out and gradually drift from their initial positions (Fig. 1b). The time-dependent normalized signal E(t) is given by

$$\begin{array}{}\text{(1)}& E\left(t\right)=〈\exp \left(\mathrm{i}\mathit{\varphi }\left(t\right)\right)〉,\end{array}$$

where 〈…〉 denotes an ensemble mean and the time-dependent phase ϕ(t) of a single spin with gyromagnetic ratio γ is determined by the time integral of the scalar product between g(t) the time-dependent position r(t) according to

$$\begin{array}{}\text{(2)}& \mathit{\varphi }\left(t\right)=-\mathit{\gamma }{\int }_{\mathrm{0}}^t\mathbf{g}\left(t^{\prime }\right)\cdot \mathbf{r}\left(t^{\prime }\right)\mathrm{d}{t}^{\prime }.\end{array}$$

The interplay between g(t) and r(t) results in ϕ(t) evolving from zero for all spins at t=0 to periodic patterns with varying directions and spatial wavelengths at intermediate times and, finally, an overall phase shift superposed on partially randomized values. From Eq. (1), it follows that the latter phase dispersion leads to a decrease in the magnitude of the signal.

The evolution of ϕ(t) may be rationalized by partial integration of Eq. (2) into

$$\begin{array}{}\text{(3)}& \mathit{\varphi }\left(t\right)=-\mathbf{q}\left(t\right)\cdot \mathbf{r}\left(t\right)+{\int }_{\mathrm{0}}^t\mathbf{q}\left(t^{\prime }\right)\cdot \mathbf{v}\left(t^{\prime }\right)\mathrm{d}{t}^{\prime },\end{array}$$

where q(t) is the dephasing vector, defined as

$$\begin{array}{}\text{(4)}& \mathbf{q}\left(t\right)=\mathit{\gamma }{\int }_{\mathrm{0}}^t\mathbf{g}\left(t^{\prime }\right)\mathrm{d}{t}^{\prime },\end{array}$$

and v(t) = dr(t) $/$ dt is the time-dependent velocity. The spatial periodicity in ϕ(t) at intermediate times is, according to the first term in Eq. (3), given by the scalar product between q(t) and r(t), which is utilized to obtain the spatial resolution in MRI where the dephasing vector is usually denoted k(t). The plots of q(t) in Fig. 1a show smooth changes in both the magnitude and direction of the vector with time.

Figure 1 Principles of motion encoding by general gradient waveforms. (a) Time-dependent gradient g(t) and dephasing vector q(t) illustrated as 2D graphs of Cartesian components vs. t (left) and 3D plots of trajectories through space (right). The average values of the q(t) components within the time interval from t=0 to τ are indicated with arrows labeled ${\mathbf{q}}_{\mathrm{v}}/\mathit{\tau }$, where the velocity-encoding vector q_v is defined in Eq. (15). (b) Sequence of snapshots from a random walk simulation of an ensemble of spins (spheres) undergoing restricted diffusion and flow within an infinite cylinder (black section) during application of the waveforms in panel (a). Each spin is color-coded by its time-dependent phase ϕ(t) given by the interplay between g(t) and the spin positions r(t) according to Eq. (2). (c) Phase distribution P(ϕ(τ)) for the ensemble of spins at t=τ (black histogram) and a Gaussian (smooth gray line) with mean 〈ϕ(τ)〉 and standard deviation std[ϕ(τ)] (spacing between vertical gray lines) which give the phase shift α and attenuation factor β of the signal E(τ) via Eqs. (8), (9), and (11). (d) Frequency-dependent elements (color-coded) of the tensor-valued encoding spectrum b(ω) and diffusion spectrum D(ω) as well as their tensor dot product b(ω) : D(ω) which gives β via Eq. (21). The centroid frequency ω_cent (arrow) is obtained from the trace of b(ω) (gray line) via Eq. (32). (e) The 3D plots of the ensemble mean velocity 〈v〉, velocity-encoding vector q_v, and encoding tensor b, with the latter being defined in Eq. (33). The scalar product of 〈v〉 and q_v gives α via Eq. (14).

Focusing on translational displacements rather than absolute positions, we select a time τ where

$$\begin{array}{}\text{(5)}& \mathbf{q}\left(\mathit{\tau }\right)=\mathrm{0}\end{array}$$

and the first term in Eq. (3) vanishes while the second one remains:

$$\begin{array}{}\text{(6)}& \mathit{\varphi }\left(\mathit{\tau }\right)={\int }_{\mathrm{0}}^{\mathit{\tau }}\mathbf{q}\left(t\right)\cdot \mathbf{v}\left(t\right)\mathrm{d}t.\end{array}$$

The value of ϕ(τ) in Eq. (6) is insensitive to r(τ) but depends on the history of q(t) and v(t) in the interval from t=0 to τ.

2.2 Gaussian phase distribution approximation

As shown in Fig. 1c, the selected gradient waveform and random walk simulation parameters yield a phase distribution P(ϕ(t)) that is well approximated at t=τ as a Gaussian function with mean 〈ϕ(τ)〉 and standard deviation std[ϕ(τ)] = (〈ϕ(τ)²〉 – 〈ϕ(τ)〉²)${}^{\mathrm{1}/\mathrm{2}}$. The Gaussian function can be expressed as

$$\begin{array}{}\text{(7)}& P\left(\mathit{\varphi }\left(\mathit{\tau }\right)\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\sqrt{\mathit{\pi }\mathit{\beta }}}}\exp \left(-{\displaystyle \frac{{\left(\mathit{\varphi }\left(\mathit{\tau }\right)-\mathit{\alpha }\right)}^{\mathrm{2}}}{\mathrm{4}\mathit{\beta }}}\right),\end{array}$$

where

$$\begin{array}{}\text{(8)}& \mathit{\alpha }=〈\mathit{\varphi }\left(\mathit{\tau }\right)〉\end{array}$$

and

$$\begin{array}{}\text{(9)}& \mathit{\beta }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(〈\mathit{\varphi }{\left(\mathit{\tau }\right)}^{\mathrm{2}}〉-{〈\mathit{\varphi }\left(\mathit{\tau }\right)〉}^{\mathrm{2}}\right).\end{array}$$

After rewriting Eq. (1) as an integral,

$$\begin{array}{}\text{(10)}& E\left(\mathit{\tau }\right)={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}P\left(\mathit{\varphi }\left(\mathit{\tau }\right)\right)\exp \left(\mathrm{i}\mathit{\varphi }\left(\mathit{\tau }\right)\right)\mathrm{d}\mathit{\varphi }\left(\mathit{\tau }\right),\end{array}$$

the insertion of Eq. (7) into Eq. (10) may be evaluated to

$$\begin{array}{}\text{(11)}& E\left(\mathit{\tau }\right)=\exp \left(\mathrm{i}\mathit{\alpha }-\mathit{\beta }\right),\end{array}$$

where α and β can be identified as quantitative measures of the overall phase shift and attenuation of the signal, respectively, as previously deduced from visual inspection of the phases of the spin ensemble in Fig. 1b. The Gaussian phase distribution approximation has been applied for the cases of free (Carr and Purcell, 1954; Douglass and McCall, 1958) and restricted (Neuman, 1974) diffusion, and investigations of its ranges of validity can be found in the literature (Balinov et al., 1993; Stepišnik, 1999).

2.3 Mean velocity and velocity correlation function

Insertion of Eq. (6) into Eq. (8) yields

$$\begin{array}{}\text{(12)}& \mathit{\alpha }=〈{\int }_{\mathrm{0}}^{\mathit{\tau }}\mathbf{q}\left(t\right)\cdot \mathbf{v}\left(t\right)\mathrm{d}t〉,\end{array}$$

which by separating v(t) into the ensemble mean 〈v(t)〉 = 〈v〉 and fluctuating part u(t), defined by

$$\begin{array}{}\text{(13)}& \mathbf{u}\left(t\right)=\mathbf{v}\left(t\right)-〈\mathbf{v}〉,\end{array}$$

can be evaluated to

$$\begin{array}{}\text{(14)}& \mathit{\alpha }={\mathbf{q}}_{\mathrm{v}}\cdot 〈\mathbf{v}〉,\end{array}$$

where the flow encoding vector q_v is defined as

$$\begin{array}{}\text{(15)}& {\mathbf{q}}_{\mathrm{v}}={\int }_{\mathrm{0}}^{\mathit{\tau }}\mathbf{q}\left(t\right)\mathrm{d}t.\end{array}$$

Correspondingly, the insertion of Eq. (6) into Eq. (9) gives

$$\begin{array}{}\text{(16)}& \mathit{\beta }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(〈{\left[{\int }_{\mathrm{0}}^{\mathit{\tau }}\mathbf{q}\left(t\right)\cdot \mathbf{v}\left(t\right)\mathrm{d}t\right]}^{\mathrm{2}}〉-{〈{\int }_{\mathrm{0}}^{\mathit{\tau }}\mathbf{q}\left(t\right)\cdot \mathbf{v}\left(t\right)\mathrm{d}t〉}^{\mathrm{2}}\right),\end{array}$$

which by reordering the time integrals and ensemble means as well as noting that 〈u(t) ⋅ 〈v〉〉=0, can be expressed as

$$\begin{array}{}\text{(17)}& \mathit{\beta }={\int }_{\mathrm{0}}^{\mathit{\tau }}{\int }_{\mathrm{0}}^t\mathbf{q}{\left(t\right)}^{\mathrm{T}}\cdot 〈\mathbf{u}\left(t\right)\mathbf{u}{\left(t{}^{\prime }\right)}^{\mathrm{T}}〉\cdot \mathbf{q}\left(t{}^{\prime }\right)\mathrm{d}t\mathrm{d}t{}^{\prime },\end{array}$$

where 〈u(t)u(t^′)^T〉 is the tensor-valued velocity correlation function.

2.4 Transformation to the frequency domain

After introducing the dephasing spectrum q(ω) and diffusion spectrum D(ω) by Fourier transformations to the frequency (ω) domain according to

$$\begin{array}{}\text{(18)}& \mathbf{q}\left(\mathit{\omega }\right)={\int }_{\mathrm{0}}^{\mathit{\tau }}\mathbf{q}\left(t\right)\exp \left(\mathrm{i}\mathit{\omega }t\right)\mathrm{d}t\end{array}$$

and

$$\begin{array}{}\text{(19)}& \mathbf{D}\left(\mathit{\omega }\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}〈\mathbf{u}\left(t\right)\mathbf{u}{\left(t{}^{\prime }\right)}^{\mathrm{T}}〉\exp \left(\mathrm{i}\mathit{\omega }\left(t{}^{\prime }-t\right)\right)\mathrm{d}\left(t{}^{\prime }-t\right),\end{array}$$

Eq. (17) can be recast into

$$\begin{array}{}\text{(20)}& \mathit{\beta }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\mathbf{q}{\left(\mathit{\omega }\right)}^{\mathrm{T}}\cdot \mathbf{D}\left(\mathit{\omega }\right)\cdot \mathbf{q}\left(-\mathit{\omega }\right)\mathrm{d}\mathit{\omega },\end{array}$$

which can be expressed more compactly as

$$\begin{array}{}\text{(21)}& \mathit{\beta }={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\mathbf{b}\left(\mathit{\omega }\right):\mathbf{D}\left(\mathit{\omega }\right)\mathrm{d}\mathit{\omega },\end{array}$$

where b(ω) is the tensor-valued encoding spectrum defined as (Topgaard, 2019a; Lundell and Lasič, 2020)

$$\begin{array}{}\text{(22)}& \mathbf{b}\left(\mathit{\omega }\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\mathbf{q}\left(\mathit{\omega }\right)\mathbf{q}{\left(-\mathit{\omega }\right)}^{\mathrm{T}},\end{array}$$

and “:” denotes a tensor dot product (Basser et al., 1994):

$$\begin{array}{}\text{(23)}& \mathbf{b}\left(\mathit{\omega }\right):\mathbf{D}\left(\mathit{\omega }\right)={\sum }_i{\sum }_j{b}_{ij}\left(\mathit{\omega }\right)D_{ij}\left(\mathit{\omega }\right).\end{array}$$

Combining Eqs. (11), (14), and (21) yields

$$\begin{array}{}\text{(24)}& E\left(\mathit{\tau }\right)=\exp \left(\mathrm{i}{\mathbf{q}}_{\mathrm{v}}\cdot 〈\mathbf{v}〉-{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\mathbf{b}\left(\mathit{\omega }\right):\mathbf{D}\left(\mathit{\omega }\right)\mathrm{d}\mathit{\omega }\right),\end{array}$$

where the motion-encoding properties of g(t) are summarized in q_v and b(ω), as illustrated in Fig. 1d and e. While q_v and 〈v〉 are (ω-independent) vectors, both b(ω) and D(ω) are symmetric second-order tensors at each value of ω.

2.5 Diffusion spectra for some simple cases

For a liquid with bulk diffusivity D₀ confined in d dimensions in planar (d=1), cylindrical (d=2), or spherical (d=3) compartments with radius r, the diffusion spectrum D_rest(ω) in the restricted dimensions can be expressed as follows (Stepišnik, 1993):

$$\begin{array}{}\text{(25)}& D_{\text{rest}}\left(\mathit{\omega }\right)=D_{\mathrm{0}}-\sum _k{w}_k{\displaystyle \frac{D_{\mathrm{0}}-D_{\mathrm{\infty }}}{\mathrm{1}+{\mathit{\omega }}^{\mathrm{2}}/{\mathrm{\Gamma }}_k^{\mathrm{2}}}},\end{array}$$

where

$$\begin{array}{}\text{(26)}& {\mathrm{\Gamma }}_k={\displaystyle \frac{{\mathit{\xi }}_k^{\mathrm{2}}{D}_{\mathrm{0}}}{r^{\mathrm{2}}}}\end{array}$$

and

$$\begin{array}{}\text{(27)}& w_k={\displaystyle \frac{\mathrm{2}}{{\mathit{\xi }}_k^{\mathrm{2}}+\mathrm{1}-d}}.\end{array}$$

Equation (25) includes the long-range diffusivity D_∞, allowing for finite permeability of the compartment walls (Lasič et al., 2009), and can be recognized as a sum of Lorentzians with widths Γ_k and weights w_k. In Eqs. (26) and (27), ξ_k is the kth solution of

$$\begin{array}{}\text{(28)}& \mathit{\xi }{J}_{d/\mathrm{2}-\mathrm{1}}\left(\mathit{\xi }\right)-\left(d-\mathrm{1}\right)J_{d/\mathrm{2}}\left(\mathit{\xi }\right)=\mathrm{0},\end{array}$$

where J_ν is the νth-order Bessel function of the first kind. Using the cylindrical case in Fig. 1b as an example, the tensor-valued diffusion spectrum D(ω) is given by

$$\begin{array}{}\text{(29)}& \mathbf{D}\left(\mathit{\omega }\right)=\mathbf{R}\left(\mathit{\theta },\mathit{\varphi }\right)\left(\begin{array}{ccc}{D}_{\text{rest}}\left(\mathit{\omega }\right)& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& D_{\text{rest}}\left(\mathit{\omega }\right)& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& D_{\mathrm{0}}\end{array}\right){\mathbf{R}}^{-\mathrm{1}}\left(\mathit{\theta },\mathit{\varphi }\right),\end{array}$$

where θ and ϕ are polar and azimuthal angles, giving the orientation of the cylinder in the lab frame, and R(θ, ϕ) is a rotation matrix. Figure 1d includes a plot of D(ω) for the case θ=0 and ϕ=0 where all off-diagonal elements are zero. At high values of ω, the diagonal elements converge towards D₀, corresponding to isotropic diffusion. Conversely, the effects of anisotropy reach a maximum in the low-ω limit where D_rest(ω) approaches D_∞, which equals zero in the example in Fig. 1d. The planar version of Eq. (29) is obtained by exchanging D_rest(ω) and D₀. In the low-ω limit, the planar and cylindrical cases are often combined into a single expression:

$$\begin{array}{}\text{(30)}& \mathbf{D}=\mathbf{R}\left(\mathit{\theta },\mathit{\varphi }\right)\left[\begin{array}{ccc}{D}_{\perp }& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& D_{\perp }& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& D_{\parallel }\end{array}\right]{\mathbf{R}}^{-\mathrm{1}}\left(\mathit{\theta },\mathit{\varphi }\right),\end{array}$$

where D_∥ and D_⊥ are the eigenvalues parallel and perpendicular to the main symmetry axis of the compartment, respectively. For completeness, we note that D_rest(ω) and D in Eqs. (25) and (30), respectively, reduce to an ω-independent scalar diffusion coefficient D for the special case of isotropic Gaussian diffusion where $D=D_{\mathrm{0}}=D_{\mathrm{\infty }}=D_{\perp }=D_{\parallel }$.

2.6 Key properties of the tensor-valued diffusion encoding spectrum

While the signal expression in Eq. (24) takes the ω dependence and tensorial properties of both b(ω) and D(ω) into account and may be numerically evaluated as a single matrix multiplication after discretization in the ω dimension and appropriate reordering of the tensor elements, the common occurrence of systems exhibiting approximately Gaussian (ω-independent) and/or isotropic diffusion has led to the introduction of simplified descriptions focusing on some specific aspects. In the absence of diffusion anisotropy – which is obviously not the case for the example in Fig. 1 – it is sufficient to use the dephasing power spectrum b(ω) (Stepišnik, 1981) obtained from b(ω) by

$$\begin{array}{}\text{(31)}& b\left(\mathit{\omega }\right)=\text{trace}\left\{\mathbf{b}\left(\mathit{\omega }\right)\right\}.\end{array}$$

The sensitivity to restriction can be summarized by the centroid frequency ω_cent (Arbabi et al., 2020), defined as

$$\begin{array}{}\text{(32)}& {\mathit{\omega }}_{\text{cent}}={\displaystyle \frac{\mathrm{1}}{b}}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\left|\mathit{\omega }\right|b\left(\mathit{\omega }\right)\mathrm{d}\mathit{\omega }.\end{array}$$

In addition to all the tensor elements of b(ω), Fig. 1d includes a plot of b(ω) with an arrow indicating ω_cent. The example of b(ω) covers both low- and high-ω features of D(ω) and is, thus, less well suited for exploring ω-dependent diffusion processes than gradient modulation schemes comprising trains of rectangular pulses (Callaghan and Stepišnik, 1995), multiple smooth oscillations (Parsons et al., 2003), or a Carr–Purcell–Meiboom–Gill sequence in the presence of a constant gradient (Lasič et al., 2006), where the encoding power is concentrated in a narrow frequency range and the single value ω_cent captures most of the relevant information about the spectral content. Even when the peaks in b(ω) are broader than the features in D(ω), the ω_cent metric has some value as a bookkeeping tool but is less suitable for quantitative analysis.

For anisotropic systems with Gaussian diffusion, it is useful to introduce the b-matrix b (Basser et al., 1994) by

$$\begin{array}{}\text{(33)}& \mathbf{b}={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\mathbf{b}\left(\mathit{\omega }\right)\mathrm{d}\mathit{\omega }.\end{array}$$

The sensitivity to anisotropy is given by the “shape” of the tensor which can be quantified with the encoding anisotropy b_Δ and asymmetry b_η (Eriksson et al., 2015), defined as

$$\begin{array}{}\text{(34)}& b_{\mathrm{\Delta }}={\displaystyle \frac{\mathrm{1}}{b}}\left(b_{ZZ}-{\displaystyle \frac{b_{YY}+b_{XX}}{\mathrm{2}}}\right)\end{array}$$

and

$$\begin{array}{}\text{(35)}& b_{\mathit{\eta }}={\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}{\displaystyle \frac{b_{YY}-b_{XX}}{b{b}_{\mathrm{\Delta }}}},\end{array}$$

respectively. Here, b_XX, b_YY, and b_ZZ are the eigenvalues of b ordered according to the Haeberlen convention $|b_{ZZ}-b/\mathrm{3}|\mathit{>}|b_{XX}-b/\mathrm{3}|\mathit{>}|b_{YY}-b/\mathrm{3}|$ (Haeberlen, 1976) and

$$\begin{array}{}\text{(36)}& b=\text{trace}\left\{\mathbf{b}\right\}\end{array}$$

is the conventional b value (Le Bihan et al., 1986) that gives the overall magnitude of the diffusion encoding. Figure 1e shows a superquadric tensor glyph (Kindlmann, 2004) representation of b where all eigenvalues and both shape parameters b_Δ and b_η are nonzero.

From the definitions of q_v, q(ω), and b(ω) in Eqs. (15), (18), and (22), it follows that the ω=0 values of b(ω) are proportional to ${\mathbf{q}}_{\mathrm{v}}{\mathbf{q}}_{\mathrm{v}}^{\mathrm{T}}$ and, thus, report on the sensitivity to flow, albeit with some ambiguity with respect to the directionality: the vectors q_v and −q_v give the same b(ω=0). It should be noted that q_v is not necessarily colinear with any of the eigenvectors of b.

2.7 Special cases for data analysis

For data-fitting purposes, it is convenient to write the normalized signal E as the ratio

$$\begin{array}{}\text{(37)}& E=S/S_{\mathrm{0}},\end{array}$$

where S is detected signal and S₀ is the signal obtained in a reference measurement with the amplitudes of the motion-encoding gradients set to zero. Equation (24) can then be expressed as

$$\begin{array}{}\text{(38)}& S=S_{\mathrm{0}}\exp \left(\mathrm{i}{\mathbf{q}}_{\mathrm{v}}\cdot 〈\mathbf{v}〉-{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\mathbf{b}\left(\mathit{\omega }\right):\mathbf{D}\left(\mathit{\omega }\right)\mathrm{d}\mathit{\omega }\right).\end{array}$$

For the special cases of (i) isotropic restricted, (ii) anisotropic Gaussian, and (iii) isotropic Gaussian diffusion in the absence of net flow (〈v〉=0), Eq. (38) is simplified to

$$\begin{array}{}\text{(39)}& {\displaystyle }S=S_{\mathrm{0}}\exp \left(-{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}b\left(\mathit{\omega }\right)D\left(\mathit{\omega }\right)\mathrm{d}\mathit{\omega }\right),\text{(40)}& {\displaystyle }S=S_{\mathrm{0}}\exp \left(-\mathbf{b}:\mathbf{D}\right),\end{array}$$

and

$$\begin{array}{}\text{(41)}& S=S_{\mathrm{0}}\exp \left(-bD\right),\end{array}$$

respectively, where D(ω) is the (isotropic) diffusion spectrum, D is the (ω-independent) diffusion tensor, and D is the (isotropic and ω-independent) diffusion coefficient introduced in Sect. 2.5 above.

For a heterogeneous system including multiple sub-ensembles i with individual signals S_i, each of which is given by one of the equations (Eqs. 38–41) above, the total signal S is obtained by the sum

$$\begin{array}{}\text{(42)}& S=\sum _i{S}_i.\end{array}$$

An important type of heterogeneity refers to the orientations of anisotropic objects with the extreme case of completely random orientations as in a “powder”. In the special case of axial symmetry of both b and D, powder averaging of Eq. (40) yields (Eriksson et al., 2015)

$$\begin{array}{}\text{(43)}& S=S_{\mathrm{0}}\exp \left(-b{D}_{\mathrm{iso}}\right){\displaystyle \frac{\sqrt{\mathit{\pi }}}{\mathrm{2}}}{\displaystyle \frac{\exp \left(A/\mathrm{3}\right)}{\sqrt{A}}}\mathrm{erf}\left(\sqrt{A}\right),\end{array}$$

where

$$\begin{array}{}\text{(44)}& A=\mathrm{3}b{D}_{\mathrm{iso}}{b}_{\mathrm{\Delta }}{D}_{\mathrm{\Delta }}.\end{array}$$

In Eq. (44), D_iso is the isotropic diffusivity and D_Δ is the normalized diffusion anisotropy; these terms are defined as

$$\begin{array}{}\text{(45)}& D_{\mathrm{iso}}={\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\text{trace}\left\{\mathbf{D}\right\}={\displaystyle \frac{D_{\parallel }+\mathrm{2}{D}_{\perp }}{\mathrm{3}}}\end{array}$$

and

$$\begin{array}{}\text{(46)}& D_{\mathrm{\Delta }}={\displaystyle \frac{D_{\parallel }-D_{\perp }}{\mathrm{3}{D}_{\mathrm{iso}}}},\end{array}$$

respectively. Here, D_∥ and D_⊥ were introduced in Eq. (30). The definitions of b_Δ and b can be found in Eqs. (34) and (36), respectively.

3 Design of gradient waveforms by double rotation of the q vector

Expanding on previous magic-angle spinning (Andrew et al., 1959; Eriksson et al., 2013; Topgaard, 2013) and variable-angle spinning (Frydman et al., 1992; Topgaard, 2016, 2017) approaches for generating motion-encoding gradient waveforms, we apply the double-rotation (DOR) technique (Samoson et al., 1998; Topgaard, 2019a) to probe the 2D acquisition space spanned by the variables ω_cent and b_Δ defined in Eqs. (32) and (34), respectively. The q-vector trajectory q(t) is expressed in terms of its time-dependent magnitude q(t) and unit vector u(t) as

$$\begin{array}{}\text{(47)}& \mathbf{q}\left(t\right)=q\left(t\right)\mathbf{u}\left(t\right).\end{array}$$

For the special case of DOR, the unit vector is written as

$$\begin{array}{}\text{(48)}& \mathbf{u}\left(t\right)={\mathbf{R}}_z\left({\mathit{\psi }}_{\mathrm{2}}\left(t\right)\right){\mathbf{R}}_y\left({\mathit{\zeta }}_{\mathrm{2}}\right){\mathbf{R}}_z\left({\mathit{\psi }}_{\mathrm{1}}\left(t\right)\right){\mathbf{R}}_y\left({\mathit{\zeta }}_{\mathrm{1}}\right){\left[\begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]}^{\mathrm{T}},\end{array}$$

where R_z and R_y are Euler rotation matrices, ζ₁ and ζ₂ are the inclinations of the two rotation axes, and ψ₁(t) and ψ₂(t) are the time-dependent angles of rotation. The rotations in Eq. (48) are applied from right to left and follow a Z–Y active rotation matrix convention.

Starting from a conventional 1D gradient waveform g_1D(t) – for instance, a pair of rectangular or sine-bell pulses of opposite polarity – the time-dependent functions q(t) and ψ₂(t) are given by Topgaard (2016):

$$\begin{array}{}\text{(49)}& q\left(t\right)=\mathit{\gamma }{\int }_{\mathrm{0}}^t{g}_{\mathrm{1}\mathrm{D}}\left(t^{\prime }\right)\mathrm{d}{t}^{\prime }\end{array}$$

and

$$\begin{array}{}\text{(50)}& {\mathit{\psi }}_{\mathrm{2}}\left(t\right)={\displaystyle \frac{\mathrm{\Delta }{\mathit{\psi }}_{\mathrm{2}}}{b}}{\int }_{\mathrm{0}}^t{q}^{\mathrm{2}}\left(t{}^{\prime }\right)\mathrm{d}t{}^{\prime },\end{array}$$

where Δψ₂ is the total angle of rotation during the encoding interval from time t=0 to τ and

$$\begin{array}{}\text{(51)}& b={\int }_{\mathrm{0}}^{\mathit{\tau }}{q}^{\mathrm{2}}\left(t\right)\mathrm{d}t\end{array}$$

is the conventional b value. After some exercises in trigonometry, combination of Eqs. (47)–(51) and the relation between g(t) and q(t) in Eq. (4) yields

$$\begin{array}{}\text{(52)}& \begin{array}{rl}{\mathbf{g}}_{\mathrm{DOR}}\left(t\right)& =g_{\mathrm{1}\mathrm{D}}\left(t\right)\left(\begin{array}{c}{a}_{+}\cos {\mathit{\psi }}_{+}\left(t\right)+a_{-}\cos {\mathit{\psi }}_{-}\left(t\right)\\ +a_{\mathrm{2}}\cos {\mathit{\psi }}_{\mathrm{2}}\left(t\right)\\ a_{+}\sin {\mathit{\psi }}_{+}\left(t\right)-a_{-}\sin {\mathit{\psi }}_{-}\left(t\right)\\ +a_{\mathrm{2}}\sin {\mathit{\psi }}_{\mathrm{2}}\left(t\right)\\ a_{\mathrm{0}}-a_{\mathrm{1}}\cos {\mathit{\psi }}_{\mathrm{1}}\left(t\right)\end{array}\right)\\ & +g_{\mathrm{rot}}\left(t\right)\left(\begin{array}{c}-\left(n+\mathrm{1}\right)a_{+}\sin {\mathit{\psi }}_{+}\left(t\right)\\ -\left(n-\mathrm{1}\right)a_{-}\sin {\mathit{\psi }}_{-}\left(t\right)\\ -a_{\mathrm{2}}\sin {\mathit{\psi }}_{\mathrm{2}}\left(t\right)\\ \left(n+\mathrm{1}\right)a_{+}\cos {\mathit{\psi }}_{+}\left(t\right)\\ -\left(n-\mathrm{1}\right)a_{-}\cos {\mathit{\psi }}_{-}\left(t\right)\\ +a_{\mathrm{2}}\cos {\mathit{\psi }}_{\mathrm{2}}\left(t\right)\\ n{a}_{\mathrm{1}}\sin {\mathit{\psi }}_{\mathrm{1}}\left(t\right)\end{array}\right),\end{array}\end{array}$$

where

$$\begin{array}{}\text{(53)}& g_{\mathrm{rot}}\left(t\right)={\displaystyle \frac{\mathrm{\Delta }{\mathit{\psi }}_{\mathrm{2}}q{\left(t\right)}^{\mathrm{3}}}{\mathit{\gamma }b}}\end{array}$$

is the time-dependent magnitude of the rotating gradient vector;

$$\begin{array}{}\text{(54)}& \begin{array}{rl}{\mathit{\psi }}_{\mathrm{1}}\left(t\right)& =n{\mathit{\psi }}_{\mathrm{2}}\left(t\right)\text{and}\\ {\mathit{\psi }}_{\pm }\left(t\right)& =\left(n\pm \mathrm{1}\right){\mathit{\psi }}_{\mathrm{2}}\left(t\right)\end{array}\end{array}$$

are time-dependent rotation angles; and

$$\begin{array}{}\text{(55)}& \begin{array}{rl}{a}_{\mathrm{0}}& =\cos {\mathit{\zeta }}_{\mathrm{1}}\cos {\mathit{\zeta }}_{\mathrm{2}},\\ a_{\mathrm{1}}& =\sin {\mathit{\zeta }}_{\mathrm{1}}\sin {\mathit{\zeta }}_{\mathrm{2}},\\ a_{\mathrm{2}}& =\cos {\mathit{\zeta }}_{\mathrm{1}}\sin {\mathit{\zeta }}_{\mathrm{2}},\text{and}\\ a_{\pm }& =\sin {\mathit{\zeta }}_{\mathrm{1}}{\displaystyle \frac{\cos {\mathit{\zeta }}_{\mathrm{2}}\pm \mathrm{1}}{\mathrm{2}}}\end{array}\end{array}$$

are amplitudes of the oscillating terms.

In the solid-state NMR field, DOR is applied with the inclinations ${\mathit{\zeta }}_{\mathrm{1}}=\mathrm{54.7}{}^{\circ }$ and ${\mathit{\zeta }}_{\mathrm{2}}=\mathrm{30.6}{}^{\circ }$, corresponding to zeros of the second and fourth Legendre polynomial, to eliminate both first- and second-order quadrupolar broadening of the NMR spectra of nuclei such as ²³Na (Samoson et al., 1998). While the encoding anisotropy b_Δ in Eq. (34) is closely related to the second Legendre polynomial (Eriksson et al., 2015), we are not yet aware of any diffusion analog of the quadrupolar interactions involving the fourth Legendre polynomial. Instead, we found that DOR with the inclinations ${\mathit{\zeta }}_{\mathrm{1}}=\mathrm{90}{}^{\circ }$ and ${\mathit{\zeta }}_{\mathrm{2}}=-\mathrm{54.7}{}^{\circ }$ yields desirable properties for diffusion encoding, namely spectral content concentrated to a narrow frequency window for all of the elements of b(ω). For n>1 and the special case of g_1D(t) ∝ [δ(t) − δ(t − τ)], where δ(x) is the Dirac delta function, these angles yield an isotropic b tensor, corresponding to b_Δ=0. Waveforms for any values of b_Δ and b_η are then conveniently obtained by scaling the components of g_DOR(t) according to

$$\begin{array}{}\text{(56)}& \mathbf{g}\left(t\right)=\left[\begin{array}{c}{g}_X\left(t\right)\\ g_Y\left(t\right)\\ g_Z\left(t\right)\end{array}\right]=\left[\begin{array}{c}{g}_{\mathrm{DOR},X}\left(t\right)\sqrt{\mathrm{1}-b_{\mathrm{\Delta }}\left(\mathrm{1}+b_{\mathit{\eta }}\right)}\\ g_{\mathrm{DOR},Y}\left(t\right)\sqrt{\mathrm{1}-b_{\mathrm{\Delta }}\left(\mathrm{1}-b_{\mathit{\eta }}\right)}\\ g_{\mathrm{DOR},Z}\left(t\right)\sqrt{\mathrm{1}+\mathrm{2}{b}_{\mathrm{\Delta }}}\end{array}\right].\end{array}$$

Figure 2 Flowchart for calculating a double-rotation gradient waveform g_DOR(t) given a 1D dephasing/rephasing waveform g_1D(t), rotation axis inclinations ζ₁ and ζ₂, double-rotation ratio n, and total angle of rotation Δψ₂ during the waveform duration τ. The waveform g_1D(t), containing dephasing and rephasing pulses with sinusoidal ramps of durations ε_up and ε_down, gives the time-dependent magnitude of the dephasing vector q(t) via Eq. (49), which yields the time-dependent rotation angles ψ_i(t) via Eqs. (50) and (54) as well as the time-dependent magnitude of the rotating and oscillating gradient g_rot(t) via Eq. (53). Combining g_1D(t), g_rot(t), and ψ_i(t) via Eq. (52) gives g_DOR(t), which, if ${\mathit{\zeta }}_{\mathrm{1}}=\mathrm{90}{}^{\circ }$ and ${\mathit{\zeta }}_{\mathrm{2}}=-\mathrm{54.7}{}^{\circ }$, n is an integer above 1, and Δψ₂ is a multiple of 360^∘, achieves isotropic encoding tensors b where the anisotropy b_Δ and asymmetry b_η are both equal to zero. Finally, waveforms g(t) for any values of b_Δ and b_η are obtained by scaling of the Cartesian components of g_DOR(t) according to Eq. (56). The shown example was generated with the MATLAB code provided in the Supplement using ε_up=0.015τ, ε_down=0.06τ, $\mathrm{\Delta }{\mathit{\psi }}_{\mathrm{2}}=\mathrm{360}{}^{\circ }$, n=4, b_Δ=0.5, and b_η=0.25.

At the selected inclinations, the a₀ and a₂ terms in Eq. (52) equal zero, while the remaining amplitudes evaluate to $a_{\mathrm{1}}\approx -\mathrm{0.816}$, ${\mathit{\text{a}}}_{+}\approx \mathrm{0.789}$, and $a_{-}\approx -\mathrm{0.211}$. For the special case g_1D(t) ∝ [δ(t) − δ(t − τ)], the main frequency components of b(ω) are thus given by

$$\begin{array}{}\text{(57)}& \begin{array}{rl}{\mathit{\omega }}_{\pm }& ={\displaystyle \frac{{\mathit{\psi }}_{\pm }\left(\mathit{\tau }\right)}{\mathit{\tau }}}=\left(n\pm \mathrm{1}\right){\displaystyle \frac{\mathrm{\Delta }{\mathit{\psi }}_{\mathrm{2}}}{\mathit{\tau }}}\mathrm{and}\\ {\mathit{\omega }}_{\mathrm{1}}& ={\displaystyle \frac{{\mathit{\psi }}_{\mathrm{1}}\left(\mathit{\tau }\right)}{\mathit{\tau }}}=n{\displaystyle \frac{\mathrm{\Delta }{\mathit{\psi }}_{\mathrm{2}}}{\mathit{\tau }}},\end{array}\end{array}$$

where, according to Eq. (52), ω_± and ω₁ are cleanly separated into the respective transverse (X, Y) and longitudinal (Z) directions. The mean frequency content, as quantified by the centroid frequency ω_cent defined in Eq. (32), can be estimated by weighting the contributions from the main frequency components by the corresponding amplitudes in Eq. (52) but is more accurately calculated by the numerical evaluation of Eq. (32), which also takes the finite durations of the sinusoidal oscillations into account. For rough prediction of ω_cent, it is useful to note that $a_{+}^{\mathrm{2}}\gg a_{-}^{\mathrm{2}}$, implying that the ω₊ component will dominate the spectra in the X and Y directions. The scaling of the waveforms according to Eq. (56) preserves the frequency content in each of the eigendirections of the b tensor but shifts the value of ω_cent between the approximate extremes ω₊ and ω₁ for $b_{\mathrm{\Delta }}=-\mathrm{1}/\mathrm{2}$ and 1, respectively. The differences in ω_± and ω₁ may give rise to a directional dependence of the sensitivity to restriction as investigated for the b_Δ=0 case by de Swiet and Mitra (1996).

Figure 2 illustrates the series of calculations required to convert a conventional 1D waveform g_1D(t) and given values of ζ₁, ζ₂, Δψ₂, n, b_Δ, and b_η to a 3D waveform g(t) by numerical evaluation of Eqs. (49)–(56). Following previous works to generate families of smooth gradient waveforms to explore the b_Δ and b_η dimensions of diffusion encoding (Topgaard, 2016, 2017), we construct g_1D(t) from a dephasing lobe with quarter-sine ramp-up of duration ε_up and half-cosine ramp-down of duration ε_down as well as a rephasing lobe obtained by inversion and time-reversal of the dephasing one. The corresponding MATLAB code is provided in the Supplement.

Figure 3 compiles waveforms and encoding spectra for an array of n and b_Δ at constant g_1D(t), τ, and Δψ₂, yielding constant b. Increasing n leads to larger rotation angles ψ₁(t) and ψ_±(t) and frequencies ω₁ and ω_± according to Eqs. (54) and (57), respectively, at the expense of overall higher gradient amplitudes on account of the terms including n in Eq. (52). For most waveforms, vanishing values of b(ω) at ω=0 correspond to q_v=0 and insensitivity to flow. Many of the examples in Fig. 3 are familiar from the literature, such as conventional Stejskal–Tanner encoding at (n=0, b_Δ=1), basic flow-compensated encoding (Caprihan and Fukushima, 1990) at (n=1, b_Δ=1), and magic-angle spinning of the q vector (Eriksson et al., 2013) at (n=0, b_Δ=0). The series of b_Δ=1 and $-\mathrm{1}/\mathrm{2}$ waveforms with varying n resemble the cosine-modulated oscillating gradients of Parsons et al. (2003) and the circularly polarized version introduced by Lundell et al. (2015), respectively. Correspondingly, the series of waveforms with n=0 and varying b_Δ has previously been introduced as a diffusion version of the variable-angle spinning technique to correlate isotropic and anisotropic chemical shifts in solid-state NMR (Topgaard, 2016, 2017). The approach for joint investigation of restricted and anisotropic diffusion proposed by Lundell et al. (2019), combining isotropic encoding with “tuned” and “detuned” directional encodings, can be recognized as measurements at the three discrete points (n=0, b_Δ=0), (n=0, b_Δ=1), and (n=1, b_Δ=1) of the 2D plane in Fig. 3. For completeness, we note that the elliptically polarized oscillating gradients by Nielsen et al. (2018) can be reproduced with $b_{\mathrm{\Delta }}=-\mathrm{1}/\mathrm{2}$ and b_η in the range from 0 to 3 (not included in Fig. 3).

4 Proof-of-principle experiments

Magnesium nitrate hexahydrate, cobalt nitrate hexahydrate, and 1-decanol were purchased from Sigma-Aldrich Sweden AB and sodium octanoate was obtained from J&K Scientific via Th. Geyer in Sweden. Water was purified with a Milli-Q system. A sample with two-component isotropic diffusion was prepared by inserting a 4 mm NMR tube containing an aqueous solution saturated with magnesium nitrate (Wadsö et al., 2009) into a 10 mm NMR tube with water (Mills, 1973). The magnesium nitrate solution was spiked with a small amount of cobalt nitrate (0.27 wt % saturated solution) to reduce T₁ and T₂ to approx. 500 and 50 ms, respectively. An anisotropic sample was prepared by mixing 85.79 wt % Milli-Q, 9.17 wt % 1-decanol, and 5.04 wt% sodium octanoate giving a lamellar liquid crystal (Persson et al., 1975). Investigation of isotropic restricted diffusion was performed with a sediment of fresh baker's yeast (trade name: Kronjäst; obtained from a local supermarket) prepared by dispersing yeast in tap water (1:1 volume ratio) in a glass vial, transferring it with a syringe to a 10 mm NMR tube to a sample height of 40 mm, and keeping the tube in an upright position overnight at 4 ^∘C to allow the cells to settle under the force of gravity into a 20 mm high pellet (Malmborg et al., 2006).

Figure 3 Gradient waveforms g(t) for comprehensive exploration of the 2D space of centroid frequency ω_cent and anisotropy b_Δ of the tensor-valued encoding spectrum b(ω). Magic-angle spinning (MAS) and double rotation (DORn) with variable frequency ratio n give q-vector trajectories q(t) shown as 3D plots for the b_Δ=0 cases. The waveforms are generated according to the scheme in Fig. 2 using ε_up=0.03τ, ε_down=0.12τ, $\mathrm{\Delta }{\mathit{\psi }}_{\mathrm{2}}=\mathrm{360}{}^{\circ }$, b_η=0, and identical b values for a 2D array of n=0, 1, …, 5 and $b_{\mathrm{\Delta }}=-\mathrm{0.5}$, 0, 0.5, and 1 with the angles ${\mathit{\zeta }}_{\mathrm{1}}=\mathrm{0}{}^{\circ }$ and ${\mathit{\zeta }}_{\mathrm{2}}=\mathrm{54.7}{}^{\circ }$ for n=0 and ${\mathit{\zeta }}_{\mathrm{1}}=\mathrm{90}{}^{\circ }$ and ${\mathit{\zeta }}_{\mathrm{2}}=-\mathrm{54.7}{}^{\circ }$ for n>0. Superquadric tensor glyphs (Kindlmann, 2004) along the vertical axis illustrate b for the chosen values of b_Δ. The main maxima in b(ω) are located at the frequencies ω₁ and ω₊ given in Eq. (57). Values of ω_cent and b_Δ, including non-idealities originating from the finite durations of the dephasing and rephasing lobes of g_1D(t), are obtained via Eqs. (32) and (34), respectively, using numerically evaluated b(ω) according to Eqs. (4), (18), and (22).

MRI was performed on a Bruker AVANCE NEO 500 MHz spectrometer equipped with an 11.7 T vertical bore magnet and a MIC-5 microimaging probe fitted with a 10 mm radiofrequency insert for observation of ¹H. Images were acquired with a TopSpin 4.0 implementation of a spin-echo prepared single-shot RARE (Rapid Imaging with Refocused Echoes) sequence (available at https://github.com/daniel-topgaard/md-dmri, last access: 1 October 2022) using a 0.6 × 0.6 mm² resolution in a plane perpendicular to the tube axis, a 1 mm slice thickness, and a 16 × 16 × 1 matrix size. Diffusion encoding employed pairs of identical gradient waveforms bracketing the 180^∘ pulse in the preparation block (Lasič et al., 2014). Data were acquired for 8 b values up to 6.44⋅10⁹ s m⁻² and 15 orientations (Θ, Φ) for each of the 24 waveforms spanning the 2D ω_cent–b_Δ plane in Fig. 3 using a maximum gradient amplitude of 3 T m⁻¹ and a waveform duration of τ=25 ms, giving values of ω_cent in the range of 20–260 Hz. With a 5 s recycle delay, the total measurement time was approximately 4 h for each sample. The sample temperature was controlled with a Bruker VT unit: 278 K for the yeast and 291 K for the isotropic solutions and liquid crystal. For the yeast, the image slice was placed in the middle of the pellet and 10 mm below the bottom of the supernatant. Image reconstruction, definition of regions of interest, and curve fitting were performed in MATLAB using in-house code available at https://github.com/daniel-topgaard/md-dmri, last access: 1 October 2022 (Nilsson et al., 2018).

Figure 4 Experimental (markers) and fitted (lines) normalized powder-averaged signal $S/S_{\mathrm{0}}$ vs. b value for phantoms with well-defined diffusion properties. (a) Tube-in-tube assembly of pure water and a concentrated solution of magnesium nitrate in water (“brine”) giving rise to two isotropic Gaussian (ω-independent) components. A two-component fit based on Eq. (41) gives diffusion spectra D(ω), as shown in the inset, with percentages indicating the relative contributions. The legend shows the 24 investigated values in the 2D ω_cent (the gray scale) and b_Δ (marker style) space corresponding to the gradient waveforms in Fig. 3 with a duration of τ=25 ms, a maximum gradient strength of 3 Tm⁻¹, and pairs of waveforms bracketing the 180^∘ pulse in the spin-echo preparation. (b) Polydomain lamellar liquid crystal giving Gaussian parallel and perpendicular diffusivities, D_∥ and D_⊥, as estimated by a fit of Eq. (43). (c) Sediment of yeast cells with intra- and extracellular compartments, with the former exhibiting restricted (ω-dependent) diffusion. The inset shows D(ω) resulting from a two-component fit with one spherically restricted (solid) and one Gaussian (dashed) component. For the former component, the signal was obtained by numerical integration of Eq. (39) with D(ω) given by Eq. (25) with d=3 and D_∞ constrained to zero.

Figure 4 compiles experimental data and fits for all investigated samples. To facilitate visual inspection of the highly multidimensional data acquired as a function of (b, ω_cent, b_Δ, Θ, Φ), the signal data were averaged over b-tensor orientations (Θ, Φ) and are displayed as conventional Stejskal–Tanner plots of log₁₀(S) vs. b with the ω_cent and b_Δ dimensions coded using different marker styles and a gray scale. For the isotropic Gaussian sample in Fig. 4a, all data points collapse onto a single master curve, thereby verifying that all 24 waveforms spanning the 2D ω_cent–b_Δ plane in Fig. 3 indeed give the same b value. The pronounced nonlinearity of the log₁₀(S) vs. b plot indicates the presence of multiple species with different diffusivities, and the bi-exponential fit yields diffusivities consistent with pure water (fast) and water in the saturated magnesium nitrate solution (slow). The anisotropic Gaussian phantom Fig. 4b yields data points stratified into one master curve for each of the four values of b_Δ, verifying independence of ω_cent. The data are well fitted by the expression for randomly oriented axisymmetric diffusion tensors in Eq. (43), giving estimates of the diffusivities D_∥ and D_⊥ parallel and perpendicular to the cylindrical symmetry axis of the crystallites. The observations $D_{\parallel }\ll D_{\perp }$ and $D_{\parallel }\approx \mathrm{0}$ are consistent with diffusion in a lamellar liquid crystal with planar surfactant bilayers being nearly impermeable to water (Callaghan and Söderman, 1983). For the isotropic restriction phantom in Fig. 4c, the signal depends strongly on ω_cent. In this case, there is no clear stratification of data points into separate master curves on account of the interplay between b_Δ and ω_cent, as reported in the legend in Fig. 4a and explained in detail below Eq. (57). The minor dependence of ω_cent on b_Δ is admittedly a drawback of our current approach for generating waveforms; however, we believe our approach is justified by the simplicity and transparency of the mathematical expressions in Eqs. (49)–(57). The sum of isotropic restricted and Gaussian components, given by Eqs. (41) and (39), yields an excellent fit to the acquired data, showing that the data feature no dependence on the value of b_Δ. To account for the fact that b(ω) for (in particular) the b_Δ=0 waveforms cannot be well approximated with a delta function at a single value of ω, the signal for the restricted component was obtained by numerical evaluation of the integral in Eq. (39) using the diffusion spectrum D(ω) for spherical compartments in Eq. (25). The obtained diffusivities are consistent with previous results for intra- and extracellular water in yeast cell sediments (Åslund and Topgaard, 2009). Taken together, the data in Fig. 4 verify that the set of waveforms allows detailed exploration of the 2D ω_cent–b_Δ plane of multidimensional diffusion encoding.

5 Conclusions and outlook

The proposed family of double-rotation gradient waveforms enables comprehensive sampling of both the frequency and “shape” dimensions of diffusion encoding, as required for detailed characterization of restriction and anisotropy in heterogeneous materials such as brain tissues. The present waveforms, deriving from simple geometrical considerations and generated by compact mathematical expressions, are suitable for preclinical investigations of tissue samples or small animals on high-gradient systems. By numerical optimizations to maximize the b value for given gradient strength (Topgaard, 2013; Sjölund et al., 2015), mitigating image artifacts from eddy currents (Yang and McNab, 2019) and concomitant gradients (Szczepankiewicz et al., 2019), and further minimizing side lobes in the encoding spectra (Hennel et al., 2020), we anticipate that the waveforms may be adapted for human in vivo studies. The merging of oscillating gradients (Aggarwal, 2020) and tensor-valued encoding (Reymbaut, 2020) into a common acquisition protocol encourages further development of a joint analysis framework, for instance, by augmenting current nonparametric diffusion tensor distributions (Topgaard, 2019b) with a Lorentzian frequency dimension (Narvaez et al., 2021, 2022) or building on the concept of confinement tensors (Yolcu et al., 2016; Boito et al., 2022).