The solid effect of dynamic nuclear polarization in liquids

The solid-state effect of dynamic nuclear polarization (DNP) is operative also in viscous liquids where the dipolar interaction between the electronic and nuclear spins is partially averaged. The proper way to quantify the degree of averaging, and thus calculate the efficiency of the effect, should be based on the time-correlation function of the dipolar interaction. Here we use the stochastic Liouville equation formalism to develop a general theoretical description of the solid effect in liquids. The derived expressions can be used with different dipolar correlations functions depending on the assumed motional model. At high magnetic fields, the theory predicts DNP enhancements at small offsets, far from the classical solid-effect positions that are displaced by one nuclear Larmor frequency from the electronic resonance. The predictions are in quantitative agreement with such enhancement peaks observed at 9.4 T (Kuzhelev et al., 2022). These non-canonical peaks are not due to thermal mixing or the cross effect but exactly follow the dispersive component of the EPR line.

1 Introduction

The last 2 decades have witnessed an overarching development in nuclear hyperpolarization techniques across the entire spectrum of mechanisms, from the classical Overhauser and solid-state effects to (photo)chemically induced and parahydrogen-based polarization (Eills et al., 2023). While the majority of the reported applications have relied on polarization transfer in the solid state at cryogenic temperatures (Ni et al., 2013; Pinon et al., 2021), transfer in the liquid state at elevated (room or physiological) temperatures has also been actively explored (Prisner et al., 2016; Denysenkov et al., 2022). In the liquid state, the mechanism of polarization has almost exclusively been the Overhauser effect (Jakdetchai et al., 2014). One notable exception is the work of Stapf and coworkers in which the solid effect of dynamic nuclear polarization (DNP) has been employed in combination with field-cycling relaxometry to characterize the molecular dynamics in ionic liquids and polymer melts at ambient temperatures (Neudert et al., 2017; Gizatullin et al., 2019, 2021a, b, 2022).

At X band (9.6 GHz/0.35 T), where the DNP measurements of Stapf and colleagues have been carried out, the nuclear Larmor frequencies of ¹H and ¹⁹F (ω_I≈15 MHz) are less than the EPR spectral width of a nitroxide free radical and even comparable to the spectral width of the single-line radical BDPA (Gizatullin et al., 2021b). As a result, the negative and positive solid-effect enhancements overlap and partly cancel each other, complicating the quantitative analysis of the effect. An additional difficulty for quantification is that, in many instances, the Overhauser and solid effects coexist (Leblond et al., 1971b; Neudert et al., 2017; Gizatullin et al., 2022). Although the contributions of these two effects can generally be distinguished on the basis of their even (Overhauser effect) and odd (solid effect) parity with the offset from the electronic resonance, this identification could be complicated when the EPR spectrum is broad and asymmetrical.

To quantify the field profile of the DNP enhancement (i.e., the DNP spectrum), Stapf and colleagues use a weighted sum of (i) the EPR line shape (for the Overhauser effect) and (ii) the same line shape shifted by ±ω_I, with one of the shifted copies flipped around the vertical axis (for the solid effect) (Neudert et al., 2017; Gizatullin et al., 2021b). The relative contribution of the two effects is then treated as a fitting parameter. In general, the resulting fits are in good overall agreement with the experimental DNP spectra, but oftentimes there are quantitative deviations. Recently, the remaining discrepancy between the experimental and calculated DNP spectra was interpreted as evidence for the simultaneous occurrence of a third DNP mechanism, in addition to the Overhauser and solid effects (Gizatullin et al., 2022).

In the companion paper (Sezer, 2023a) we showed that when the microwave (mw) nutation frequency approaches the nuclear Larmor frequency, as could be the case at X band under the high mw powers used in the DNP experiments (Neudert et al., 2016), the forbidden transitions of the solid effect are no longer shifted by ±ω_I with respect to the electronic resonance but come closer together and may even coalesce. In such cases, the theoretical justification for modeling the solid-effect field profile by shifting the EPR line shape by ±ω_I becomes questionable. Unfortunately, the analytical expressions of Sezer (2023a), which remain valid in this regime, are not applicable to liquids since they do not account for the modulation of the dipolar interaction by molecular diffusion. In the current paper, the time-domain description of the solid effect from Sezer (2023a) is extended to liquids.

Recently, Kuzhelev et al. (2022) reported proton DNP enhancements in the liquid phase of lipid bilayers at 9.4 T (260 GHz) using the free radical BDPA as a polarizing agent. The large nuclear Larmor frequency (400 MHz) enabled a clear spectral separation of the acyl-chain protons and the polar protons of water. For the non-polar protons, maximal enhancements were observed at the canonical resonance positions of the zero- and double-quantum forbidden transitions, characteristic of the solid effect, while enhancement due to the Overhauser effect was missing. The large spectral separation of the positive and negative enhancements and the narrow line of BDPA provided a uniquely “clean” access to the solid effect in a viscous liquid environment (Kuzhelev et al., 2022).

In addition to the maximum enhancements at the canonical solid-effect offsets, the DNP spectrum of Kuzhelev et al. (2022) revealed additional enhancement peaks at much smaller offsets. These were postulated to arise due to the DNP mechanism known as thermal mixing. However, thermal mixing is commonly associated with a broad EPR spectrum (Wenckebach, 2021), while the spectrum of BDPA was extremely narrow at the elevated DNP concentration used in the experiment, and the observed spurious enhancement peaks lay outside this narrow spectrum (Kuzhelev et al., 2022, Fig. 2). Here we explain the entire DNP spectrum, including the puzzling features at low offsets, considering only one electronic and one nuclear spin.

The rest of the paper is organized as follows. To account for molecular diffusion, in Sect. 3 we transform the equations of motion of Sezer (2023a) into stochastic Liouville equations (SLEs) (Kubo, 1969). After taking into account that all relevant timescales are orders of magnitude shorter than the nuclear spin-lattice relaxation time, the SLE formalism yields the time-correlation function of the dipolar interaction. In Sect. 4 we show that the solid-effect lines in the DNP spectrum (i.e., those shifted by ±ω_I) experience additional motional broadening compared to the homogeneous EPR line width. As a result, the tails of these lines around the position of the electronic resonance increase substantially. Under favorable conditions, the product of these tails with the dispersive EPR component may become sufficiently large to be visible as separate enhancement peaks in the DNP field profile. We attribute the non-canonical peaks in the DNP spectrum of Kuzhelev et al. (2022) to this phenomenon. Our conclusions are presented in Sect. 5. The next section summarizes the needed background.

2 Motivation and background

2.1 Dynamic nuclear polarization in liquids

The transfer of polarization in solids involves two mechanistically different steps (Hovav et al., 2010; Smith et al., 2012). The first one is the direct polarization of the nuclear spins that are sufficiently close to the free radical and have appreciable dipole–dipole interaction with the electronic spin. Being closest to the unpaired electron, the nuclei on the free radical itself benefit most from this first step of direct polarization (Tan et al., 2019; Delage-Laurin et al., 2021). Polarizing the intramolecular nuclei in this way, however, is not particularly useful unless the polarization can spread to the rest of the sample. This is where the second step comes in. In this step, the polarization spreads from the directly polarized nuclei to the distant nuclei by spin diffusion. Because it relies on the relatively weak dipole–dipole interaction between the nuclear spins, spin diffusion is slow and is often the bottleneck for efficient polarization transfer in the solid state (Hovav et al., 2010; Wiśniewski et al., 2016).

In liquids, spin diffusion is not efficient because the nuclei constantly change their positions due to random thermal motions. However, since molecular diffusion moves the nuclei across nanometer distances in nanoseconds and thus rapidly spreads the polarization of the directly polarized nuclei across the sample, spin diffusion is also not needed. Taking glycerol as an example, with a self-diffusion coefficient of $\mathrm{6.6}\times {\mathrm{10}}^{-\mathrm{3}}$ nm² ns⁻¹ at 40 ^∘C (Tomlinson, 1973), which is 500 times less than the self-diffusion coefficient of water at the same temperature (Holz et al., 2000), it is a rather viscous liquid. Nevertheless, at the relatively small radical concentration of 1 mM, a molecule of glycerol covers the average distance between two radicals in less than 400 ns. This is at least 5 orders of magnitude less than the nuclear T₁ of protons, even after accounting for paramagnetic relaxation. Any given solvent nucleus will thus encounter the electronic spins a million times during its T₁ relaxation time. Even in viscous liquids, therefore, molecular diffusion is expected to homogenize the nuclear polarization across the sample during times that are orders of magnitude shorter than the nuclear T₁. This advantage of liquids over solids, however, comes at a price: the polarization of the nuclei on the free radical is no longer accessible to the solvent, and proximal solvent nuclei have to be polarized directly in the first step of polarization transfer (Fig. 1).

Figure 1 Two conceptually different steps of the polarization transfer process in liquids. (1) Direct transfer from the electronic spin on the free radical to the proximate nuclear spins on the solvent molecules due to dipolar interaction. (2) Diffusion of the proximate solvent molecules to the bulk.

Given that every solvent molecule gets directly polarized and also spreads the polarization, the distinction between the two steps of polarization transfer in liquids (Fig. 1) is conceptual and does not reflect fundamental differences in the mechanisms of the two steps. In fact, both steps are enabled by molecular diffusion which sometimes brings a solvent molecule closer to the radical and sometimes takes it further away. Since the analytical description of translational diffusion in simple liquids is well developed (Ayant et al., 1975; Hwang and Freed, 1975), a unified theoretical treatment of the two steps of polarization transfer becomes possible, as we demonstrate in the present paper.

From the six terms of the dipolar alphabet, the part that contributes to the solid effect is $A_{\mathrm{1}}{S}_z{I}_{+}+A_{\mathrm{1}}^{*}{S}_z{I}_{-}$ (Wenckebach, 2016), where

$$\begin{array}{}\text{(1)}& A_{\mathrm{1}}=D_{\mathrm{dip}}{\displaystyle \frac{-\mathrm{3}\cos \mathit{\theta }\sin \mathit{\theta }{e}^{-i\mathit{\varphi }}}{r^{\mathrm{3}}}}.\end{array}$$

Here $D_{\mathrm{dip}}=\left({\mathit{\mu }}_{\mathrm{0}}/\mathrm{4}\mathit{\pi }\right)\mathrm{\hslash }{\mathit{\gamma }}_S{\mathit{\gamma }}_I$ is the dipolar constant, γ_S and γ_I are the gyromagnetic ratios of the electronic and nuclear spins, and $(r,\mathit{\theta },\mathit{\varphi })$ are the spherical polar coordinates of the vector pointing from one of the spins to the other. The angular dependence of A₁ is that of a second-degree spherical harmonic of order m=1, as implied by the subscript. The need for direct polarization of the solvent nuclei in liquids increases the shortest possible distance r in Eq. (1) and thus reduces the largest achievable dipolar coupling. This requirement for interaction across a larger distance, however, does not explain why the solid effect works in solids but is compromised in liquids.

To understand the difference between solids and liquids one should consider the time-correlation function of the dipolar interaction:

$$\begin{array}{}\text{(2)}& C_{\mathrm{11}}\left(t\right)=\langle \langle A_{\mathrm{1}}^{*}(t^{\prime }+t)A_{\mathrm{1}}\left(t^{\prime }\right){\rangle }_{t^{\prime }}\rangle .\end{array}$$

Here the inner angular brackets with the subscript t^′ denote averaging with respect to the time point t^′ along the random trajectory of a single nuclear spin. Because every nucleus encounters the electronic spins millions of times during its T₁ relaxation time, this average should be the same for all nuclei in the liquid. Thus, in addition to the time averaging, in Eq. (2) we also average over the ensemble of identical nuclear spins in the sample (outer angular brackets).

Now, if the dipolar correlation function (Eq. 2) decays on timescales that are much longer than some relevant characteristic time, then the experiment essentially detects the initial value $C_{\mathrm{11}}\left(\mathrm{0}\right)=\langle \langle A_{\mathrm{1}}^{*}\left(t^{\prime }\right)A_{\mathrm{1}}\left(t^{\prime }\right){\rangle }_{t^{\prime }}\rangle =\langle A_{\mathrm{1}}^{*}{A}_{\mathrm{1}}\rangle$. The last ensemble average over all electron–nucleus pairs requires integration over the spatial variables $(r,\mathit{\theta },\mathit{\varphi })$ and multiplication by the concentration N of the unpaired electrons:

$$\begin{array}{}\text{(3)}& \langle A_{\mathrm{1}}^{*}{A}_{\mathrm{1}}\rangle =N{D}_{\mathrm{dip}}^{\mathrm{2}}{\displaystyle \frac{\mathrm{24}\mathit{\pi }}{\mathrm{5}}}\underset{b}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\mathrm{d}r}{r^{\mathrm{4}}}}=D_{\mathrm{dip}}^{\mathrm{2}}{\displaystyle \frac{\mathrm{24}\mathit{\pi }}{\mathrm{5}}}{\displaystyle \frac{N}{\mathrm{3}{b}^{\mathrm{3}}}}.\end{array}$$

(The factor $\mathrm{24}\mathit{\pi }/\mathrm{5}$ comes from the normalization of the spherical harmonic $Y_{\mathrm{2}}^{\mathrm{1}}$.) This slow-motional limit corresponds to the situation in solids under the (unrealistic) assumption of fast and efficient spin diffusion. If, on the other hand, the decay time of the correlation function is much shorter than the relevant characteristic timescale, then the experiment detects the long-time limit $C_{\mathrm{11}}\left(\mathrm{\infty }\right)=\langle \langle A_{\mathrm{1}}^{*}\left(\mathrm{\infty }\right)A_{\mathrm{1}}\left(t^{\prime }\right){\rangle }_{t^{\prime }}\rangle =\langle A_{\mathrm{1}}^{*}\rangle \langle A_{\mathrm{1}}\rangle$. The solid effect vanishes because the average of the spherical harmonic $Y_{\mathrm{2}}^{\mathrm{1}}$ over the angles gives 〈A₁〉=0. This fast-motional limit corresponds to low-viscosity liquids in which the dipolar interaction is averaged out. To the extent that they exhibit the solid effect, viscous liquids must lie somewhere between these two extremes.

The interpolation between these two limiting cases on the basis of the dipolar correlation function is formally developed in Sect. 3. This task requires a time-domain description of the solid effect, similar to the treatment of relaxation by random motion where the correlation function arises from second-order, time-dependent perturbation theory (Abragam, 1961, chap. VIII). In principle, there are two such time-domain descriptions that we can utilize for the treatment of the solid effect in liquids. The first is the rate-equation formalism, which models the dynamics of the electronic and nuclear polarizations, and the second is the description developed in Sezer (2023a), which additionally accounts for the dynamics of the coherences. Both of these options will be explored in Sect. 3. When modeling the stochastic dynamics of the dipolar interaction, we resort to the stochastic Liouville equation (SLE) of Kubo (1954) and Anderson (1954), rather than to second-order perturbation theory. In agreement with previous work (Papon et al., 1968; Leblond et al., 1971a), our analysis shows that the characteristic timescale against which the dipolar correlation time should be compared is the electronic T₂ relaxation time.

In the next two subsections we summarize the time-domain analysis of Sezer (2023a).

2.2 Rate equations

The dynamics of the electronic polarization, P_S, is justifiably taken to be independent of the dipolar interaction with the nuclear spins, as other mechanisms relax the electrons more efficiently. With R_1S denoting the rate of electronic T₁ relaxation, and v₁ the rate constant of the mw excitation of the (allowed) EPR transition, the rate equation of the electronic polarization is

$$\begin{array}{}\text{(4)}& {\dot{P}}_S=-R_{\mathrm{1}S}(P_S-P_S^{\mathrm{eq}})-\mathrm{2}{v}_{\mathrm{1}}{P}_S.\end{array}$$

Here $P_S^{\mathrm{eq}}$ is the equilibrium (Boltzmann) electronic polarization and the dot over the symbol denotes differentiation with respect to time. Solving this equation for $P_S^{\mathrm{ss}}$ at steady state, we arrive at the ratio

$$\begin{array}{}\text{(5)}& p={\displaystyle \frac{P_S^{\mathrm{ss}}}{P_S^{\mathrm{eq}}}}={\displaystyle \frac{R_{\mathrm{1}S}}{R_{\mathrm{1}S}+\mathrm{2}{v}_{\mathrm{1}}}}=\mathrm{1}-s,\end{array}$$

where s is the familiar electronic saturation factor. We refer to p as the electronic polarization factor, since p=0 indicates that the steady-state polarization has vanished, and p=1 indicates that it is identical to the Boltzmann polarization (i.e., maximally polarized).

The rate equation of the nuclear polarization, P_I, is

$$\begin{array}{}\text{(6)}& {\dot{P}}_I=-R_{\mathrm{1}I}\left(P_I-P_I^{\mathrm{eq}}\right)-v_{+}{P}_I-v_{-}{P}_S,\end{array}$$

where R_1I is the nuclear T₁ relaxation rate, and the phenomenological rate constants v_± quantify the mw excitation of the forbidden transitions. The steady state of Eq. (6) is

$$\begin{array}{}\text{(7)}& R_{\mathrm{1}I}\left(P_I^{\mathrm{ss}}-P_I^{\mathrm{eq}}\right)=-v_{+}{P}_I^{\mathrm{ss}}-p{v}_{-}{P}_S^{\mathrm{eq}},\end{array}$$

where $P_S^{\mathrm{ss}}=p{P}_S^{\mathrm{eq}}$ was used in the last term.

As the derivations in Sect. 3 consider only the effect of mw excitation, we have written Eq. (7) such that the relaxation contribution is on the left and the mw contribution on the right of the equality. Subsequently, to identify the phenomenological rate constants v_±, we will match the terms on the right-hand side with the predictions of the proper analysis in liquids.

From Eq. (7) we find the DNP enhancement

$$\begin{array}{}\text{(8)}& \mathit{ϵ}=P_I^{\mathrm{ss}}/P_I^{\mathrm{eq}}-\mathrm{1}={\mathit{ϵ}}_{\mathrm{SE}}-\left(\mathrm{1}-p_X\right),\end{array}$$

where the first equality is the definition of ϵ and

$$\begin{array}{}\text{(9)}& {\mathit{ϵ}}_{\mathrm{SE}}={\displaystyle \frac{p{v}_{-}}{R_{\mathrm{1}I}+v_{+}}}{\displaystyle \frac{\left|{\mathit{\gamma }}_S\right|}{{\mathit{\gamma }}_I}},p_X={\displaystyle \frac{R_{\mathrm{1}I}}{R_{\mathrm{1}I}+v_{+}}}.\end{array}$$

From ϵ_SE it is clear that the solid effect benefits from large pv₋ and small $R_{\mathrm{1}I}+v_{+}$. The ratio p_X in Eq. (9) is analogous to the electronic polarization factor Eq. (5), and we call it the nuclear cross-polarization factor. In liquids, v₊ is typically negligible compared to the nuclear spin-lattice relaxation rate, and p_X≈1. Then,

$$\begin{array}{}\text{(10)}& {\mathit{ϵ}}_{\mathrm{SE}}\approx \left(p{v}_{-}\right)T_{\mathrm{1}I}\left|{\mathit{\gamma }}_S\right|/{\mathit{\gamma }}_I(p_X\approx \mathrm{1}).\end{array}$$

2.3 Spin dynamics

The dynamics of the quantum-mechanical expectation values s_n of the electronic spin operators S_n ($n=x,y,z$) is described by the classical Bloch equations

$$\begin{array}{}\text{(11)}& \left[\begin{array}{c}{\dot{s}}_x\\ {\dot{s}}_y\\ {\dot{s}}_z\end{array}\right]=-\left[\begin{array}{ccc}{R}_{\mathrm{2}S}& \mathrm{\Omega }& \mathrm{0}\\ -\mathrm{\Omega }& R_{\mathrm{2}S}& {\mathit{\omega }}_{\mathrm{1}}\\ \mathrm{0}& -{\mathit{\omega }}_{\mathrm{1}}& R_{\mathrm{1}S}\end{array}\right]\left[\begin{array}{c}{s}_x\\ s_y\\ s_z\end{array}\right]+R_{\mathrm{1}S}\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ s_z^{\mathrm{eq}}\end{array}\right].\end{array}$$

The matrix in Eq. (11) contains the electronic transverse relaxation rate R_2S, the mw nutation frequency ω₁ and the offset $\mathrm{\Omega }={\mathit{\omega }}_S-\mathit{\omega }$ between the electronic Larmor frequency ω_S and the mw frequency ω.

In Sezer (2023a) we visualized such coupled differential equations diagrammatically. In our visual depiction, the time derivative of a dynamical variable, like s_n, is represented by an oval node. The contributions to this time derivative, which are on the right-hand side of the differential equation, are represented by arrows that flow into that node (Fig. 2a). The contribution of a given arrow is obtained by multiplying the weight of the arrow by the variable from which the arrow originates. The self-arrows that exit from an oval node and enter the same node correspond to the relaxation terms along the diagonal of the Bloch matrix. The negative sign of the weight of a self-arrow is written separately inside the loop formed by the arrow. The constant variable $s_z^{\mathrm{eq}}$ in the inhomogeneous term of the Bloch equations (Eq. 11) is represented by a gray rectangular node. With this notation, the Bloch equations are depicted by the four nodes in the top row of Fig. 2a and by the black arrows connecting these nodes.

Figure 2 Diagrammatic representation of (a) the equations of motion of the spin dynamics relevant to the solid effect, namely Eq. (11) (purple rectangle), Eq. (12) (cyan rectangle) and Eq. (14), and (b) the rate equations of the electronic and nuclear polarizations (Eqs. 4 and 6). (c) Steady-state relationships between the input and output variables.

The lower half of Fig. 2a shows the dynamics of the electron–nuclear coherences that are relevant for the solid effect. In particular, the quantum-mechanical expectation values of the operators S_nI₊ ($n=x,y,z$), which we denote interchangeably by g_n and s_ni₊, evolve according to the following coupled differential equations (Sezer, 2023a):

$$\begin{array}{}\text{(12)}& \left[\begin{array}{c}{\dot{g}}_x\\ {\dot{g}}_y\\ {\dot{g}}_z\end{array}\right]=-\mathfrak{B}\left[\begin{array}{c}{g}_x\\ g_y\\ g_z\end{array}\right]-i{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{A}_{\mathrm{1}}\left[\begin{array}{c}-i{s}_y\\ i{s}_x\\ i_z\end{array}\right],\end{array}$$

where

$$\begin{array}{}\text{(13)}& \mathfrak{B}=\left[\begin{array}{ccc}{R}_{\mathrm{2}S}+i{\mathit{\omega }}_I& \mathrm{\Omega }& \mathrm{0}\\ -\mathrm{\Omega }& R_{\mathrm{2}S}+i{\mathit{\omega }}_I& {\mathit{\omega }}_{\mathrm{1}}\\ \mathrm{0}& -{\mathit{\omega }}_{\mathrm{1}}& R_{\mathrm{1}S}+i{\mathit{\omega }}_I\end{array}\right].\end{array}$$

The matrix 𝔅 is essentially the Bloch matrix but with the nuclear Larmor frequency ω_I added as an imaginary part to its main diagonal. The time derivatives of g_n in Eq. (12) are represented by the three oval nodes enclosed in the cyan rectangle in Fig. 2a. The arrows between these nodes are seen to exactly replicate the classical Bloch equations in the rectangle above them. The inhomogeneous term in Eq. (12) couples the dynamics of the variables g_n to the transverse components of the electronic magnetization, on the one hand, and to the longitudinal component of the nuclear magnetization, on the other. All these couplings scale with the dipolar interaction A₁. They play an essential role in the solid effect, as they connect the Boltzmann electronic polarization to the nuclear polarization (labeled “Input” and “Output” in Fig. 2a).

Lastly, the coherent dynamics of the operator I_z, whose expectation value is denoted by i_z, is

$$\begin{array}{}\text{(14)}& {\dot{i}}_z{|}_{\mathrm{coh}}=-\mathrm{Re}\mathit{\left\{}i{A}_{\mathrm{1}}^{*}{g}_z\mathit{\right\}},\end{array}$$

where Re takes the real part of a complex number. In Fig. 2a this equation is represented by the oval node i_z (labeled “Output”) and the red arrow flowing into it.

In liquids, the weights A₁ fluctuate randomly due to molecular diffusion. When extending the formalism to liquids (Sect. 3.3), we will transform Eq. (12), which constitutes a system of coupled differential equations, to an SLE (Kubo, 1969) that describes the spin dynamics under random modulation of A₁.

For comparison, Fig. 2b shows the dynamics of the longitudinal components implied by the rate equations (Eqs. 4 and 6). Visual inspection of Fig. 2a and b makes it clear that the rate constant v₋ provides a reduced description of the complicated network connecting s_z to i_z. Similarly, the rate constant v₊ accounts for the self-influence of i_z mediated by the coherences in the second set of Bloch equations (enclosed in the cyan rectangle). In (Sezer, 2023a), we identified the rates v_± and v₁ by requiring that the dynamics in Fig. 2a and b reached identical steady states.

At steady state, the three dynamical variables of the classical Bloch equations (Eq. 11) were related to each other and to the electronic Boltzmann polarization as follows:

$$\begin{array}{}\text{(15)}& s_x^{\mathrm{ss}}=\left({\mathit{\omega }}_{\mathrm{1}}{f}_x\right)s_z^{\mathrm{ss}},s_y^{\mathrm{ss}}=-\left({\mathit{\omega }}_{\mathrm{1}}{f}_y\right)s_z^{\mathrm{ss}},s_z^{\mathrm{ss}}=\left(R_{\mathrm{1}S}{f}_z\right)s_z^{\mathrm{eq}},\end{array}$$

where

$$\begin{array}{}\text{(16)}& f_y={\displaystyle \frac{R_{\mathrm{2}S}}{R_{\mathrm{2}S}^{\mathrm{2}}+{\mathrm{\Omega }}^{\mathrm{2}}}},f_x={\displaystyle \frac{\mathrm{\Omega }}{R_{\mathrm{2}S}}}{f}_y,f_z={\displaystyle \frac{\mathrm{1}}{R_{\mathrm{1}S}+{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}{f}_y}}.\end{array}$$

(Comparing the last equality in Eq. (15) with Eq. (5) we found that p=R_1Sf_z and that the rate constant of the allowed EPR transition was $v_{\mathrm{1}}={\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}{f}_y/\mathrm{2}$.)

Solving the second set of Bloch equations (Eq. 12) at steady state, and substituting $g_z^{\mathrm{ss}}$ in Eq. (14), we obtain

$$\begin{array}{}\text{(17)}& {\dot{i}}_z{|}_{\mathrm{coh}}^{\mathrm{ss}}=-{\mathit{\delta }}^{\mathrm{2}}\mathrm{Re}\left\{\left[\begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]{\mathfrak{B}}^{-\mathrm{1}}\left[\begin{array}{c}-i{s}_y^{\mathrm{ss}}\\ i{s}_x^{\mathrm{ss}}\\ i_z^{\mathrm{ss}}\end{array}\right]\right\},\end{array}$$

where the dipolar interaction is isolated in

$$\begin{array}{}\text{(18)}& {\mathit{\delta }}^{\mathrm{2}}=\left(A_{\mathrm{1}}^{*}{A}_{\mathrm{1}}\right)/\mathrm{4}.\end{array}$$

Since the transverse components $s_{x,y}^{\mathrm{ss}}$ are related to the Boltzmann polarization (Eq. 15), the right-hand side of Eq. (17) is of the form

$$\begin{array}{}\text{(19)}& {\dot{i}}_z{|}_{\mathrm{coh}}^{\mathrm{ss}}=-{\mathit{\delta }}^{\mathrm{2}}\left(T_i{i}_z^{\mathrm{ss}}+T_s{s}_z^{\mathrm{eq}}\right),\end{array}$$

where

$$\begin{array}{}\text{(20)}& \begin{array}{rl}{T}_i& =\left[\begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]\mathrm{Re}\mathit{\left\{}{\mathfrak{B}}^{-\mathrm{1}}\mathit{\right\}}{\left[\begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]}^{\mathrm{T}}=\mathrm{Re}\mathit{\left\{}{\mathfrak{B}}_{\mathrm{33}}^{-\mathrm{1}}\mathit{\right\}}\\ T_s& =\left[\begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]\mathrm{Re}\mathit{\left\{}i{\mathfrak{B}}^{-\mathrm{1}}\mathit{\right\}}{\left[\begin{array}{ccc}p{\mathit{\omega }}_{\mathrm{1}}{f}_y& p{\mathit{\omega }}_{\mathrm{1}}{f}_x& \mathrm{0}\end{array}\right]}^{\mathrm{T}}.\end{array}\end{array}$$

(The superscript “T” indicates transpose, and ${\mathfrak{B}}_{ij}^{-\mathrm{1}}$ is the ijth element of the inverse matrix 𝔅⁻¹.) Because 𝔅 has units of inverse time, T_i and T_s have units of time.

We note that T_s receives contributions from both $s_y^{\mathrm{ss}}$ and $s_x^{\mathrm{ss}}$. As was shown in Sezer (2023a), it is possible to rewrite the contribution of the former as if it also came from $s_x^{\mathrm{ss}}$. In other words, the entire contribution of $s_z^{\mathrm{eq}}$ to the derivative of i_z at steady state can be expressed as if it is mediated only through the dispersive component $s_x^{\mathrm{ss}}$, as depicted in Fig. 2c. (Dashed arrows represent mathematical relationships between the variables that hold at steady state. Differently from the solid arrows, which correspond to causal dependencies governing the dynamics at all times, the dashed arrows need not reflect direct causal dependence. A rectangular node indicates that the inflowing arrows contribute directly to the value of the variable and not to its time derivative. The gray shade of the nodes signals that the variables remain constant in time, as they should at steady state.)

In addition to ${\mathfrak{B}}_{\mathrm{33}}^{-\mathrm{1}}$, in Eq. (20) we also need ${\mathfrak{B}}_{\mathrm{31}}^{-\mathrm{1}}$ and ${\mathfrak{B}}_{\mathrm{32}}^{-\mathrm{1}}$. These are ${\mathfrak{B}}_{\mathrm{31}}^{-\mathrm{1}}=\left({\mathit{\omega }}_{\mathrm{1}}{F}_x\right)F_z$, ${\mathfrak{B}}_{\mathrm{32}}^{-\mathrm{1}}=\left({\mathit{\omega }}_{\mathrm{1}}{F}_y\right)F_z$ and ${\mathfrak{B}}_{\mathrm{33}}^{-\mathrm{1}}=F_z$, where the functions

$$\begin{array}{}\text{(21)}& \begin{array}{rl}{F}_y& ={\displaystyle \frac{R_{\mathrm{2}S}+i{\mathit{\omega }}_I}{{\left(R_{\mathrm{2}S}+i{\mathit{\omega }}_I\right)}^{\mathrm{2}}+{\mathrm{\Omega }}^{\mathrm{2}}}},F_x={\displaystyle \frac{\mathrm{\Omega }}{R_{\mathrm{2}S}+i{\mathit{\omega }}_I}}{F}_y\\ F_z& ={\displaystyle \frac{\mathrm{1}}{R_{\mathrm{1}S}+i{\mathit{\omega }}_I+{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}{F}_y}}\end{array}\end{array}$$

play an analogous role in the steady-state analysis of the second set of Bloch equations as the functions f_y, f_x and f_z (Eq. 16) in the classical Bloch equations. In terms of these,

$$\begin{array}{}\text{(22)}& T_i=\mathrm{Re}\mathit{\left\{}{F}_z\mathit{\right\}},T_s=\left(p{\mathit{\omega }}_{\mathrm{1}}{f}_x\right)T_x,\end{array}$$

where

$$\begin{array}{}\text{(23)}& T_x=\mathrm{Re}\mathit{\left\{}i{F}_z\right({\mathit{\omega }}_{\mathrm{1}}{F}_y^{\prime }\left)\mathit{\right\}},F_y^{\prime }={\displaystyle \frac{\mathrm{2}{R}_{\mathrm{2}S}+i{\mathit{\omega }}_I}{R_{\mathrm{2}S}+i{\mathit{\omega }}_I}}{F}_y.\end{array}$$

(The functions in Eq. 23 emerge from lumping the contribution of $s_y^{\mathrm{ss}}$ to that of $s_x^{\mathrm{ss}}$, as mentioned above.)

With T_i and T_s determined, the forbidden-transition rates on the right-hand side of Eq. (7) become

$$\begin{array}{}\text{(24)}& v_{+}={\mathit{\delta }}^{\mathrm{2}}\left(T_i-T_i^{\mathrm{0}}\right),p{v}_{-}={\mathit{\delta }}^{\mathrm{2}}{T}_s=\left(p{\mathit{\omega }}_{\mathrm{1}}{f}_x\right){\mathit{\delta }}^{\mathrm{2}}{T}_x,\end{array}$$

where the mw-independent part of T_i, namely

$$\begin{array}{}\text{(25)}& T_i^{\mathrm{0}}=\mathrm{Re}\left\{F_z\left({\mathit{\omega }}_{\mathrm{1}}=\mathrm{0}\right)\right\}={\left(R_{\mathrm{1}S}+i{\mathit{\omega }}_I\right)}^{-\mathrm{1}},\end{array}$$

is subtracted in the first equality of Eq. (24) because it contributes to the nuclear T₁ relaxation rate. (In Fig. 2a this mw-independent part corresponds to the loop formed by the green arrow from i_z to g_z and the red arrow in the opposite direction.)

2.4 The solid effect

Using the rate constants in Eq. (24), we rewrite the solid-effect DNP enhancement (Eq. 9) as

$$\begin{array}{}\text{(26)}& {\mathit{ϵ}}_{\mathrm{SE}}={\displaystyle \frac{\left(p{\mathit{\omega }}_{\mathrm{1}}{f}_x\right)T_x}{R_{\mathrm{1}I}/{\mathit{\delta }}^{\mathrm{2}}+\left(T_i-T_i^{\mathrm{0}}\right)}}{\displaystyle \frac{\left|{\mathit{\gamma }}_S\right|}{{\mathit{\gamma }}_I}}.\end{array}$$

The functions pω₁f_x, T_x and $(T_i-T_i^{\mathrm{0}})$ in this expression are visualized in, respectively, the first, second and third rows of Fig. 3. The product of the first two rows, which appears in the numerator of Eq. (26), is shown in the fourth row of the figure. In the right margin of the figure, we have included the flow diagram from Fig. 2c, which has been straightened here so that the weights of the arrows correspond to the respective rows. Note that the dipolar interaction strength, δ, and the nuclear spin-lattice relaxation rate were not needed to calculate the properties in the first four rows of the figure. (They will be needed for the last two rows.)

Figure 3 Comparison between solids (with fast spin diffusion) and liquids. First row: dispersive component of the power-broadened EPR line. Second and third rows: relevant transfer functions of the new Bloch equations in Fig. 2a. Fourth row: the product of the first and second rows, which relates the input $s_z^{\mathrm{eq}}$ to the output $i_z^{\mathrm{ss}}$. Fifth row: DNP enhancement calculated from the third and fourth rows using $R_{\mathrm{1}I}/\langle {\mathit{\delta }}^{\mathrm{2}}\rangle$. Sixth row: nuclear cross-polarization factor calculated from the third row using $R_{\mathrm{1}I}/\langle {\mathit{\delta }}^{\mathrm{2}}\rangle$. Simulation parameters: T_2S=60 ns, T_1S=9 T_2S, B₁=6 G (converted to ω₁ assuming g=2), contact distance b=1 nm, radical concentration N=0.1 M, and T_1I of 4.7 ms (X band), 27.4 ms (Q band) and 50 ms (W and J bands). The dipolar correlation time of the liquid simulation (dashed black lines) is $\mathit{\tau }=T_{\mathrm{2}S}/\mathrm{5}=\mathrm{12}$ ns.

While different magnetic fields B₀ yield different relaxation times, for illustrative purposes we used the same electronic T₁ and T₂ times for X, Q, W and J bands. We additionally used the same mw field (B₁=6 G) at all bands. Hence, the steady state of the classical Bloch equations (Fig. 3, first row) is identical across the four columns of the figure. The solid blue lines, which correspond to the dispersive component of the power-broadened EPR line, are identical but appear different due to the different scalings of the horizontal axes. The absorptive component is much smaller under the power-broadening conditions considered here and is not shown. However, its contribution is exactly accounted for in the analysis. (This was the reason for introducing the functions in Eq. 23.) Anticipating the liquid state, we observe that the classical Bloch equations are independent of the dipolar coupling A₁ (Fig. 2a). Hence, the first row of Fig. 3 will not change when going to liquids because we use identical relaxation rates.

The transfer functions T_x and $(T_i-T_i^{\mathrm{0}})$ encapsulate all relevant steady-state properties of the second set of Bloch equations, as well as their coupling to the classical Bloch equations and to i_z through the dipolar interaction (Fig. 2a). These functions are visualized in the second and third rows of Fig. 3 (orange and red lines). The solid colored lines are calculated using the equations given above and correspond to solids, subject to the (unrealistic) assumption of very fast spin diffusion. The dashed black lines are calculated as described in the next section and correspond to liquids. Clearly, the time-dependent modulation of the dipolar interaction in liquids has a dramatic effect on these functions.

The fourth row in Fig. 3 shows the product of the blue lines in the first row and the orange lines in the second row and corresponds to the total transfer function from the primary input, $s_z^{\mathrm{eq}}$, to the ultimate output, i_z (Fig. 2a). With the exception of X band, going from solids to liquids substantially reduces the peaks of pv₋. (We used identical relaxation parameters for liquids and solids to highlight the role of the dipolar correlation time.)

The transfer functions in the first four rows of Fig. 3 depend only on the electronic relaxation times (assuming $T_{\mathrm{2}I}\gg T_{\mathrm{1}S},T_{\mathrm{2}S}$). To calculate the enhancement ϵ_SE and the nuclear cross-polarization factor, p_X, which are shown in the last two rows of Fig. 3, we had to select specific values for R_1I and δ². For the latter, we used the ensemble-averaged static value from Eq. (3),

$$\begin{array}{}\text{(27)}& \langle {\mathit{\delta }}^{\mathrm{2}}\rangle ={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}〈A_{\mathrm{1}}^{*}{A}_{\mathrm{1}}〉=D_{\mathrm{dip}}^{\mathrm{2}}{\displaystyle \frac{\mathrm{6}\mathit{\pi }}{\mathrm{5}}}{\displaystyle \frac{N}{\mathrm{3}{b}^{\mathrm{3}}}},\end{array}$$

which applies to solids with fast and efficient spin diffusion. The numerical calculations in Fig. 3 are for contact distance b_ref=1 nm and radical concentration N_ref=0.1 M. These are realistic but otherwise arbitrary values.

For the purposes of illustration we wanted to use the same R_1I for all four mw bands in the figure. In this way, by comparing the four columns with each other, one would be able to assess the effect of changing only B₀. This strategy worked for solids, at least for the numerical values that were used, but failed for liquids due to the very different contributions of T_i to the nuclear relaxation rate (denoted by $T_i^{\mathrm{0}}$ in Eq. 25). This part of T_i is shown in the third row of Fig. 3 with horizontal dotted lines. The dotted red lines for solids are very close to zero. The more visible dotted black lines for liquids change dramatically with the mw band. Since the total relaxation rate R_1I must be larger than the contribution of $T_i^{\mathrm{0}}$, the nuclear T₁ had to be only a few ms at X band. Using such small T₁ at J band, however, gave tiny liquid-state DNP enhancements.

Even if, admittedly, our calculated enhancements are only illustrative, in an effort to have somewhat realistic nuclear T₁ relaxation times, we used T_1I=50 ms when the corresponding rate R_1I was larger than $\langle {\mathit{\delta }}^{\mathrm{2}}\rangle T_i^{\mathrm{0}}$ and used $R_{\mathrm{1}I}/\langle {\mathit{\delta }}^{\mathrm{2}}\rangle =\mathrm{2}{T}_i^{\mathrm{0}}$ otherwise. This resulted in the following nuclear T₁ relaxation times: 4.7 ms (X band), 27.4 ms (Q band) and 50 ms (W and J bands). These were used for both liquids and solids. Naturally, the choice of different nuclear T₁ times has a direct influence on the calculated enhancements. For example, the peak DNP enhancements at X and W bands differ by about 1 order of magnitude (purple lines in the fifth row of Fig. 3) mostly because the nuclear T₁ times at these two bands also differ by 1 order of magnitude.

The theory behind the liquid calculations in Fig. 3 is presented in the next section.

3 Liquids

3.1 Molecular motion as a random process

Let us denote the components $(r,\mathit{\theta },\mathit{\varphi })$ of the inter-spin vector collectively by ζ. To describe the solid effect in liquids we consider a random process that changes ζ and thus modulates the dipolar interaction between the two types of spins.

When the random dynamics of ζ is modeled as a discrete-state process, the probabilities of observing the different discrete states are collected in the vector p(t). This probability vector evolves in time as $\dot{\mathit{p}}\left(t\right)=-\mathbf{K}\mathit{p}\left(t\right)$, where the matrix K contains the rate constants of the random transitions between the states. All eigenvalues of such stochastic matrices are non-negative, and, for an ergodic chain of states, only one of the eigenvalues equals zero. In general, the stochastic matrix K is not symmetric, which means that there are a right eigenvector and a left eigenvector associated with each eigenvalue. The right eigenvector of the zero eigenvalue corresponds to the vector of equilibrium probabilities, p^eq, and the left eigenvector of the zero eigenvalue corresponds to the vector 1, which contains ones in all of its entries. Note that 1^Tp^eq=1.

When the random dynamics of ζ is modeled as a continuous-state diffusion process, then the time evolution of the probability density p(ζ,t) is described by a Fokker–Planck equation of the form

$$\begin{array}{}\text{(28)}& {\displaystyle \frac{\partial p(\mathit{\zeta },t)}{\partial t}}=-K_{\mathit{\zeta }}p(\mathit{\zeta },t),\end{array}$$

where K_ζ is a linear differential operator acting on the ζ dependence of p(ζ,t). As in the discrete case, the eigenvalues of K_ζ would be non-negative, and one eigenvalue would equal zero. The corresponding right eigenfunction is the equilibrium probability density p^eq(ζ), and the left eigenfunction is constant in ζ.

For brevity, we will also adopt the discrete notation for the continuous case. In particular, we will use italic bold symbols to indicate the dependence on ζ and will denote operators that act on the ζ dependence with non-italicized capital bold symbols. With this understanding,

$$\begin{array}{}\text{(29)}& \dot{\mathit{p}}\left(t\right)=-\mathbf{K}\mathit{p}\left(t\right)\end{array}$$

will apply to both the continuous and discrete cases. Similarly, 1^Tf will imply integration over the ζ dependence of the function f(ζ) in the continuous case and summation over all different states in the discrete case.

Note that in the probabilistic description of the random process by the Fokker–Planck equation (Eq. 28), the probability density p(ζ,t) characterizes an ensemble of nuclei, and ζ is treated as an independent variable which is not a function of t. In contrast, when a single nucleus is followed in time (e.g., through molecular dynamics simulations), ζ is a random function of t. Although this second picture of random trajectories was invoked when writing the dipolar correlation function in Eq. (2), in the following pages we only work with the probabilistic description of an ensemble of identical nuclei.

Below, we will use the dynamical rule (Eq. 29) when combining the stochastic dynamics of ζ with the spin dynamics from Sect. 2. The combined dynamics will be described by a stochastic Liouville equation (SLE) for a ζ-conditioned spin variable. In the case of the nuclear polarization, for example, the SLE will describe the dynamics of P_I(t), which stands for P_I(ζ,t) in the continuous case. For a detailed explanation of SLE the reader is referred to the literature (Kubo, 1969; Gamliel and Levanon, 1995). A more recent discussion can be found in Kuprov (2016).

3.2 Rate equations in liquids

3.2.1 Electronic polarization

The electronic polarization was assumed to be insensitive to the dipolar coupling with the nuclear spins. Hence, the ζ-conditioned electronic polarization P_S(t) is of the following separable form:

$$\begin{array}{}\text{(30)}& {\mathit{P}}_S\left(t\right)={\mathit{p}}^{\mathrm{eq}}{P}_S\left(t\right),\end{array}$$

in which all ζ dependence is isolated in the equilibrium probability of the stochastic process. From P_S(t) we obtain the averaged (over ζ) electronic polarization by summing/integrating over the ζ dependence. This is done with the help of the constant vector/function 1 as follows:

$$\begin{array}{}\text{(31)}& {\mathbf{1}}^{\mathrm{T}}{\mathit{P}}_S\left(t\right)={\mathbf{1}}^{\mathrm{T}}{\mathit{p}}^{\mathrm{eq}}{P}_S\left(t\right)=P_S\left(t\right).\end{array}$$

(In the last equality we used the normalization of the probability, 1^Tp^eq=1, which reads $\int \mathrm{d}\mathit{\zeta }{p}^{\mathrm{eq}}\left(\mathit{\zeta }\right)=\mathrm{1}$ in the continuous case.) Note that P_S(t) in Eq. (30) is, in fact, the electronic polarization averaged over the stochastic variable (Eq. 31).

In this description, the experimentally accessible polarizations correspond to the averaged values, while the ζ-dependent variables, like P_S(t), serve only an intermediate, book-keeping role. In other words, at the end we will always average over ζ by using the constant vector/function 1.

Since Eq. (30) holds in general for the electronic polarization, it also holds at steady state and at equilibrium:

$$\begin{array}{}\text{(32)}& {\mathit{P}}_S^{\mathrm{ss}}={\mathit{p}}^{\mathrm{eq}}{P}_S^{\mathrm{ss}},{\mathit{P}}_S^{\mathrm{eq}}={\mathit{p}}^{\mathrm{eq}}{P}_S^{\mathrm{eq}}.\end{array}$$

Here $P_S^{\mathrm{ss}}$ and $P_S^{\mathrm{eq}}$ are the averaged (over ζ) values which correspond to the macroscopic polarization.

Lastly, we point out that at equilibrium all joint spin-ζ properties are of the above separable form. In other words, the last equality in Eq. (32) is not limited to the electronic polarization but applies to all other equilibrium properties.

3.2.2 Nuclear polarization

To illustrate the SLE formalism and to introduce further notation, we start by transforming the rate equation of the nuclear polarization (Eq. 6) to an SLE:

$$\begin{array}{}\text{(33)}& \begin{array}{rl}{\dot{\mathit{P}}}_I\left(t\right)& =-\mathbf{K}{\mathit{P}}_I\left(t\right)-{\mathbb{R}}_{\mathrm{1}I}\left({\mathit{P}}_I\left(t\right)-{\mathit{p}}^{\mathrm{eq}}{P}_I^{\mathrm{eq}}\right)\\ & -{\mathbb{V}}_{+}{\mathit{P}}_I\left(t\right)-{\mathbb{V}}_{-}{\mathit{p}}^{\mathrm{eq}}{P}_S\left(t\right).\end{array}\end{array}$$

There are several different things going on here, so let us examine them one by one.

First, following the convention introduced above, P_I(t) stands for P_I(ζ,t), which is the nuclear polarization conditional on the random state ζ. In this case the dot indicates partial derivative with respect to the time dependence, at fixed ζ. Second, the term KP_I drives the dynamics in the ζ space by providing “off-diagonal” elements that mix the different random states. All remaining terms on the right-hand side of the SLE are “diagonal” in the ζ space and act only on the spin degree(s) of freedom (which are conditioned on ζ). Third, the mw excitation rates v_± and the relaxation rate R_1I have acquired ζ dependence, turning into operators in ζ space that act on P_I(t) or p^eq. In the discrete case, these would be matrices with different v_± and R_1I values for each discrete state ζ along their main diagonals. We use hollow capital letters to denote such “diagonal” operators in ζ space, also in the continuous case. Fourth, as all equilibrium properties, the nuclear Boltzmann polarization is separable, with the ζ dependence confined to the equilibrium probability of the random process.

The steady state of Eq. (33) is

$$\begin{array}{}\text{(34)}& \left(\mathbf{K}+{\mathbb{R}}_{\mathrm{1}I}+{\mathbb{V}}_{+}\right){\mathit{P}}_I^{\mathrm{ss}}=\left(P_I^{\mathrm{eq}}{\mathbb{R}}_{\mathrm{1}I}-P_S^{\mathrm{eq}}p{\mathbb{V}}_{-}\right){\mathit{p}}^{\mathrm{eq}},\end{array}$$

where we used $P_S^{\mathrm{ss}}=p{P}_S^{\mathrm{eq}}$ (Eq. 5). Our aim is to solve Eq. (34) for ${\mathit{P}}_I^{\mathrm{ss}}$ and then obtain the macroscopic nuclear polarization by calculating the average $P_I^{\mathrm{ss}}={\mathbf{1}}^{\mathrm{T}}{\mathit{P}}_I^{\mathrm{ss}}$.

Clearly, solving Eq. (34) consists of calculating the inverse of the operator $(\mathbf{K}+{\mathbb{R}}_{\mathrm{1}I}+{\mathbb{V}}_{+})$. This is a daunting task in general and requires the matrix representation of K in some basis set. Here we will limit the discussion to random motions that are orders of magnitude faster than the nuclear T₁ relaxation rate, which we concluded to be the case even in viscous liquids like glycerol. This assumption ensures that, at steady state, all nuclear spins in the sample are equivalent and have the same polarization. Hence, we will look for a separable steady-state solution of the form

$$\begin{array}{}\text{(35)}& {\mathit{P}}_I^{\mathrm{ss}}={\mathit{p}}^{\mathrm{eq}}{P}_I^{\mathrm{ss}}\left(\text{ansatz for liquids}\right).\end{array}$$

With this ansatz, Eq. (34) becomes

$$\begin{array}{}\text{(36)}& \left(\mathbf{K}+{\mathbb{R}}_{\mathrm{1}I}+{\mathbb{V}}_{+}\right){\mathit{p}}^{\mathrm{eq}}{P}_I^{\mathrm{ss}}=\left(P_I^{\mathrm{eq}}{\mathbb{R}}_{\mathrm{1}I}-P_S^{\mathrm{eq}}p{\mathbb{V}}_{-}\right){\mathit{p}}^{\mathrm{eq}}.\end{array}$$

While the difference between Eqs. (34) and (36) appears to be minor, in fact we have achieved a tremendous simplification since Kp^eq=0, and thus the dynamical aspect of the random process is gone; only its equilibrium (i.e., time-independent) properties remain. Indeed, since in Eq. (36) all ζ operators act on the equilibrium probability p^eq, integration over the ζ dependence brings the average values:

$$\begin{array}{}\text{(37)}& \langle {\mathbb{R}}_{\mathrm{1}I}\rangle ={\mathbf{1}}^{\mathrm{T}}{\mathbb{R}}_{\mathrm{1}I}{\mathit{p}}^{\mathrm{eq}},\langle {\mathbb{V}}_{\pm }\rangle ={\mathbf{1}}^{\mathrm{T}}{\mathbb{V}}_{\pm }{\mathit{p}}^{\mathrm{eq}}.\end{array}$$

(These are the average values of the ζ-dependent functions R_1I(ζ) and v_±(ζ), respectively.) After averaging, Eq. (36) becomes

$$\begin{array}{}\text{(38)}& \left(\langle {\mathbb{R}}_{\mathrm{1}I}\rangle +\langle {\mathbb{V}}_{+}\rangle \right)P_I^{\mathrm{ss}}=\langle {\mathbb{R}}_{\mathrm{1}I}\rangle P_I^{\mathrm{eq}}-p\langle {\mathbb{V}}_{-}\rangle P_S^{\mathrm{eq}}.\end{array}$$

Comparison of Eqs. (38) and (7) shows that the phenomenological rates R_1I and v_± in the rate equation should be identified with the macroscopic averages 〈R_1I〉 and 〈𝕍_±〉 over the liquid sample. This is the familiar regime of fast motion, where one observes the averaged values of the magnetic parameters. We have thus provided a formal justification of why the averaged δ² in Eq. (27) corresponds to fast spin diffusion in the case of solids.

We observe that the static averages over 3D space, which are implied by Eq. (37), do not allow for the partial dynamical averaging of the dipolar interaction. As discussed above, such averaging should be based on the time-correlation function of A₁. However the rate constants v_± always contain the square of the dipolar interaction and do not provide access to A₁ itself. We thus conclude that the partial averaging of the dipolar interaction in liquids is inaccessible to modeling by rate equations.

3.3 Spin dynamics in liquids

To gain access to the dipolar interaction before it is squared, we turn to the equations of motion of the coherences from Sect. 2.3. We first transform the equation of i_z (Eq. 14) to an SLE:

$$\begin{array}{}\text{(39)}& {\dot{\mathit{i}}}_z\left(t\right){|}_{\mathrm{coh}}=-\mathbf{K}{\mathit{i}}_z\left(t\right)-\mathrm{Re}\mathit{\left\{}i{\mathbb{A}}_{\mathrm{1}}^{*}{\mathit{g}}_z\right(t\left)\mathit{\right\}}.\end{array}$$

As before, i_z(t) and g_z(t) stand for i_z(ζ,t) and g_z(ζ,t) in the continuous case, K acts on the ζ dependence of i_z and $A_{\mathrm{1}}^{*}$ has become a “diagonal” operator in ζ space. Averaging Eq. (39) over ζ, we obtain the macroscopic equation

$$\begin{array}{}\text{(40)}& {\mathbf{1}}^{\mathrm{T}}{\dot{\mathit{i}}}_z\left(t\right){|}_{\mathrm{coh}}=-\mathrm{Re}\mathit{\left\{}i{\mathbf{1}}^{\mathrm{T}}{\mathbb{A}}_{\mathrm{1}}^{*}{\mathit{g}}_z\right(t\left)\mathit{\right\}}.\end{array}$$

At steady state, using the ansatz for liquids (Eq. 35) in the form ${\mathit{i}}_z^{\mathrm{ss}}={\mathit{p}}^{\mathrm{eq}}{i}_z^{\mathrm{ss}}$, we have

$$\begin{array}{}\text{(41)}& {\dot{i}}_z{|}_{\mathrm{coh}}^{\mathrm{ss}}=-\mathrm{Re}\mathit{\left\{}i{\mathbf{1}}^{\mathrm{T}}{\mathbb{A}}_{\mathrm{1}}^{*}{\mathit{g}}_z^{\mathrm{ss}}\mathit{\right\}}.\end{array}$$

Since we accounted only for the coherent contribution to the time derivative of i_z, the right-hand side of Eq. (41) corresponds to the right-hand side of Eq. (7). Our aim is to identify the phenomenological rate constants v_± that should be used in Eq. (7) by analyzing $\mathrm{Re}\mathit{\left\{}i{\mathbf{1}}^{\mathrm{T}}{\mathbb{A}}_{\mathrm{1}}^{*}{\mathit{g}}_z^{\mathrm{ss}}\mathit{\right\}}$. However, because the random modulation of A₁ additionally contributes to the nuclear T₁ relaxation, we have the equality

$$\begin{array}{}\text{(42)}& \mathrm{Re}\mathit{\left\{}i{\mathbf{1}}^{\mathrm{T}}{\mathbb{A}}_{\mathrm{1}}^{*}{\mathit{g}}_z^{\mathrm{ss}}\mathit{\right\}}=R_{\mathrm{1}I}^A{i}_z^{\mathrm{ss}}+v_{+}{i}_z^{\mathrm{ss}}+p{v}_{-}{s}_z^{\mathrm{eq}}\end{array}$$

from which we will read out the desired rates.

3.3.1 Contribution to nuclear T₁ relaxation

The contribution of $\mathrm{Re}\mathit{\left\{}i{\mathbf{1}}^{\mathrm{T}}{\mathbb{A}}_{\mathrm{1}}^{*}{\mathit{g}}_z^{\mathrm{ss}}\mathit{\right\}}$ to the nuclear spin-lattice relaxation can be identified by its value in the absence of mw irradiation. From Eqs. (12) and (13) we see that for ω₁=0 the dynamics of g_z completely decouples from g_x and g_y. The SLE of g_z in this case becomes

$$\begin{array}{}\text{(43)}& {\dot{\mathit{g}}}_z\left(t\right)=-\mathbf{K}{\mathit{g}}_z\left(t\right)-(R_{\mathrm{1}S}+i{\mathit{\omega }}_I){\mathit{g}}_z\left(t\right)-i{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\mathbb{A}}_{\mathrm{1}}{\mathit{i}}_z\left(t\right).\end{array}$$

Technically, R_1S and ω_I should be operators that act on the ζ dependence of g_z. However, we take the electronic T₁ relaxation rate and the nuclear Larmor frequency to be independent of the dipolar coupling, which is parametrized by ζ. The corresponding operators are then R_1S𝕀 and ω_I𝕀, where 𝕀 is the identity operator in ζ space. This identity operator will not be written explicitly.

The steady-state solution of Eq. (43) is

$$\begin{array}{}\text{(44)}& {\mathit{g}}_z^{\mathrm{ss}}=-i{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}(\mathbf{K}+R_{\mathrm{1}S}+i{\mathit{\omega }}_I{)}^{-\mathrm{1}}{\mathbb{A}}_{\mathrm{1}}{\mathit{p}}^{\mathrm{eq}}{i}_z^{\mathrm{ss}},\end{array}$$

where we again used the ansatz for liquids. Substituting this ${\mathit{g}}_z^{\mathrm{ss}}$ on the left-hand side of Eq. (42), we find that the nuclear relaxation rate due to A₁ is

$$\begin{array}{}\text{(45)}& R_{\mathrm{1}I}^A={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\mathrm{Re}\left\{{\mathbf{1}}^{\mathrm{T}}{\mathbb{A}}_{\mathrm{1}}^{*}{\left(\mathbf{K}+R_{\mathrm{1}S}+i{\mathit{\omega }}_I\right)}^{-\mathrm{1}}{\mathbb{A}}_{\mathrm{1}}{\mathit{p}}^{\mathrm{eq}}\right\}.\end{array}$$

To express this relaxation rate in a more intelligible manner, we observe that the inverse of a matrix M whose eigenvalues have strictly positive real parts can be written as

$$\begin{array}{}\text{(46)}& M^{-\mathrm{1}}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{e}^{-Mt}\mathrm{d}t.\end{array}$$

Applying this identity to the matrix $(\mathbf{K}+R_{\mathrm{1}S}+i{\mathit{\omega }}_I)$ in Eq. (45), we find

$$\begin{array}{}\text{(47)}& R_{\mathrm{1}I}^A={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\mathrm{Re}\left\{\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\mathrm{d}t{e}^{-(R_{\mathrm{1}S}+i{\mathit{\omega }}_I)t}{C}_{\mathrm{11}}\left(t\right)\right\},\end{array}$$

where

$$\begin{array}{}\text{(48)}& \begin{array}{rl}{C}_{\mathrm{11}}\left(t\right)& ={\mathbf{1}}^{\mathrm{T}}{\mathbb{A}}_{\mathrm{1}}^{*}{e}^{-\mathbf{K}t}{\mathbb{A}}_{\mathrm{1}}{\mathit{p}}^{\mathrm{eq}}\\ & =\int \mathrm{d}\mathit{\zeta }{A}_{\mathrm{1}}^{*}\left(\mathit{\zeta }\right)e^{-K_{\mathit{\zeta }}t}{A}_{\mathrm{1}}\left(\mathit{\zeta }\right)p^{\mathrm{eq}}\left(\mathit{\zeta }\right)\end{array}\end{array}$$

is the time-correlation function of the dipolar interaction (Eq. 2). Since the integral in Eq. (47) corresponds to the Laplace transform

$$\begin{array}{}\text{(49)}& J_{\mathrm{11}}\left(s\right)=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\mathrm{d}t{e}^{-st}{C}_{\mathrm{11}}\left(t\right),\end{array}$$

we have

$$\begin{array}{}\text{(50)}& R_{\mathrm{1}I}^A={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\mathrm{Re}\left\{J_{\mathrm{11}}(R_{\mathrm{1}S}+i{\mathit{\omega }}_I)\right\}.\end{array}$$

The real part of the Laplace transform is known as spectral density. Here the spectral density is evaluated at a the complex argument R_1S+iω_I, which contains both the nuclear Larmor frequency and the electronic T₁ relaxation rate.

Let us examine Eq. (50) in the solid-state limit where A₁ does not change with time. Then $C_{\mathrm{11}}=\langle A_{\mathrm{1}}^{*}{A}_{\mathrm{1}}\rangle$ and

$$\begin{array}{}\text{(51)}& R_{\mathrm{1}I,\mathrm{solid}}^A=\langle {\mathit{\delta }}^{\mathrm{2}}\rangle {\displaystyle \frac{T_{\mathrm{1}S}}{\mathrm{1}+T_{\mathrm{1}S}^{\mathrm{2}}{\mathit{\omega }}_I^{\mathrm{2}}}}.\end{array}$$

In Abragam's nomenclature (Abragam, 1961) this is relaxation of the second kind, meaning that it is due to the relaxation of the electronic spins and not due to the modulation of the dipolar interaction by motion. This $R_{\mathrm{1}I,\mathrm{solid}}^A/\langle {\mathit{\delta }}^{\mathrm{2}}\rangle$ was shown with horizontal, dotted red lines in the third row of Fig. 3.

In the case of liquids, we expect C₁₁(t) to decay with time. Assuming a mono-exponential decay with correlation time τ,

$$\begin{array}{}\text{(52)}& \begin{array}{rl}{R}_{\mathrm{1}I,\exp }^A& =\langle {\mathit{\delta }}^{\mathrm{2}}\rangle \mathrm{Re}\left\{\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\mathrm{d}t{e}^{-\left(R_{\mathrm{1}S}+i{\mathit{\omega }}_I+{\mathit{\tau }}^{-\mathrm{1}}\right)t}\right\}\\ & =\langle {\mathit{\delta }}^{\mathrm{2}}\rangle \mathrm{Re}\left\{{\left(R_{\mathrm{1}S}+{\mathit{\tau }}^{-\mathrm{1}}+i{\mathit{\omega }}_I\right)}^{-\mathrm{1}}\right\}.\end{array}\end{array}$$

This $R_{\mathrm{1}I,\exp }^A/\langle {\mathit{\delta }}^{\mathrm{2}}\rangle$ was shown with horizontal, dotted black lines in the third row of Fig. 3. To understand why it increases with decreasing ω_I, let us examine the case of motion that is faster than the electronic T₁ time, i.e., ${\mathit{\tau }}^{-\mathrm{1}}\gg R_{\mathrm{1}S}$. The result,

$$\begin{array}{}\text{(53)}& R_{\mathrm{1}I,\exp }^A\approx \langle {\mathit{\delta }}^{\mathrm{2}}\rangle {\displaystyle \frac{\mathit{\tau }}{\mathrm{1}+{\mathit{\tau }}^{\mathrm{2}}{\mathit{\omega }}_I^{\mathrm{2}}}},\end{array}$$

is relaxation of the first kind with Lorentzian spectral density. Clearly, smaller ω_I implies larger dipolar contribution to the nuclear T₁ relaxation rate.

Having identified the relaxation rate $R_{\mathrm{1}I}^A$ on the right-hand side of Eq. (42), we now proceed with the analysis of the rates v_± characterizing the forbidden transitions.

3.3.2 Contribution to forbidden transitions

Combining the dynamics of the coherences (Eq. 12) with the stochastic dynamics (Eq. 29), we arrive at the following SLE:

$$\begin{array}{}\text{(54)}& \left[\begin{array}{c}{\dot{\mathit{g}}}_x\\ {\dot{\mathit{g}}}_y\\ {\dot{\mathit{g}}}_z\end{array}\right]=-\left(\mathbf{K}+\mathfrak{B}\right)\left[\begin{array}{c}{\mathit{g}}_x\\ {\mathit{g}}_y\\ {\mathit{g}}_z\end{array}\right]-i{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\mathbb{A}}_{\mathrm{1}}\left[\begin{array}{c}-i{s}_y{\mathit{p}}^{\mathrm{eq}}\\ i{s}_x{\mathit{p}}^{\mathrm{eq}}\\ {\mathit{i}}_z\end{array}\right].\end{array}$$

Note that the matrix 𝔅 does not depend on time as the electronic relaxation properties were taken to be insensitive to the dipolar interaction between the electronic and nuclear spins. Recall that the operators written as upright bold letters (including the hollow ones) act on the ζ dependence of the variables, which is encoded by the italic bold symbols. The script uppercase letters denote 3×3 matrices, which act on the column vectors that are shown explicitly.

Although not shown explicitly in Eq. (54), we imply the tensor products of the operators K, 𝔅 and 𝔸₁ with the identity operators in the spaces on which K, 𝔅 and 𝔸₁ do not act (i.e., K and 𝔸₁ are in fact ℐ⊗K and ℐ⊗𝔸₁ where ℐ is the 3×3 identity matrix, and 𝔅 is 𝔅⊗𝕀 where 𝕀 is the identity operator in the ζ space).

Before solving Eq. (54) at steady state, let us introduce the (right) eigenvalue problem of 𝔅,

$$\begin{array}{}\text{(55)}& \mathfrak{B}\mathfrak{U}=\mathfrak{U}\mathcal{D},\end{array}$$

where the diagonal matrix

$$\begin{array}{}\text{(56)}& \mathcal{D}=\text{diag}\left({\mathit{\lambda }}_{\mathrm{1}},{\mathit{\lambda }}_{\mathrm{2}},{\mathit{\lambda }}_{\mathrm{3}}\right)\end{array}$$

contains the eigenvalues of 𝔅 along its main diagonal, and the columns of the 3×3 matrix 𝔘 contain the corresponding right eigenvectors. Then, the steady state of Eq. (54) is

$$\begin{array}{}\text{(57)}& \mathfrak{U}\left(\mathbf{K}+\mathcal{D}\right){\mathfrak{U}}^{-\mathrm{1}}\left[\begin{array}{c}{\mathit{g}}_x^{\mathrm{ss}}\\ {\mathit{g}}_y^{\mathrm{ss}}\\ {\mathit{g}}_z^{\mathrm{ss}}\end{array}\right]=-i{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\mathbb{A}}_{\mathrm{1}}{\mathit{p}}^{\mathrm{eq}}\left[\begin{array}{c}-i{s}_y^{\mathrm{ss}}\\ i{s}_x^{\mathrm{ss}}\\ i_z^{\mathrm{ss}}\end{array}\right],\end{array}$$

which, after inverting the matrices, yields

$$\begin{array}{}\text{(58)}& \left[\begin{array}{c}{\mathit{g}}_x^{\mathrm{ss}}\\ {\mathit{g}}_y^{\mathrm{ss}}\\ {\mathit{g}}_z^{\mathrm{ss}}\end{array}\right]=-i{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\mathfrak{U}{\left(\mathbf{K}+\mathcal{D}\right)}^{-\mathrm{1}}{\mathbb{A}}_{\mathrm{1}}{\mathit{p}}^{\mathrm{eq}}{\mathfrak{U}}^{-\mathrm{1}}\left[\begin{array}{c}-i{s}_y^{\mathrm{ss}}\\ i{s}_x^{\mathrm{ss}}\\ i_z^{\mathrm{ss}}\end{array}\right].\end{array}$$

Plugging this solution for ${\mathit{g}}_z^{\mathrm{ss}}$ into the left-hand side of Eq. (42), and defining the matrix

$$\begin{array}{}\text{(59)}& \mathcal{L}={\mathbf{1}}^{\mathrm{T}}{\mathbb{A}}_{\mathrm{1}}^{*}{\left(\mathbf{K}+\mathcal{D}\right)}^{-\mathrm{1}}{\mathbb{A}}_{\mathrm{1}}{\mathit{p}}^{\mathrm{eq}},\end{array}$$

we find

$$\begin{array}{}\text{(60)}& \mathrm{Re}\mathit{\left\{}i{\mathbf{1}}^{\mathrm{T}}{\mathbb{A}}_{\mathrm{1}}^{*}{\mathit{g}}_z^{\mathrm{ss}}\mathit{\right\}}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\mathrm{Re}\mathit{\left\{}\left[\begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]\mathfrak{U}\mathcal{L}{\mathfrak{U}}^{-\mathrm{1}}\left[\begin{array}{c}-i{s}_y^{\mathrm{ss}}\\ i{s}_x^{\mathrm{ss}}\\ i_z^{\mathrm{ss}}\end{array}\right]\mathit{\right\}}.\end{array}$$

Comparison with the right-hand side of Eq. (42) yields

$$\begin{array}{}\text{(61)}& \begin{array}{rl}& R_{\mathrm{1}I}^A+v_{+}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left[\begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]\mathrm{Re}\left\{\mathfrak{U}\mathcal{L}{\mathfrak{U}}^{-\mathrm{1}}\right\}{\left[\begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]}^{\mathrm{T}}\\ & v_{-}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left[\begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]\mathrm{Re}\left\{i\mathfrak{U}\mathcal{L}{\mathfrak{U}}^{-\mathrm{1}}\right\}{\left[\begin{array}{ccc}{\mathit{\omega }}_{\mathrm{1}}{f}_y& {\mathit{\omega }}_{\mathrm{1}}{f}_x& \mathrm{0}\end{array}\right]}^{\mathrm{T}},\end{array}\end{array}$$

where we used the relationships between $s_{x,y}^{\mathrm{ss}}$ and $s_z^{\mathrm{eq}}$ (Eq. 15) to arrive at v₋.

We observe that ℒ is a 3×3 diagonal matrix without any ζ dependence since the right-hand side of Eq. (59) is averaged over ζ. With $\mathcal{L}=\text{diag}(L_{\mathrm{1}},L_{\mathrm{2}},L_{\mathrm{3}})$, we have

$$\begin{array}{}\text{(62)}& L_n={\mathbf{1}}^{\mathrm{T}}{\mathbb{A}}_{\mathrm{1}}^{*}{\left(\mathbf{K}+{\mathit{\lambda }}_n\right)}^{-\mathrm{1}}{\mathbb{A}}_{\mathrm{1}}{\mathit{p}}^{\mathrm{eq}}(n=\mathrm{1},\mathrm{2},\mathrm{3}).\end{array}$$

Using Eq. (46), these diagonal elements can be written as

$$\begin{array}{}\text{(63)}& \begin{array}{rl}{L}_n& ={\mathbf{1}}^{\mathrm{T}}{\mathbb{A}}_{\mathrm{1}}^{*}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{e}^{-\left(\mathbf{K}+{\mathit{\lambda }}_n\right)t}\mathrm{d}t{\mathbb{A}}_{\mathrm{1}}{\mathit{p}}^{\mathrm{eq}}\\ & =\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{e}^{-{\mathit{\lambda }}_{n}t}\left({\mathbf{1}}^{\mathrm{T}}{\mathbb{A}}_{\mathrm{1}}^{*}{e}^{-\mathbf{K}t}{\mathbb{A}}_{\mathrm{1}}{\mathit{p}}^{\mathrm{eq}}\right)\mathrm{d}t\\ & =\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{e}^{-{\mathit{\lambda }}_{n}t}{C}_{\mathrm{11}}\left(t\right)\mathrm{d}t=J_{\mathrm{11}}\left({\mathit{\lambda }}_n\right),\end{array}\end{array}$$

where we used Eq. (48) in the third equality and Eq. (49) in the last one. Hence, each L_n is the Laplace transform of the time-correlation function C₁₁(t) evaluated at the eigenvalue λ_n of 𝔅. The matrix ℒ to be used in Eq. (61) is thus

$$\begin{array}{}\text{(64)}& \mathcal{L}=\text{diag}\left(J_{\mathrm{11}}\left({\mathit{\lambda }}_{\mathrm{1}}\right),J_{\mathrm{11}}\left({\mathit{\lambda }}_{\mathrm{2}}\right),J_{\mathrm{11}}\left({\mathit{\lambda }}_{\mathrm{3}}\right)\right).\end{array}$$

In summary, for any given set of parameters, we form the 3×3 matrix 𝔅 (Eq. 13) and numerically calculate its eigenvalues and right eigenvectors. The former are used in Eq. (64) to calculate ℒ. Sandwiching ℒ by the eigenvectors, as required in Eq. (61), we arrive at the desired rates v_±. This prescription applies to any motional model describing the stochastic dynamics of the inter-spin vector. Different models will differ only in their spectral densities J₁₁.

From a mathematical point of view, the simplest case is a model with exponential dipolar correlation function, $C_{\mathrm{11}}^{\exp }\left(t\right)=\langle {\mathit{\delta }}^{\mathrm{2}}\rangle e^{-t/\mathit{\tau }}$, where τ is the correlation time. Then

$$\begin{array}{}\text{(65)}& J_{\mathrm{11}}^{\exp }\left(s\right)=〈{\mathit{\delta }}^{\mathrm{2}}〉{\displaystyle \frac{\mathrm{1}}{s+\frac{\mathrm{1}}{\mathit{\tau }}}}=〈{\mathit{\delta }}^{\mathrm{2}}〉\mathit{\tau }{\displaystyle \frac{\mathrm{1}}{s\mathit{\tau }+\mathrm{1}}},\end{array}$$

and Eq. (64) becomes

$$\begin{array}{}\text{(66)}& {\mathcal{L}}_{\exp }=〈{\mathit{\delta }}^{\mathrm{2}}〉\text{diag}\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{1}}+\frac{\mathrm{1}}{\mathit{\tau }}}},{\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{2}}+\frac{\mathrm{1}}{\mathit{\tau }}}},{\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{3}}+\frac{\mathrm{1}}{\mathit{\tau }}}}\right).\end{array}$$

All dashed black lines labeled “liquid” in Fig. 3 were calculated using Eq. (66) with τ=12 ns.

Comparing v₊ and pv₋ (Fig. 3, third and fourth rows) between the solid and liquid cases, we see that at Q, W and J bands, the fluctuations of the dipolar interaction have substantially broadened the lines centered at the canonical solid-effect offsets $\mathrm{\Omega }\approx \pm {\mathit{\omega }}_I$ and have reduced the peak enhancements in liquids compared to solids (fifth row). At X band, where the two lines had already merged in the solid case, the effect of fluctuations is qualitatively different, although line broadening is also visible. Most strikingly, the rate v₊ is seen to become negative at offsets larger than ω_I, which leads to a nuclear polarization factor (Eq. 9) that exceeds 1 (Fig. 3, bottom row).

4 Closer look at liquids

4.1 Translational diffusion of hard spheres

A mono-exponential dipolar correlation function is a poor model of translational diffusion in liquids. The so-called force-free hard-sphere (FFHS) model, which assumes spherical molecules that contain the spins at their centers, is a more realistic yet analytically tractable model (Ayant et al., 1975; Hwang and Freed, 1975). It is universally employed in the analysis of diverse magnetic-resonance measurements, including nuclear relaxation by paramagnetic impurities (Okuno et al., 2022) and DNP via the Overhauser effect (Franck et al., 2013).

Because the spins are taken to be at the centers of the spherical molecules, the FFHS model has only two parameters: the coefficient of translational diffusion, D, and the distance of the spins upon contact of the spherical molecules, b. These two parameters form the characteristic motional timescale of the model (Ayant et al., 1975):

$$\begin{array}{}\text{(67)}& \mathit{\tau }=b^{\mathrm{2}}/D.\end{array}$$

The Laplace transform of the dipolar correlation function of this model is (Ayant et al., 1975, Eqs. 51 and 55)

$$\begin{array}{}\text{(68)}& J_{\mathrm{11}}^{\mathrm{ffhs}}\left(s\right)=\langle {\mathit{\delta }}^{\mathrm{2}}\rangle \mathit{\tau }{\displaystyle \frac{(s\mathit{\tau }{)}^{\frac{\mathrm{1}}{\mathrm{2}}}+\mathrm{4}}{(s\mathit{\tau }{)}^{\frac{\mathrm{3}}{\mathrm{2}}}+\mathrm{4}(s\mathit{\tau })+\mathrm{9}(s\mathit{\tau }{)}^{\frac{\mathrm{1}}{\mathrm{2}}}+\mathrm{9}}}.\end{array}$$

Using $J_{\mathrm{11}}^{\mathrm{ffhs}}$ in Eq. (64) we calculated numerically the same properties as in Fig. 3 but for the FFHS model. The results are shown with colored solid lines in Fig. 4. For comparison, the model with mono-exponential correlation function from Fig. 3 is also reproduced in Fig. 4 with dashed black lines.

Figure 4 Same as Fig. 3 for the model with exponential time-correlation function (dashed black lines) and the FFHS model (colored solid lines) both with $\mathit{\tau }=T_{\mathrm{2}S}/\mathrm{5}=\mathrm{12}$ ns.

The general observation from Fig. 3 that the fluctuation of the dipolar interaction broadens the solid-effect lines at $\mathrm{\Omega }\approx \pm {\mathit{\omega }}_I$ is even more relevant for the FFHS model. Indeed, for the same dipolar timescale τ, the FFHS lines are much broader and, correspondingly, much smaller in peak amplitude than the lines of the exponential model. Hence, the FFHS model predicts significantly smaller DNP enhancements (Fig. 4, second last row) compared to the exponential model with the same timescale τ. At X band, the negative values of v₊ are still present, but their magnitude is substantially reduced (third row). The corresponding offsets where the nuclear polarization factor, p_X, is larger than 1 are similar in the two models, but again the deviation from 1 is much smaller in the FFHS model (last row).

Overall, p_X in liquids is very close to 1 (last row of Fig. 4, FFHS model), which indicates that v₊ is very small compared to R_1I. In such cases, the solid-effect DNP enhancement (Eq. 9) is well approximated by Eq. (10). This explains why the enhancement in the fifth row of Fig. 4 is essentially a rescaled version of the row directly above it.

The substantial reduction of the peak intensities at the solid-effect offsets $\mathrm{\Omega }\approx \pm {\mathit{\omega }}_I$ is accompanied by a smaller but still appreciable increase of the intensities at small offsets (Ω≈0). This trend is visible both in the transition from the solid case to a mono-exponential correlation function (Fig. 3, second row) and in the further transition to the FFHS model (Fig. 4). The significance of this observation will become clear in Sect. 4.4, where we compare our calculations with the experiments of Kuzhelev et al. (2022).

4.2 Approximate matrix inversion

Since 𝔅 is a 3×3 matrix, its eigenvalues and eigenvectors are easily determined numerically, as we did when calculating the exponential and FFHS models in Fig. 4. Nevertheless, to gain insight into the eigenvalue problem that is being solved, here we analyze Eq. (55) using perturbation theory. The analysis reveals that the eigenvalue problem is related to the effective magnetic field and the associated “tilted” coordinate frame (Wenckebach, 2016).

Let us introduce the matrix

$$\begin{array}{}\text{(69)}& {\mathfrak{B}}_{\mathrm{0}}=\left[\begin{array}{ccc}{R}_{\mathrm{2}S}+i{\mathit{\omega }}_I& \mathrm{\Omega }& \mathrm{0}\\ -\mathrm{\Omega }& R_{\mathrm{2}S}+i{\mathit{\omega }}_I& {\mathit{\omega }}_{\mathrm{1}}\\ \mathrm{0}& -{\mathit{\omega }}_{\mathrm{1}}& R_{\mathrm{2}S}+i{\mathit{\omega }}_I\end{array}\right],\end{array}$$

where R_1S in the lower right corner of 𝔅 (Eq. 13) has been replaced by R_2S. The three eigenvalues of 𝔅₀ are

$$\begin{array}{}\text{(70)}& {\mathit{\lambda }}_{\mathrm{0},\mathrm{0}}=R_{\mathrm{2}S}+i{\mathit{\omega }}_I,{\mathit{\lambda }}_{\mathrm{0},\mp }=R_{\mathrm{2}S}+i\left({\mathit{\omega }}_I\mp {\mathit{\omega }}_{\mathrm{eff}}\right),\end{array}$$

where the frequency

$$\begin{array}{}\text{(71)}& {\mathit{\omega }}_{\mathrm{eff}}=\sqrt{{\mathrm{\Omega }}^{\mathrm{2}}+{\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}}\end{array}$$

corresponds to the effective magnetic field in the rotating frame. This field is tilted away from the z axis by an angle α such that

$$\begin{array}{}\text{(72)}& \cos \mathit{\alpha }=\mathrm{\Omega }/{\mathit{\omega }}_{\mathrm{eff}}=\mathrm{c},\sin \mathit{\alpha }={\mathit{\omega }}_{\mathrm{1}}/{\mathit{\omega }}_{\mathrm{eff}}=\mathrm{s}.\end{array}$$

With the sine and cosine of α abbreviated as s and c, the right eigenvectors of 𝔅₀ are

$$\begin{array}{}\text{(73)}& {\mathfrak{U}}_{\mathrm{0}}=\left[\begin{array}{ccc}\mathrm{s}& -\mathrm{c}/\sqrt{\mathrm{2}}& -\mathrm{c}/\sqrt{\mathrm{2}}\\ \mathrm{0}& i/\sqrt{\mathrm{2}}& -i/\sqrt{\mathrm{2}}\\ \mathrm{c}& \mathrm{s}/\sqrt{\mathrm{2}}& \mathrm{s}/\sqrt{\mathrm{2}}\end{array}\right],\end{array}$$

where the first column corresponds to λ_0,0, the second to ${\mathit{\lambda }}_{\mathrm{0},-}$, and the third to ${\mathit{\lambda }}_{\mathrm{0},+}$. By inspection, ${\mathfrak{U}}_{\mathrm{0}}^{-\mathrm{1}}={\mathfrak{U}}_{\mathrm{0}}^{\mathrm{H}}$, where the superscript “H” denotes Hermitian conjugation.

We treat the difference 𝔅−𝔅₀ as a perturbation to 𝔅₀. To first order in the perturbation, the eigenvalues of the original matrix 𝔅 are ${\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_n={\mathit{\lambda }}_{\mathrm{0},n}+u_{\mathrm{0},n}^{\mathrm{H}}(\mathfrak{B}-{\mathfrak{B}}_{\mathrm{0}})u_{\mathrm{0},n}$, where u_0,n is the nth column of 𝔘₀. Using this expression we find the corrected eigenvalues

$$\begin{array}{}\text{(74)}& {\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{\mathrm{0}}={\stackrel{\mathrm{̃}}{R}}_{\mathrm{1}}+i{\mathit{\omega }}_I,{\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{\mp }={\stackrel{\mathrm{̃}}{R}}_{\mathrm{2}}+i\left({\mathit{\omega }}_I\mp {\mathit{\omega }}_{\mathrm{eff}}\right),\end{array}$$

with

$$\begin{array}{}\text{(75)}& \begin{array}{rl}{\stackrel{\mathrm{̃}}{R}}_{\mathrm{1}}& =R_{\mathrm{1}S}(\cos \mathit{\alpha }{)}^{\mathrm{2}}+R_{\mathrm{2}S}(\sin \mathit{\alpha }{)}^{\mathrm{2}}\\ {\stackrel{\mathrm{̃}}{R}}_{\mathrm{2}}& =R_{\mathrm{2}S}[\mathrm{1}-(\sin \mathit{\alpha }{)}^{\mathrm{2}}/\mathrm{2}]+R_{\mathrm{1}S}(\sin \mathit{\alpha }{)}^{\mathrm{2}}/\mathrm{2}.\end{array}\end{array}$$

Collecting the eigenvalues (Eq. 74) in the diagonal matrix $\stackrel{\mathrm{̃}}{\mathcal{D}}=\text{diag}({\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{\mathrm{0}},{\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{-},{\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{+})$, we have ${\mathfrak{B}}^{-\mathrm{1}}\approx {\mathfrak{U}}_{\mathrm{0}}{\stackrel{\mathrm{̃}}{\mathcal{D}}}^{-\mathrm{1}}{\mathfrak{U}}_{\mathrm{0}}^{\mathrm{H}}$. As an example, the element in the lower right corner of the inverse matrix is

$$\begin{array}{}\text{(76)}& {\mathfrak{B}}_{\mathrm{33}}^{-\mathrm{1}}=F_z\approx (\cos \mathit{\alpha }{)}^{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{\mathrm{0}}^{-\mathrm{1}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\sin \mathit{\alpha }{)}^{\mathrm{2}}\left({\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{-}^{-\mathrm{1}}+{\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{+}^{-\mathrm{1}}\right).\end{array}$$

(This approximation of F_z was used in Sezer (2023a) without proof.)

From Eq. (76), and using the first equality in Eq. (22), we immediately find

$$\begin{array}{}\text{(77)}& T_i\approx \mathrm{Re}\left\{{\mathrm{c}}^{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{\mathrm{0}}^{-\mathrm{1}}+{\mathrm{s}}^{\mathrm{2}}\left({\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{-}^{-\mathrm{1}}+{\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{+}^{-\mathrm{1}}\right)/\mathrm{2}\right\}.\end{array}$$

To obtain v₊, we need to subtract $T_i^{\mathrm{0}}$ from T_i (Eq. 24). Since ω₁=0 implies s=0 and c=1,

$$\begin{array}{}\text{(78)}& T_i^{\mathrm{0}}=\mathrm{Re}\left\{{\left(R_{\mathrm{1}S}+i{\mathit{\omega }}_I\right)}^{-\mathrm{1}}\right\},\end{array}$$

which is identical to the exact result in Eq. (25). Hence,

$$\begin{array}{}\text{(79)}& v_{+}/{\mathit{\delta }}^{\mathrm{2}}\approx \mathrm{Re}\left\{{\mathrm{s}}^{\mathrm{2}}\left({\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{-}^{-\mathrm{1}}+{\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{+}^{-\mathrm{1}}\right)/\mathrm{2}+{\mathrm{c}}^{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{\mathrm{0}}^{-\mathrm{1}}-{\left(R_{\mathrm{1}S}+i{\mathit{\omega }}_I\right)}^{-\mathrm{1}}\right\}.\end{array}$$

One can similarly obtain the rate constant v₋ as a linear combination of the reciprocals of the approximate eigenvalues ${\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{\mathrm{0}}$ and ${\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{\mp }$. The result is

$$\begin{array}{}\text{(80)}& T_x={\displaystyle \frac{{\mathit{\omega }}_{\mathrm{1}}}{{\mathit{\omega }}_{\mathrm{eff}}^{\mathrm{2}}}}\mathrm{Re}\mathit{\left\{}iX\mathit{\right\}},\end{array}$$

with

$$\begin{array}{}\text{(81)}& X=R_{\mathrm{2}S}{\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{\mathrm{0}}^{-\mathrm{1}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(R_{\mathrm{2}S}+i{\mathit{\omega }}_{\mathrm{eff}}\right){\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{-}^{-\mathrm{1}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(R_{\mathrm{2}S}-i{\mathit{\omega }}_{\mathrm{eff}}\right){\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{+}^{-\mathrm{1}}.\end{array}$$

Recall that $v_{-}/{\mathit{\delta }}^{\mathrm{2}}=\left({\mathit{\omega }}_{\mathrm{1}}{f}_x\right)T_x$ (Eq. 24).

The first eigenvalue in Eq. (74) does not depend on the offset Ω. The other two eigenvalues depend on the offset through ω_eff. Let us consider sufficiently large offsets such that $\left|\mathrm{\Omega }\right|\gg {\mathit{\omega }}_{\mathrm{1}}$, and so ${\mathit{\omega }}_{\mathrm{eff}}\approx \left|\mathrm{\Omega }\right|$. This condition is satisfied at the solid-effect offset positions $\mathrm{\Omega }\approx \pm {\mathit{\omega }}_I$ at W and J bands but may be entirely inapplicable to X band at large mw powers, as discussed in Sezer (2023a). When the condition holds, s≈0 and c≈1, and the eigenvalues in Eq. (74) become ${\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{\mathrm{0}}\approx R_{\mathrm{1}S}+i{\mathit{\omega }}_I$ and ${\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{\mp }\approx R_{\mathrm{2}S}+i({\mathit{\omega }}_I\mp |\mathrm{\Omega }\left|\right)$. Thus ${\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{\mp }^{-\mathrm{1}}$ correspond to complex-valued Lorentzians centered at $\mathrm{\Omega }=\pm {\mathit{\omega }}_I$ and with widths equal to the homogeneous EPR line width (without power broadening). These are the Lorentzians that we see as narrow lines at W and J bands in the second and third rows of Fig. 3 (orange and red solid lines).

In the case of motion, assuming mono-exponential correlation function for simplicity, each eigenvalue is replaced by ${\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_n+\mathrm{1}/\mathit{\tau }$. This amounts to increasing the widths of the solid-effect Lorentzians from R_2S to $R_{\mathrm{2}S}+\mathrm{1}/\mathit{\tau }$. The resulting motional broadening is the reason for the differences between the “solid” and “liquid” lines in the second and third rows of Fig. 3.

Figure 5 Comparison between the exact (solid colored lines) and perturbative (dashed black lines) calculations of the FFHS model. All parameters as in Fig. 4. The eigenvalue problem of the generalized Bloch matrix 𝔅 has a simple closed-form solution when T_1S=T_2S. The perturbative approximation uses these analytical eigenvectors and corrects the eigenvalues to first order in the difference $\mathrm{1}/T_{\mathrm{1}S}-\mathrm{1}/T_{\mathrm{2}S}$.

For a general motional model, we have the approximate $\stackrel{\mathrm{̃}}{\mathcal{L}}=\text{diag}\left(J_{\mathrm{11}}\right({\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{\mathrm{0}}),J_{\mathrm{11}}({\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{-}),J_{\mathrm{11}}({\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{+}\left)\right)$, which yields the approximation $\mathfrak{U}\mathcal{L}{\mathfrak{U}}^{-\mathrm{1}}\approx {\mathfrak{U}}_{\mathrm{0}}\stackrel{\mathrm{̃}}{\mathcal{L}}{\mathfrak{U}}_{\mathrm{0}}^{\mathrm{H}}$ to be used in Eq. (61). For the spectral density of the FFHS model, the perturbative expressions are compared with the exact numerical calculation in Fig. 5. The former are plotted with dashed–dotted black lines and the latter with colored solid lines, like in Fig. 4. We see that the perturbative analysis is satisfactory in general, at least for the specific choice of parameters that were used. It gives excellent predictions for v₋ (Fig. 5, fourth row) and, because the two are related by a global scaling factor (Eq. 10), also for the DNP enhancement (fifth row). At the same time, it is seen to consistently fail for the rate v₊ at small offsets in the vicinity of the origin (third row).

We should mention that the perturbative approximation becomes progressively better when R_2S approaches R_1S (not shown), as it is exact for R_2S=R_1S.

Leaving the approximation quality of the perturbative analysis aside, we observe that the enhancement profiles in the fifth row of Fig. 5 reveal the emergence of a novel feature at small offsets. At W band, this feature appears as a shoulder in the broadened lines, and at J band it is already separated from the canonical solid-effect peaks. Comparison with the lines in the first row of Fig. 5 makes it clear that this new feature in the DNP spectrum coincides with the extrema of the dispersive component of the power-broadened EPR line. For saturating mw powers, where ω₁≫R_1SR_2S, these extrema are at

$$\begin{array}{}\text{(82)}& {\mathrm{\Omega }}_{\mathrm{1}/\mathrm{2}}={\mathit{\omega }}_{\mathrm{1}}\sqrt{T_{\mathrm{1}S}/T_{\mathrm{2}S}}.\end{array}$$

(The subscript $\mathrm{1}/\mathrm{2}$ was selected because these are also the offset positions where the electronic saturation factor equals one half.) The factor pv₋ in the fourth row of Fig. 5 is obtained as the product of the first and second rows, as elaborated in Sezer (2023a). When the solid-effect lines at $\mathrm{\Omega }\approx \pm {\mathit{\omega }}_I$ (second row) become sufficiently broad, their amplitude at ${\mathrm{\Omega }}_{\mathrm{1}/\mathrm{2}}$ gets large enough for the peak of the dispersive EPR line (first row) to be visible in the DNP spectrum.

4.3 Motional suppression and broadening

Let us examine more closely the suppression of the lines at the solid-effect offsets and the concurrent increase of their intensity at ${\mathrm{\Omega }}_{\mathrm{1}/\mathrm{2}}$. We will limit the discussion to J band where the condition ω_I≫ω₁ holds, and the solid-effect offsets are $\mathrm{\Omega }=\pm {\mathit{\omega }}_I$. Using the perturbative eigenvalues in Eq. (74), we see that only the real parts of ${\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{\mp }$ survive at these offsets. The peak amplitudes of the solid-effect lines are then proportional to J₁₁(R_2S).

The limits τ≪T_2S and τ≫T_2S correspond to, respectively, very fast and very slow diffusive motion relative to the electronic T₂. In the slow limit τ≫T_2S, we have

$$\begin{array}{}\text{(83)}& \underset{R_{\mathrm{2}S}\mathit{\tau }\to \mathrm{\infty }}{lim}{J}_{\mathrm{11}}\left(R_{\mathrm{2}S}\right)/〈{\mathit{\delta }}^{\mathrm{2}}〉\to T_{\mathrm{2}S}\end{array}$$

for both the mono-exponential and FFHS motional models. This means that, in solids, the peaks increase with the electronic T₂. In the opposite limit of very fast motion, i.e., τ≪T_2S, we find

$$\begin{array}{}\text{(84)}& \underset{R_{\mathrm{2}S}\mathit{\tau }\to \mathrm{0}}{lim}{\displaystyle \frac{J_{\mathrm{11}}\left(R_{\mathrm{2}S}\right)}{〈{\mathit{\delta }}^{\mathrm{2}}〉}}\to \left\{\begin{array}{cc}\mathit{\tau }& \mathrm{exponential}\\ \frac{\mathrm{4}}{\mathrm{9}}\mathit{\tau }& \mathrm{FFHS}\end{array}\right.\end{array}$$

which means that, in liquids, the peak amplitudes are proportional to the dipolar correlation time τ. In other words, faster fluid diffusion (i.e., smaller τ) corresponds to smaller peaks and thus smaller solid-effect enhancement at the canonical offsets. We also see that for the same τ the peaks of the FFHS model are less than half of the peaks of the exponential model, which is in agreement with Fig. 4 (second and third rows).

To describe the transition between the fast and slow limits, we define the reduction factor

$$\begin{array}{}\text{(85)}& \mathit{\rho }\left(\mathit{\tau }\right)={\displaystyle \frac{\mathrm{1}}{T_{\mathrm{2}S}}}{\displaystyle \frac{J_{\mathrm{11}}(R_{\mathrm{2}S};\mathit{\tau })}{\langle {\mathit{\delta }}^{\mathrm{2}}\rangle }}\end{array}$$

which equals one in the solid limit and approaches zero for small τ. Since it quantifies how much smaller the peaks are compared to the solid case, ρ is a measure of how “solid-like” the liquid is.

Figure 6 Multiplicative deviation of T_x(Ω) from the solid limit at (a) the solid-effect offsets $\mathrm{\Omega }=\pm {\mathit{\omega }}_I$ and (b) the offsets ${\mathrm{\Omega }}_{\mathrm{1}/\mathrm{2}}=\pm {\mathit{\omega }}_{\mathrm{1}}\sqrt{T_{\mathrm{1}S}/T_{\mathrm{2}S}}$ for the exponential and FFHS models.

Figure 6a shows the reduction factors of the exponential and FFHS models against the relative motional timescale $\mathit{\tau }/T_{\mathrm{2}S}$. For the exponential model, the solid-effect peaks drop to half of their maximum values at τ=T_2S. In the case of the FFHS model, this happens already at τ≈4T_2S (see inset). In other words, appreciable reduction compared to the solid limit occurs even for exceedingly long diffusive timescales, several-fold compared to the electronic T₂. For identical τ's the exponential model is seen to be more solid-like than the diffusive FFHS model across the entire motional range. Hence, realistic translational diffusion suppresses the solid-effect peaks more effectively than mono-exponential decay.

As a quantitative measure of the motional broadening, let us consider the magnitude of T_x (Fig. 5, second row) at the locations of the extrema of the dispersive EPR line (first row). Since the intensity at these small offsets increases when going from the solid to the liquid case, we define the magnification factor

$$\begin{array}{}\text{(86)}& \mathit{\mu }\left(\mathit{\tau }\right)={\displaystyle \frac{T_x^{\mathrm{liquid}}({\mathrm{\Omega }}_{\mathrm{1}/\mathrm{2}};\mathit{\tau })}{T_x^{\mathrm{solid}}\left({\mathrm{\Omega }}_{\mathrm{1}/\mathrm{2}}\right)}}.\end{array}$$

This factor is shown in Fig. 6b. For the FFHS model, the intensity at ${\mathrm{\Omega }}_{\mathrm{1}/\mathrm{2}}$ is 2 to 3 orders of magnitude larger than the solid case across a broad range of motional timescales, between τ=T_2S and τ=0.01T_2S. Hence the peak of the dispersive EPR line should be magnified 100- to 1000-fold in liquids compared to the solid limit. It is also magnified for the mono-exponential dipolar correlation function, although not to the same extent.

In the light of these observations, next we analyze the DNP field profile of recent experiments with the free radical BDPA in DMPC lipid bilayers at 320 K (Kuzhelev et al., 2022).

4.4 Comparison with experiment

The DNP experiments of Kuzhelev et al. (2022) were carried out at J band (260 GHz/400 MHz). For the acyl chain protons of the DMPC lipids, the peak DNP enhancements at the canonical solid-effect offsets were ±12 (Kuzhelev et al., 2022). Two additional enhancement peaks of ±8 were also observed at much smaller offsets. These were attributed to thermal mixing. Here we argue that they correspond to the extrema of the dispersive component of the EPR line.

The enhancements in Kuzhelev et al. (2022) were for a BDPA-to-lipid ratio of 1:10 at a temperature of about 320 K. The room-temperature EPR spectrum of BDPA at J band for this relatively high radical concentration was very narrow (Kuzhelev et al., 2022, Fig. 2). The transverse relaxation time implied by this narrow line is T_2S=215 ns. For the same radical concentration, the nuclear spin-lattice relaxation time at J band was 50 ms at 298 K (Kuzhelev et al., 2022). Although the experimental T_1I and T_2S are for 298 K, below we use these values to fit the DNP spectrum at 320 K. We also use B₁=6 G, as estimated in Kuzhelev et al. (2022).

In addition to these parameters with experimental support, three more parameters are needed for the calculation of the DNP enhancement: T_1S, τ and $N/\mathrm{3}{b}^{\mathrm{3}}$. We will treat these as fitting parameters. Let us introduce the ratios

$$\begin{array}{}\text{(87)}& r_{\mathrm{1}}={\displaystyle \frac{T_{\mathrm{1}S}}{T_{\mathrm{2}S}}},r_{\mathrm{2}}={\displaystyle \frac{T_{\mathrm{2}S}}{\mathit{\tau }}},r_{\mathrm{3}}={\displaystyle \frac{N/\mathrm{3}{b}^{\mathrm{3}}}{N_{\mathrm{ref}}/\mathrm{3}{b}_{\mathrm{ref}}^{\mathrm{3}}}}.\end{array}$$

The ratio r₁ expresses the unknown electronic T₁ relaxation time in terms of the known electronic T₂. From physical considerations, r₁≥1. The ratio r₂ relates the diffusion timescale τ to the electronic T₂. Since T_2S is rather long, we expect τ to be shorter and hence r₂≥1. Finally, the ratio r₃ expresses the actual factor $N/\mathrm{3}{b}^{\mathrm{3}}$, which is unknown, as a multiple of this same factor for arbitrarily selected reference values N_ref and b_ref. In principle, r₃ can be any positive number.

The mean volume per particle at a concentration of 1 M is 1.66 nm³ and corresponds to a cube with side length of 1.18 nm. From molecular modeling, the “radius” of a BDPA molecule is about 0.6 nm, so it barely fits in the above cube. The partial molecular volume of a DMPC lipid in a lipid bilayer is 1.1 nm³ (Greenwood et al., 2006). Thus, the concentration of one BDPA when surrounded by 10 DMPC lipids cannot exceed N_ref=0.1 M but is also likely close to this value. Additionally taking b_ref=1 nm in the last equality of Eq. (87), we anticipate r₃>1.

From Eq. (10), the expected dependence of the DNP enhancement on the three fitting parameters can be written as

$$\begin{array}{}\text{(88)}& {\mathit{ϵ}}_{\mathrm{SE}}\approx p\left(r_{\mathrm{1}}\right)v_{-}(r_{\mathrm{1}},r_{\mathrm{2}})r_{\mathrm{3}}\langle {\mathit{\delta }}_{\mathrm{ref}}^{\mathrm{2}}\rangle T_{\mathrm{1}I}\left|{\mathit{\gamma }}_S\right|/{\mathit{\gamma }}_I,\end{array}$$

where $\langle {\mathit{\delta }}_{\mathrm{ref}}^{\mathrm{2}}\rangle$ is calculated according to Eq. (27) using N_ref and b_ref. The ratio r₁, which determines ${\mathrm{\Omega }}_{\mathrm{1}/\mathrm{2}}$ (Eq. 82), will influence the electronic polarization factor (or, equivalently, saturation factor). Together with the ratio r₂, it will also influence the forbidden transition rate v₋, although the effect of r₁ is expected to be small. From the previous discussion, we expect that r₂ will mostly be responsible for the width of the solid-effect lines that comprise v₋. Finally, r₃ will serve as a global scaling factor that will adjust the amplitude of the overall enhancement. Since the three fitting parameters are responsible for different features of the DNP spectrum, it should be possible to determine them uniquely.

The top plot in Fig. 7 shows the experimental enhancements (red circles) together with the best fit obtained using the FFHS model with B₁=6 G (solid black line). The solid green line in this plot is the difference between the experimental data and the fit. The corresponding fitting parameters are shown in the first row of Table 1. Note that the fits were performed numerically using the exact expressions of p, v_±, and the DNP enhancement (Eq. 9) and were not restricted to the dependencies on the fitting parameters r_i ($i=\mathrm{1},\mathrm{2},\mathrm{3}$) that are indicated in Eq. (88). (For example, the general dependence of the electronic polarization factor p on T_1S and T_2S is not limited to the ratio $T_{\mathrm{1}S}/T_{\mathrm{2}S}$.)

After the fits converged, we used the final values of the fitting parameters to calculate the dispersive component of the EPR line $s_x^{\mathrm{ss}}=p{\mathit{\omega }}_{\mathrm{1}}{f}_x{s}_z^{\mathrm{eq}}$, and the factor T_x, such that $\left(p{\mathit{\omega }}_{\mathrm{1}}{f}_x\right)T_x=p{v}_{-}$. These are shown in the lower plot of Fig. 7, where $s_x^{\mathrm{ss}}$ (blue) and T_x (orange) are scaled independently along the vertical axis. Their product (dotted–dashed black line) is also scaled independently along the y axis. Since p_X≈1 in our case, the product pv₋ is itself proportional to the solid-effect enhancement. Hence the black lines in the upper and lower plots of Fig. 7 are directly comparable. We can thus visually conclude that the unusual enhancement peaks at small offsets are a direct manifestation of the dispersive component of the power-broadened EPR line.

Figure 7 Experimental DNP field profile at J band (red circles) and fit with the FFHS model using B₁=6 G (solid black line). The central peaks in the enhancement profile follow the dispersive component of the power-broadened EPR line (solid blue line).

Table 1 Fitting parameters r_i ($i=\mathrm{1},\mathrm{2},\mathrm{3}$) and implied timescales (T_1S and τ) and contact distance (b). For different mw fields, r₁ changes such that $B_{\mathrm{1}}\sqrt{r_{\mathrm{1}}}$ remains constant.

^* Assuming radical concentration N=0.1 M.

Because the mw field in the experiment is not known precisely, we also attempted fits with smaller B₁. The best fits were practically identical to the one shown in Fig. 7 but with different values of the fitting parameters (Table 1). (The fit for B₁=2 G is shown in Fig. A1 as an example.)

From Table 1 we see that all fits resulted in the same value of the parameter r₂, implying τ=7.1 ns for the motional timescale of the FFHS model. This parameter is very robust because it directly reflects the width of the experimental solid-effect lines at $\mathrm{\Omega }=\pm {\mathit{\omega }}_I$. Although, normally, their line width should depend on both T_2S and τ, the exceptionally narrow EPR line puts us in the regime τ≪T_2S where the influence of T_2S is negligible. As a result, the motional broadening of the solid-effect lines in the DNP spectrum reports directly on the diffusive timescale of BDPA in the lipid environment.

The fitting parameter r₁ adjusts the extrema of the dispersive EPR line (solid blue line in Fig. 7), which are at ${\mathrm{\Omega }}_{\mathrm{1}/\mathrm{2}}=\pm {\mathit{\omega }}_{\mathrm{1}}\sqrt{r_{\mathrm{1}}}$ (Eq. 82). By monitoring the product of $\sqrt{r_{\mathrm{1}}}$ and B₁ in the last column of Table 1, we see that, each time B₁ is modified, the fitted r₁ changes such that ${\mathrm{\Omega }}_{\mathrm{1}/\mathrm{2}}$ remains unchanged, as required by the positions of the non-canonical enhancement peaks in the experimental data. However, because r₁ has to increase quadratically to compensate for the reduction of B₁, the implied electronic T₁ times become exceedingly long (tens of microseconds) at the smaller values of B₁ (1–2 G).

When B₁ is reduced, the fitting parameter r₃ also increases quadratically to compensate for the dependence of the overall enhancement on ${\mathit{\omega }}_{\mathrm{1}}^{\mathrm{2}}$. Assuming N_ref=0.1 M is a good estimate of the actual concentration of BDPA in the lipid bilayer, it is possible to calculate a contact distance, b, from the fitted value of r₃. The deduced contact distances are given in the second last column of Table 1. Only the values for large B₁ (5–6 G) are in qualitative agreement with the molecular structure of BDPA.

Figure 8 Same as Fig. 7 for the model with mono-exponential dipolar correlation function and B₁=6 G (solid black line). The fit parameters r₁=3.6, r₂=89 and r₃=7.7 correspond to T_1S=0.77 µs, τ=2.4 ns and b=0.507 nm (assuming N=0.1 M).

In Fig. 8 we show a fit to the same experimental data using a mono-exponential dipolar correlation function. While the difference between the data and the fit (green line in top panel of Fig. 8) is not much worse than what we had for the FFHS model, it is apparent that the exponential model strives to find the right balance between the broadening of the solid-effect lines and the tails of these lines at the lower offsets, ultimately producing too broad solid-effect lines and too narrow non-canonical peaks. (An exponential fit with B₁=4 G is shown in Fig. A2.) We thus see that the J-band DNP spectrum clearly differentiates between two alternative motional models, ruling out the less realistic one.

The most certain outcome of the fits with the FFHS model is the deduced motional timescale τ, as it comes directly from the width of the solid-effect lines (T_2S is too long to contribute). The deduced value of T_1S is somewhat less certain since it is accessed relative to T_2S and also depends on the mw field B₁. Nevertheless, with reasonable choices of T_2S and B₁, the fit to the non-canonical extrema in the DNP spectrum restricts T_1S to a meaningful window between 1.5 and 2.5 µs (Table 1). Least certain is the estimate of the contact distance b since, in addition to B₁, it requires precise knowledge of the radical concentration and the nuclear spin-lattice relaxation time. Although the latter is accessible experimentally, its value was measured at 298 K, while the DNP measurements are at 320 K.

In spite of the uncertainty in the estimated value of the distance parameter b, let us use τ and b in Eq. (67) to calculate the coefficient of relative translational diffusion. With b=0.482 nm (Table 1, first row) and τ=7.1 ns, we get $D=b^{\mathrm{2}}/\mathit{\tau }=\mathrm{0.033}{\mathrm{nm}}^{\mathrm{2}}{\mathrm{ns}}^{-\mathrm{1}}=\mathrm{33}\times {\mathrm{10}}^{-\mathrm{12}}{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$. Alternatively, with b=0.427 nm (Table 1, second row) and τ=7.1 ns, we find $D=\mathrm{26}\times {\mathrm{10}}^{-\mathrm{12}}{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$. The first value corresponds to B₁=6 G and the second to B₁=5 G.

For comparison, the coefficient of lateral diffusion of phospholipids in oriented DMPC bilayers, as determined from pulsed field gradient NMR, is about $\mathrm{11}\times {\mathrm{10}}^{-\mathrm{12}}{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ at 308 K, $\mathrm{20}\times {\mathrm{10}}^{-\mathrm{12}}{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ at 323 K and $\mathrm{27}\times {\mathrm{10}}^{-\mathrm{12}}{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ at 333 K (Filippov et al., 2003, Fig. 5b, 0 mol % cholesterol). As the temperature of the DNP measurements is closer to the middle value, our two estimates of D are seen to be larger by a factor of 1.65 and 1.3, respectively.

However, the D of the FFHS model corresponds to the relative translational diffusion of the electronic and nuclear spins, i.e., $D=D_S+D_I$, where D_S and D_I denote the coefficients of translational diffusion of the two spin types. Disregarding all complicating factors, one could thus take $D_I=\mathrm{20}\times {\mathrm{10}}^{-\mathrm{12}}{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ from the literature value and rationalize the values that we deduced from the width of the solid-effect DNP lines as implying either $D_S=\mathrm{13}\times {\mathrm{10}}^{-\mathrm{12}}{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ or $D_S=\mathrm{6}\times {\mathrm{10}}^{-\mathrm{12}}{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ for the diffusion coefficient of the free radical BDPA in the lipid bilayer. As the obtained numerical values are rather plausible, we conclude that the quantitative analysis of the J-band DNP spectrum leads to meaningful molecular properties. Without the theoretical framework developed in this paper, neither the molecular distance b nor the diffusion coefficient D would be accessible from a solid-effect DNP spectrum in the liquid state.

5 Concluding discussion

Erb, Motchane and Uebersfeld had the hunch that the dispersive component of the EPR line is reflected in the solid-effect DNP enhancement (Erb et al., 1958a). A theoretical justification of their intuition was provided in the companion paper (Sezer, 2023a). Here, the formalism was extended to the solid effect in liquids. Our theoretical predictions were compared with recent DNP measurements at high field (Kuzhelev et al., 2022). The comparison demonstrated that, under appropriate conditions, the dispersive component of the EPR line is literally visible in the field profile of the DNP enhancement. Provided that seeing is believing, we have thus closed the circle.

The DNP mechanism which became known as the solid-state effect due to Abragam (Abragam and Proctor, 1958) had been observed in liquids from the very beginning (Erb et al., 1958a, b). Nevertheless, perhaps because it yielded comparatively smaller absolute enhancements and often coexisted with the Overhauser effect (Leblond et al., 1971b), the solid effect has remained less explored in liquids compared to solids. The recent use of this DNP mechanism as a new modality for probing the molecular dynamics in ionic liquids (Neudert et al., 2017; Gizatullin et al., 2021b) and its first applications at high magnetic field (Kuzhelev et al., 2022) indicate that the potential of the solid effect in the liquid state is yet to be harvested. A theoretical understanding of the mechanism in liquids is clearly going to be helpful in these endeavors. Developing the needed theory has been the main aim of the companion and current papers.

Admittedly, a theoretical description of the solid effect in liquids was developed more than 50 years ago by Korringa and colleagues (Papon et al., 1968; Leblond et al., 1971a). In fact, their analysis was much more ambitious than ours, as it aimed to quantify the DNP spectrum during the transition from the Overhauser effect to the solid effect upon reduction of the experimental temperature (Leblond et al., 1971b). Thus, in addition to the secular terms of the dipolar interaction that we considered here and in Sezer (2023a), their Hamiltonian also contained the non-secular terms, which are important for the cross-relaxation rates of the Overhauser effect, as well as the orientation-dependent part of the electronic Zeeman interaction, which determines the electronic relaxation rates and thus the degree of saturation. Following the prescription of second-order time-dependent perturbation theory, Korringa et al. derived equations for the deviations of both the electronic and nuclear polarizations from their values at thermal equilibrium (Papon et al., 1968).

The analytical framework of Korringa and colleagues had two additional aspects. First, as is well known, the semi-classical description of spin-lattice relaxation relaxes the system to infinite temperature. The usual way of correcting for this shortcoming in magnetic resonance is to subtract the correct thermal equilibrium from the right-hand side of the dynamical equation of the density matrix (Abragam, 1961). Instead, Korringa (1964) imposed the correct temperature by writing the equation of motion of the density matrix for complex-valued time, whose imaginary part was proportional to the inverse temperature. This mathematical trick exploits the fact that a quantum-mechanical propagator with imaginary time becomes a Boltzmann factor. The analytical continuation to complex time modified the familiar Liouville–von Neumann equation of the density matrix to a form that is not common in magnetic resonance. Second, as an integral part of their formalism, Korringa et al. (1964) modeled the stochastic modulation of the spin Hamiltonian as rotational diffusion of one coordinate frame with respect to another, which led to an exponential correlation function with single decay time τ. It is not straightforward to see how their final analytical expressions should be modified if one were to use the FFHS model, for example.

In this context, it is worth mentioning that the mono-exponential model did not accurately fit the experimental data of Leblond et al. (1971b), and the authors took a Gaussian distribution for ln τ (Leblond et al., 1971a). In Fig. 8 we also observed that an exponential motional model did not fit the experimental DNP spectrum at J band (Kuzhelev et al., 2022), whereas the FFHS model with a single motional parameter did (Fig. 7). One should remember, however, that the analysis of Leblond et al. (1971b) was performed 4 years before the spectral density of the FFHS model was solved analytically (Ayant et al., 1975; Hwang and Freed, 1975).

In spite of the differences between the analytical framework of Korringa and colleagues (Papon et al., 1968; Leblond et al., 1971a) and our approach, which hamper a direct comparison of the results, we observe that the derivations in their first paper (Papon et al., 1968) assumed isotropic electronic relaxation, i.e., T_1S=T_2S. As we saw in Sect. 4.2, in this case the matrix 𝔅 (Eq. 13) becomes equal to 𝔅₀ (Eq. 69) and the eigenvalue problem of the latter has a simple closed-form solution. All quantities of interest then become linear combinations of the reciprocals of the eigenvalues. Indeed, the final expressions of Papon et al. (1968) are linear combinations of Lorentzian spectral densities, which contain the effective frequency ω_eff (Eq. 71).

The assumption of equal longitudinal and transverse electronic relaxation rates was relaxed in the second paper (Leblond et al., 1971a). Sadly, this second paper has been cited only five times, and although all of the citing papers report new experiments, they do not use the theoretical expressions of Leblond et al. (1971a) to analyze the experimental data. One can only hope that, by being less ambitious, the theory developed in the current paper fares differently.

Appendix A: Additional figures

Figure A1 Same as Fig. 7 (FFHS model) for B₁=2 G (solid black line). The fit parameters are given in the second last row of Table 1.

Figure A2 Same as Fig. 8 (exponential model) for B₁=4 G (solid black line). The fit parameters r₁=8.1, r₂=90, r₃=17 correspond to T_1S=1.7 µs, τ=2.4 ns, b=0.387 nm (assuming N=0.1 M).