Dynamic nuclear polarization as kinetically constrained diffusion

Dynamic nuclear polarization (DNP) is a promising strategy for generating a significantly increased non-thermal spin polarization in nuclear magnetic resonance (NMR) applications thereby circumventing the need for strong magnetic fields. Although much explored in recent experiments, a detailed theoretical understanding of the precise mechanism behind DNP is so far lacking. We address this issue by theoretically investigating solid effect DNP in a system where a single electron is coupled to an ensemble of interacting nuclei and which can be microscopically modelled by a quantum master equation. By deriving effective equations of motion that govern the polarization dynamics we show analytically that DNP can be understood as kinetically constrained spin diffusion. On the one hand this approach provides analytical insights into the mechanism and timescales underlying DNP. On the other hand it permits the numerical study of large ensembles which are typically intractable from the perspective of a quantum master equation. This paves the way for a detailed exploration of DNP dynamics which might form the basis for future NMR applications.

Dynamic nuclear polarization (DNP) is a promising strategy for generating a significantly increased non-thermal spin polarization in nuclear magnetic resonance (NMR) applications thereby circumventing the need for strong magnetic fields. Although much explored in recent experiments, a detailed theoretical understanding of the precise mechanism behind DNP is so far lacking. We address this issue by theoretically investigating solid effect DNP in a system where a single electron is coupled to an ensemble of interacting nuclei and which can be microscopically modelled by a quantum master equation. By deriving effective equations of motion that govern the polarization dynamics we show analytically that DNP can be understood as kinetically constrained spin diffusion. On the one hand this approach provides analytical insights into the mechanism and timescales underlying DNP. On the other hand it permits the numerical study of large ensembles which are typically intractable from the perspective of a quantum master equation. This paves the way for a detailed exploration of DNP dynamics which might form the basis for future NMR applications.

PACS numbers: Valid PACS appear here
The sensitivity of imaging and spectroscopy techniques based on nuclear magnetic resonance (NMR) depends on the spin polarization that arises from the Zeeman interaction of the nuclei with an externally applied magnetic field. However, even with the highest magnetic fields that can currently be generated by superconductive magnets, the nuclear spin polarization of 1 H nuclei is only in the order of less than ten parts per million at ambient temperature. A promising strategy to overcome the low sensitivity in many NMR applications is the use of dynamic nuclear polarization (DNP). Here the electronic spin polarization -that is about three orders of magnitude higher than the nuclear spin polarization -is transferred from paramagnetic centres in solid state materials to the surrounding nuclear spin ensemble by irradiating with microwave fields at an appropriate frequency. DNP was already proposed and demonstrated in the early days of nuclear magnetic resonance [1][2][3], but it has recently attracted increased attention since progress in hardware technology, mainly in the form of stable high frequency and high power microwave sources, has made it possible to substantially enhance the NMR signal in several applications [4,5] including the detection of nuclei with small gyromagnetic constants [6], of nuclei located on the surface of structured materials [7] and of nuclei in protein micro-crystallites [8] by magic angle spinning solid state NMR spectroscopy. Furthermore, DNP carried out at cryogenic temperatures combined with a fast dissolution step to rapidly bring the sample to ambient temperature, was successfully used to generate highly polarised 13 C-labelled compounds that could be detected in in vivo Magnetic Resonance Imaging applications after their injection [9,10]. A full understanding of the spin dynamics during DNP is essential for the interpretation of the experimental data and ultimately for optimising the experimental protocols. In principle this problem can be addressed by numerically solving the full quantum master equation [11,20]. Owing to its complexity, however, this approach is limited to small systems [12,13].
In this work we show how to overcome this limitation and to extend the analysis of the spin dynamics to large ensembles containing more than a thousand coupled spins. By using an adiabatic elimination procedure we derive an effective equation for the evolution of the spin polarization. This equation is classical in the sense that it describes transitions between Zeeman basis states and therefore the computational complexity of simulating the spin dynamics is drastically reduced. Moreover, the analytical form of the effective equations of motion permits a transparent interpretation of the physics underlying DNP in terms of operator-valued, i.e kinetically constrained, diffusion rates.
Solid effect DNP (SE DNP) is observed in systems that consist of interacting unpaired electrons S and nuclei I [see Fig. 1(a)] in an external magnetic field with Lamor frequencies ω S and ω I . Polarization buildup is mediated by the electrons that are irradiated at an off-resonance frequency ω 0 which is in resonance with either the socalled double or zero quantum transition ω S ± ω I [14][15][16][17][18][19][20]. For SE DNP to be effective, the electrons need to have a resonance linewidth that is smaller than the nuclear Larmor frequency. A suitable minimal model for SE DNP consists of an ensemble of n nuclei coupled to a single unpaired electron through the hyperfine interaction. The spin dynamics of SE DNP can be described by the quantum master equation dρ/dt = Lρ with the Liouvillian L = −iĤ + Γ. It describes the coherent evolution of the spin system under the action of the Hamiltonian commutation superoperatorĤ and the effect of the dissipative processes (decoherence) represented by the relaxation superoperator Γ. Note, that throughout a bold symbol is  1. (a) Sketch of the solid effect DNP setup. The electron is coherently driven by the applied microwave field and coherent excitation-exchange between the spins is mediated through dipole-dipole couplings. Furthermore the electron and nuclei are subject to dissipation leading to relaxation (see text for details). (b) The effective dynamics of the polarization is described by a classical master equation and thus permits the simulations of large systems. (c) Comparison between the exact (solid black lines) and effective (colored lines) relaxation dynamics of a system consisting of one electron and four nuclei. (d) Polarization dynamics in a linear spin chain with one electron and 30 nuclear spins. The polarization versus spin index (0 denotes the electron position) is shown for different times and the polarization versus time is shown for the electron (blue line) and four selected nuclear spins (red lines). used to represent superoperators and a doublehat for a commutation superoperator, i.e.Ô ≡ [Ô, ·].
In the following we discuss the structure of the quantum master equation in detail. The HamiltonianĤ of the system consists of four terms:Ĥ =Ĥ Z +Ĥ 0 +Ĥ + +Ĥ − , whereĤ Z = ω I (Ŝ z + kÎ kz ) is the Zeeman term. The term describes the hyperfine and nuclear dipole interaction. It depends on the secular strengths of the electron-nuclear hyperfine interaction, parameterized by A k and furthermore on the inter-nuclear dipolar coupling coefficients d jk [21]. There is an additional dependence on the so-called offset parameter λ = ω S − ω I − ω 0 which represents a resonance condition between the microwave frequency ω 0 and that of the double quantum transition ω S − ω I . The remaining terms of the system Hamiltonian arê which arise from the pseudosecular hyperfine interaction and the microwave irradiation. They are parameterized by the pseudosecular coupling strength B k± between the electron and the nuclei and the microwave amplitude ω 1 . Note, that the interaction parameters A k , B k± and d kj depend on the geometry of the spin system [21]. The relaxation superoperator is composed of two terms, Γ = Γ 1 + Γ 2 . The first one is given by where we have defined the dissipator D(X)ρ ≡XρX † − {ρ,X †X }/2. Γ 1 describes dissipative processes or longitudinal relaxation due to interaction of the spins with the environment [22]. The rates for longitudinal relaxation are R ]/2 acts as a correction term and ensures that the electron relaxes back to its thermal equilibrium state which-at typical temperatures of T > 0.1K-is up to a scaling factor approximated by ρ th = 1 − 2P 0Ŝz . Here P 0 = tanh ( ω S /2k B T ) is the thermal electronic polarization [21]. The second term, represents the loss of coherences or transverse relaxation of the spin ensemble, which depends on the rates R (S) 2 and R (I) 2 for the electron and the nuclear spins, respectively. The full quantum master equation describes the dynamics of the spin system in the Liouville space L of dimension 2 2N where N = n + 1.
However, all information that is relevant for characterizing SE DNP is typically contained in the Zeeman subspace L Z [23] which is spanned by the operators {1,Î kz ,Î kzÎk ′ z ,Î kzÎk ′ zÎk ′′ z , . . .} and has a dimension 2 N . For example the bulk nuclear polarization is given by the expectation value of nÎ kz . We have obtained a closed equation for the dynamics in the Zeeman subspace by means of a quite lengthy adiabatic approximation procedure which is detailed in the supplemental material. This resulting effective master equation is of the form dρ Z /dt = L ZρZ . It has not only a substantially reduced dimension but in addition it is of entirely classical nature -although we employ in the following a notation that is also used in the quantum case. As shown in the supplemental material the dynamics can be expressed in terms of jump operators [see Fig. 1(b)], which connect classical spin configurations, and their corresponding effective jump rates, i.e. the equation is of Lindblad form [24]. Similar approaches have been employed in the study of cold atomic systems [25][26][27][28]. When formulated in this way the problem lends itself to the use of classical kinetic Monte Carlo (kMC) algorithms which permits the simulation of large systems. To check the validity of this approximate approach we have simulated small systems with the exact quantum Master equation and kMC. Example results are shown in Fig. 1(c) as well as in the supplemental material and show excellent agreement between the two methods.
The effective dynamics in the Zeeman subspace is governed by a set of single-spin and two-spin dissipators: L Z = L single−spin + L two−spin . The single-spin terms are given by and depend on the constant rates Here Γ (S) ± are the rates at which the electron spin flips 'up' (+) or 'down' (−). They depend on the longitudinal relaxation rate R (S) 1 of the electron with a weighting pre-factor which ensures that in the absence of any perturbation by the microwave field the electronic steadystates polarization is identical to the thermal electronic polarization P 0 . The second term is a consequence of the applied microwave field. Eq. (3) describe the rates for a nuclear spin k to flip 'up' or 'down'. Since we have assumed in the derivation of these equations that the thermal nuclear polarization is negligibly small, these two rates coincide and correspond to half the longitudinal nuclear relaxation rate R (I) 1 . The additional term appearing in Eq. (3) arises from the pseudosecular interaction with the electron which leads to a tilt of the nuclear effective field axis.
The two-spin dissipators are given as They clearly describe a diffusive (flip-flop) dynamics between the electron and the nuclei as well as among the nuclei. Interestingly the diffusion rates Γ are not constant, but operator-valued. Hence the diffusion rate of the spin-polarization depends on the state of the entire spin system. In the context of the study of glassy relaxation [29,30] such operator-valued rates are often referred to as kinetic constraint. Explicitly, they . (6) A brief discussion of these expressions provides insight into the spin dynamics on which SE DNP is based on. Eq. (5) describes the rate of flip-flop transitions or jumps between a nuclear spin k and the unpaired electron. Hencê Γ (IS) k is of fundamental importance for the SE DNP effect. This rate depends on both the strength of the microwave field ω 1 and the pseudosecular interaction strength B k± . It is inversely proportional to the square of the nuclear Larmor frequency ω I and also inversely proportional to the sum of the transverse relaxation rates of the electron, R (S) 2 , and the nuclear spin, R (I) 2 . This indicates that although this process is represented by the flip-flop jump operator, the underlying quantum mechanical process is mediated by coherences which in fact decay with these relaxation rates.
Most importantly, the rate for polarization transfer between the electron and the k-th nucleus is controlled by the operatorD k whose (eigen)value depend on the polarization level of the nuclear spin ensemble. For negligible nuclear polarization in an ensemble of nuclear spins the sum s =k A sÎsz evaluates to a value close to zero as there is an equal probability for the nuclei to be in their spin 'up' or 'down' state. However, with increasing level of nuclear polarization the magnitude of this sum grows which decreases the effective rateΓ (IS) k . Note, that the effect of the off-resonance parameter λ is negligible as long as the microwave frequency ω 0 is matched to the double quantum transition at ω S − ω I .
The transport of polarization within the dipolarcoupled network of the nuclear spins takes places via spin diffusion which is governed by the operator-valued rate given in Eq. (6). The rate of polarization transport among two nuclei is proportional to the square of their nuclear dipolar coupling constant d kj , and it is indirectly proportional to the nuclear transverse relaxation rate R (I) 2 . Interestingly, there is also a dependence on the difference of the secular hyperfine interaction strength A k of the two nuclei. In the case in which this difference is significant, the internuclear polarization transfer between the spins can become strongly quenched provided the electronic polarizationŜ z is non-zero.
As the effective dynamics in the Zeeman subspace can be treated efficiently numerically SE DNP can be studied in relatively large samples. In Fig. 1(d) we already showed a simple example for 31 spins in one dimension FIG. 2. Simulation of the polarization dynamics of 1330 nuclear spins ( 13 C) arranged on a regular grid in a cube (11 × 11 × 11) around one central electron (black). The diameter and the colour of the spheres indicate the expectation value of the nuclear spin polarization at the given position. The simulation demonstrates that polarization buildup is fastest for the nuclei that are located close to the electron at an angle of π/4 with respect to direction of the applied static magnetic field. In order to show the contribution of spin diffusion to the distribution of polarization within the spin ensemble the dipolar coupling constants were set to zero in the rightmost panel. This prevents any transport of polarization through spin flip-flop transitions between pairs of nuclear spins. Note, that the front quarter of the 13 C spin ensemble is not shown for better visibility of the polarization buildup near the electron.
while Fig. 2 now shows the progressive buildup of spin polarization within a three-dimensional sample of 1330 13 C spins arranged on a regular cubic lattice with one electron spin at the centre. For the latter the distance between the nuclei is 10 ± 5%Å and the field strength of the static external magnetic B 0 field is 3.4 T. The applied microwave field has a frequency of ω 0 = ω S − ω I = 34.5 MHz and a strength ω 1 =100 kHz. The relaxation time constants are R Let us finally investigate the role the nuclear dipolar interaction plays in the distribution of polarization between the nuclear spins. This question has been addressed e.g. by Hovav et al. in Ref. [11] who conducted simulations of small chains (up to 9 nuclei) within the framework of the full quantum master equation. Their results show that there is only a weak dependence of the polarization build up rates on the nuclear dipolar couplings between nuclei further away from the electron (nuclei belonging to the bulk) suggesting that there is a direct transfer of polarization from the electron to all bulk nuclei through mixing of the multi-nuclear product states. Our data shown in Fig. 2 clearly suggests that the nuclear dipolar couplings are important. To further investigate their role in polarization transport we have simulated the polarization dynamics of a linear spin chain consisting of 30 spins as shown in Fig. 1(d) [see supplemental material for further details]. In contrast to the analysis in Ref. [11] our findings show that a decrease of the nuclear dipolar interaction constants of the bulk nuclei by a factor of two alters the polarization buildup dramatically. This is consistent with a diffusive polarization transport. The discrepancy is thus most likely caused by boundary effects, which can be severe for short chains. Moreover, further discrepancies might arise due to the way in which relaxation was introduced in the reference frame in Ref. [11] as was already previously discussed in Ref. [22].
Let us now establish a more quantitative connection between SE DNP and the diffusion of polarization. For a linear chain of n nuclear spins, as shown in Fig. 1, we can derive an average diffusion constant: Here a k is the distance between the nuclear spins k and k + 1, d k is the dipolar coupling constant between the two spins and A k is the strength of the secular hyperfine interaction of spin k. The constant D k is the rate of diffusive polarization transport among neighboring nuclei. It is quenched when the difference between the secular hyperfine interactions of the two spins, |A k − A k+1 |, is large and when the nuclear transverse relaxation rate R (I) 2 is fast. Comparing the polarization dynamics calculated via kMC with a simple diffusion model, ∂p/∂t = D av where p ≡ p(x, t) is the polarization field, indeed shows good agreement with data as the one presented in Fig.  1(d). Further details and a discussion of the relevant boundary conditions is provided in the supplemental material. We conclude that polarization transport in a SE DNP experiment is mediated by nuclear spin diffusion. The presence of the electron will hinder -or constrain -this process for nuclei close to the electron since for these nuclei the difference between the strengths of the secular hyperfine interactions is significant [11,32].
We have derived effective equations of motions which show that the mechanism underlying SE DNP is kinetically constrained diffusion. The derived set of equations enable the study polarization buildup in relatively large spin systems which has been demonstrated by investigating DNP in a system with 1330 13 C nuclei and a single electron. Finally, we have highlighted the importance of nuclear dipolar couplings in SE DNP and derived an effective rate for polarization diffusion in a one-dimensional chain. We believe that this study represents an important step in the quantitative understanding of SE DNP for practical applications in NMR. In the future we will investigate whether the method presented here is also applicable in more complex settings, which involve e.g. more than a single electron.
We are grateful to Prof Shimon Vega, Weizman Institute, Rehovot, for interesting discussions related to DNP. and introducing the projectionsσ k = π kσ , L kj = π k Lπ j , k, j = 1, 2, we can writeσ Givenσ 1 , the second equation is resolved aŝ

Substituting this into the first equation, we obtain the equation closed in H
where K(T ) = L 12 e L22T L 21 .
In terms of the Laplace transforms (defined for ζ > 0) is resolved asl where The initial valueσ 1 (0) and the Laplace transform of the kernel l K uniquely define the Laplace transform of the solution. The solution itself is then uniquely found by inversion of its Laplace transform. It follows from (9) that the solution at time t depends on its history in the time interval [0, t]. The integral kernel K is often called memory function. Besides the memory effects, K describes how dynamics in the complementary subspace H 2 affect the dynamics in H 1 .
If all eigenvalues ζ k of the superoperator L 22 have negative real parts then any solution to (9) tends to the steadystate solutionσ On the other hand, in terms of (10), (12) and the Laplace transform, Due to the formulas we have that under the conditions we obtain the approximation where the second term is negligibly small. Since the small values of ζ are responsible for the long-term slow dynamics, it means that the steady-state part M of the integral term well approximates the informative long-term asymptotic of the dynamics. Thus, under the conditions we have the asymptotics in (12) Returning to the Laplace transform (11), this means that solutions to equation (9) are well approximated by solutions to the equationσ Equation (14) does not contain any memory function. Its operator L 11 + M is time-independent. Due to (13), fast dynamics in the complementary subspace H 2 make the solution in H 1 tend to the steady-state, rapidly "forgetting" its previous history. We call equation (14) an adiabatic approximation for equation (9) and say that the complementary subspace H 2 is adiabatically eliminated. The term adiabatic means that we first separate fast motions in H 2 from slow motions in H 1 and then eliminate all information about the fast motions, not influencing the slow dynamics. Here the equation preserves its initial form of the linear system with constant coefficients. The technique is easily generalised to inhomogeneous systems and initial conditions outside the informative subspace H 1 .
Such procedures are commonly used (among other fields) in condensed matter theory, quantum optics (where the term adiabatic elimination initially comes from), quantum statistics and quantum information, in cases where either external driving or an interaction with an environment are involved in such way that a non-resonant part of the dynamics becomes non-informative. A similar effect, well-known in the NMR context, is achieved by proceeding to a rotating frame and averaging out secular terms. Another example is eliminating high spin correlation orders in multi-spin dynamics where the role of the fast motion is played by spin relaxation. See in this context, for example, [12,13,31]. If the dimension of the informative subspace H 1 is much smaller than the dimension of the total space H, this gives a significant reduction of the initial problem.
Within the method described in the previous section, let the subspace H 1 be the subspace of zero-quantum coherences L 0 and H 2 the complementary subspace, Starting from the thermal equilibrium, we havê The relations along with (15) give Due to the presence of relaxation, the eigenvalues of the superoperator L 22 have negative real parts, satisfying the requirement outlined in the previous section. Physically, at high magnetic field, so the superoperator M is found as the series converging exponentially fast Using (16), we can calculate On the other hand, it follows from (16) Hence, because of (17), (18), (19), conditions (13) are satisfied, so we can replace equation (9) by its adiabatic approximationσ = L 0σ (20) with M 1,2 are defined in (19), and in (18) we restricted to the second order approximation for M.

Elimination of non-Zeeman spin orders
In the previous section, we adiabatically proceeded from the initial Liouvillian L to the new Liouvillian L 0 , well describing the dynamics closed in the subspace L 0 of zero-quantum coherences.
We decompose now the zero-quantum subspace as is the subspace of Zeeman spin orders, and L C is the complementary subspace consisting of non-Zeeman zero-quantum coherences.
The commutation character of the notationÔ ≡ [Ô, ·] implies that the superoperator M 1 is a commutation, We have alsoĤ The superoperatorsĤ 0,0 ,Ĥ 1,0 , Γ 0 trivially act on L Z . Within the accuracy of ∼ ω −1 I , only the superoperatorŝ H ′ 0 ,Ĥ ′ 1 contribute to L 21 , L 12 . The superoperator Thus, under the conditions min 2R conditions (13) are satisfied and the subspace L C is adiabatically eliminated. We come then to the adiabatic approximationσ closed in the Zeeman subspace L Z . In accordance with the previous notations, where π ′ 1,2 are projections onto the subspaces L Z,C respectively. Using conditions (17) , (22), it can be shown that the right-hand side of (23) is well approximated as where Γ ′ is given by (21), The advantage of formulas (24) is that they reduce the inversion L 0,22 −1 in the subspace L C to the much simpler problem of inversionsĈ ′ kj ,D ′ k of Zeeman operators. The latter are defined in the 2 n+1 -dimensional Hilbert space and are diagonal in the usual Zeeman basis.
The last terms in the expressions forĈ kj ,D k are the second order corrections with respect to the inverse nuclear Larmor frequency ω −1 I . For typical random spin geometries, these terms are either quenched by the first order terms or small simultaneously with them. Hence, these terms can often be neglected.

The Lindblad form
Introducing the notation and utilizing the double-commutator character of approximation (24) as well as the thermal correction Γ th , it is straightforward to see that the right-hand side of the Zeeman projected equation (23) is written in the purely Lindblad form related to the single-spin jump operatorsŜ ± ,Î k± , and the operator rateŝ (1 +Ĉ 2 kj ) −1 , related to the two-spin jump operatorŝ

Comparison between full Liouville space simulations and kMC simulations
We have carried out a number of comparisons between simulations of the polarization dynamics using the quantum master equation and our proposed method that involves an adiabatic approximation procedure to obtain a master equation for the states contained in the Zeeman subspace L Z and the use of a kinetic Monte Carlo (kMC) algorithm to calculate the polarization dynamics. As long as the parameters are chosen in accordance with the conditions described in the supplementary material, an excellent agreement between the two simulation methods is achieved. An example of such a comparison is given in Fig. (3) which shows the polarization dynamics in a system containing one electron and four nuclear spins forming a coupled network. The agreement can be further improved by increasing the number of averaged trajectories for the kMC method, however this is done on the expense of further increasing the computational time. The parameters used for the simulation are summarised in the following table.

SE DNP and Spin Diffusion
To evaluate the role that the nuclear dipolar interaction plays in the distribution of the polarization within a larger ensemble of coupled nuclear spins, we have simulated the polarization buildup dynamics in a linear chain of 40 coupled spins with an electron located at one end. The individual buildup curves for the 40 nuclear spins are shown in Fig. (4). The simulation parameters are listed in the following table (they are identical with the parameters used in [11]) . To investigate the role of the nuclear dipolar coefficient d kj for the distribution of the polarization within the linear spin chain, we calculated the average polarization buildup and compared it to the buildup obtained when the diplar coefficient is set to half its value for all nuclei k > 1 ( fig. (4)). In further simulations we either set all dipolar coefficients to half their value or ,in addition, we decreased the pseudo-secular interaction by a factor √ 2. These changes of the interaction parameters are in analogy to the set of simulations carried out by Hovav et al. in Ref. [11]. In contrast to their observations, the biggest effect on the average polarization is already brought about by halving the dipolar interaction strength for all nuclei k > 1, indicating the importance of this interaction for the diffusive transport. Any additional changes only introduced further small decreases of the average polarization.
We then compared a simple diffusion model ∂p/∂t = D av ∂ 2 p ∂x 2 , where p ≡ p(x, t) is the polarization field and a constant source is assumed for the start of the linear spin chain. The average diffusion constant was obtained from FIG. 4. Polarization buildup in a linear chain of 40 coupled nuclear spins with an electron located at one end. The buildup curve for each individual nucleus is shown. On the right hand we show the average polarization for the spin chain for the cases that i) all dipolar constants are given by the spacing between the individual nuclei, ii) the dipolar coupling constants d kj are halved for bulk nuclei, iii) the dipolar coupling constants are halved for all nuclei and iv) the dipolar couplings are halved for all nuclei and in addition the pseudo-secular interaction strength between the electron and the first nuclei is reduced by √ 2 eq. (7). This diffusion model was solved for two different boundary conditions, either with an absorbing boundary condition or a reflective boundary condition. The time course of the polarization moving through the linear spin chain is compared in Fig. (5) with a contour plot derived from simulating the polarizationchanges within the linear spin chain using our method described in the main part of the paper. It is clear evident that the dynamics is reasonably well approximated by the simple diffusion model if a reflective boundary condition is used. The parameters for this simulation are provided in the following table.