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

Abstract 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.


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 frequencydependent (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 solidstate 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.
Published by Copernicus Publications on behalf of the Groupement AMPERE.
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).

Theoretical background
To set the stage for later sections dealing with the proposed gradient waveforms to investigate both frequencydependent 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.

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 where . . .denotes an ensemble mean and the timedependent 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 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 where q(t) is the dephasing vector, defined as 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.Focusing on translational displacements rather than absolute positions, we select a time τ where q (τ ) = 0 (5) and the first term in Eq. (3) vanishes while the second one remains: 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 τ .

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[φ(τ )] = ( φ(τ ) 2 -φ(τ ) 2 ) 1/2 .The Gaussian function can be expressed as where 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 q v /τ , 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).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). and After rewriting Eq. ( 1) as an integral, the insertion of Eq. ( 7) into Eq.( 10) may be evaluated to 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 https://doi.org/10.5194/mr-4-73-2023Magn. Reson., 4, 73-85, 2023 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).

Mean velocity and velocity correlation function
Insertion of Eq. ( 6) into Eq.( 8) yields which by separating v(t) into the ensemble mean v(t) = v and fluctuating part u(t), defined by can be evaluated to where the flow encoding vector q v is defined as Correspondingly, the insertion of Eq. ( 6) into Eq.( 9) gives which by reordering the time integrals and ensemble means as well as noting that u(t) • v = 0, can be expressed as where u(t)u(t ) T is the tensor-valued velocity correlation function.

Transformation to the frequency domain
After introducing the dephasing spectrum q(ω) and diffusion spectrum D(ω) by Fourier transformations to the frequency (ω) domain according to and Eq. ( 17) can be recast into which can be expressed more compactly as where b(ω) is the tensor-valued encoding spectrum defined as (Topgaard, 2019a;Lundell and Lasič, 2020) and ":" denotes a tensor dot product (Basser et al., 1994): Combining Eqs. ( 11), ( 14), and ( 21) yields 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 ω.

Diffusion spectra for some simple cases
For a liquid with bulk diffusivity D 0 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): where and 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 where J ν is the νth-order Bessel function of the first kind.Using the cylindrical case in Fig. 1b as an example, the tensorvalued diffusion spectrum D(ω) is given by 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 0 , 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 0 .In the low-ω limit, the planar and cylindrical cases are often combined into a single expression: 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

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 The sensitivity to restriction can be summarized by the centroid frequency ω cent (Arbabi et al., 2020), defined as 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) 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 and  , 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 q v q T v 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.

Special cases for data analysis
For data-fitting purposes, it is convenient to write the normalized signal E as the ratio where S is detected signal and S 0 is the signal obtained in a reference measurement with the amplitudes of the motionencoding gradients set to zero.Equation ( 24) can then be expressed as 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 and 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 subensembles 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 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) where In Eq. ( 44), D iso is the isotropic and D is the normalized diffusion anisotropy; these terms are defined as and 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.

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 variableangle spinning (Frydman et al., 1992;Topgaard, 2016Topgaard, , 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 For the special case of DOR, the unit vector is written as where R z and R y are Euler rotation matrices, ζ 1 and ζ 2 are the inclinations of the two rotation axes, and ψ 1 (t) and ψ 2 (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 ψ 2 (t) are given by Topgaard (2016): and where ψ 2 is the total angle of rotation during the encoding interval from time t = 0 to τ and 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 where is the time-dependent magnitude of the rotating gradient vector; ψ 1 (t) = nψ 2 (t) and At the selected inclinations, the a 0 and a 2 terms in Eq. ( 52) equal zero, while the remaining amplitudes evaluate to a 1 ≈ −0.816, a + ≈ 0.789, and a − ≈ −0.211.For the special case g 1D (t) ∝ [δ(t) − δ(t − τ )], the main frequency components of b(ω) are thus given by ψ 2 τ and where, according to Eq. ( 52), ω ± and ω 1 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 The waveform g 1D (t), containing dephasing and rephasing pulses with sinusoidal ramps of durations ε up and ε down , gives the timedependent 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 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 2 + a 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 ω 1 for b = −1/2 and 1, respectively.The differences in ω ± and ω 1 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 ζ 1 , ζ 2 , ψ 2 , 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(Topgaard, , 2017)), we construct g 1D (t) from a dephasing lobe with quarter-sine ramp-up of duration ε up and halfhttps://doi.org/10.5194/mr-4-73-2023Magn. Reson., 4, 73-85, 2023 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 ψ 2 , yielding constant b.Increasing n leads to larger rotation angles ψ 1 (t) and ψ ± (t) and frequencies ω 1 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 −1/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 solidstate NMR (Topgaard, 2016(Topgaard, , 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 = −1/2 and b η in the range from 0 to 3 (not included in Fig. 3).

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 1 and T 2 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).
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 1 H. Images were acquired with a TopSpin 4.0 implementation of a spinecho 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 2 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 9 s m −2 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 −1 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 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 10 (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 10 (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 sym-  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).metry axis of the crystallites.The observations D D ⊥ and D ≈ 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 be-lieve 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  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.

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(Narvaez et al., , 2022) ) or building on the concept of confinement tensors (Yolcu et al., 2016;Boito et al., 2022).
Figure1.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 q v /τ , 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).
respectively.Here, b XX , b Y Y , and b ZZ are the eigenvalues of b ordered according to the Haeberlen convention |b ZZ − b/3|>|b XX − b/3|>|b Y Y − b/3| (Haeberlen, 1976) and b = trace {b} (36) is the conventional b value (Le Bihan et al. ) https://doi.org/10.5194/mr-4-73-2023Magn.Reson., 4, 73-85, 2023 ) are time-dependent rotation angles; and a 0 = cos ζ 1 cos ζ 2 , a 1 = sin ζ 1 sin ζ 2 , a 2 = cos ζ 1 sin ζ 2 , anda ± = sin ζ 1 cos ζ 2 ± 1 2 (55)are amplitudes of the oscillating terms.In the solid-state NMR field, DOR is applied with the inclinations ζ 1 = 54.7 • and ζ 2 = 30.6• , 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 23 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 ζ 1 = 90 • and ζ 2 = −54.7 • 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

Figure 2 .
Figure2.Flowchart for calculating a double-rotation gradient waveform g DOR (t) given a 1D dephasing/rephasing waveform g 1D (t), rotation axis inclinations ζ 1 and ζ 2 , double-rotation ratio n, and total angle of rotation ψ 2 during the waveform duration τ .The waveform g 1D (t), containing dephasing and rephasing pulses with sinusoidal ramps of durations ε up and ε down , gives the timedependent 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 ζ 1 = 90 • and ζ 2 = −54.7 • , n is an integer above 1, and ψ 2 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τ , ψ 2 = 360 • , n = 4, b = 0.5, and b η = 0.25.
Figure2.Flowchart for calculating a double-rotation gradient waveform g DOR (t) given a 1D dephasing/rephasing waveform g 1D (t), rotation axis inclinations ζ 1 and ζ 2 , double-rotation ratio n, and total angle of rotation ψ 2 during the waveform duration τ .The waveform g 1D (t), containing dephasing and rephasing pulses with sinusoidal ramps of durations ε up and ε down , gives the timedependent 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 ζ 1 = 90 • and ζ 2 = −54.7 • , n is an integer above 1, and ψ 2 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τ , ψ 2 = 360 • , n = 4, b = 0.5, and b η = 0.25.

Figure 4 .
Figure 4. Experimental (markers) and fitted (lines) normalized powder-averaged signal S/S 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 −1 , 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.