Phase transitions in electron spin resonance under continuous microwave driving

We study an ensemble of strongly coupled electrons under continuous microwave irradiation interacting with a dissipative environment, a problem of relevance to the creation of highly polarized non-equilibrium states in nuclear magnetic resonance. We analyse the stationary states of the dynamics, described within a Lindblad master equation framework, at the mean-field approximation level. This approach allows us to identify steady state phase transitions between phases of high and low polarization controlled by the distribution of electronic interactions. We compare the mean-field predictions to numerically exact simulations of small systems and find good agreement. Our study highlights the possibility of observing collective phenomena, such as metastable states, phase transitions and critical behaviour in appropriately designed paramagnetic systems. These phenomena occur in a low-temperature regime which is not theoretically tractable by conventional methods, e.g., the spin-temperature approach.

Introduction -The control and detection of magnetization arising from a polarized ensemble of unpaired electron spins forms the basis of electron spin, or paramagnetic, resonance (ESR/EPR); a powerful spectroscopy tool for studying paramagnetic materials placed in a static external magnetic field. The underpinning key principle for this technique is the application of oscillating magnetic fields close to or at the electronic Larmor frequency (usually in the microwave regime) to generate non-equilibrium distributions of populations and coherences between quantum states that lead to detectable signals [1][2][3]. The evolution of systems of isolated or only weakly coupled paramagnetic centres under the effect of these fields is well understood. A more challenging problem is to predict the response of strongly coupled electron ensembles to such perturbations, particularly in samples in the solid state in which anisotropic components of the electronic interactions are not averaged out by thermal motion. Insight into the dynamics of strongly coupled, microwave driven electronic ensembles is also needed in order to improve our understanding of dynamic nuclear polarization (DNP), which is an out-of-equilibrium technique to enhance the sensitivity of nuclear magnetic resonance (NMR) applications by orders of magnitude (see, e.g., Ref. [4][5][6]) -in particular, this concerns the cross effect and thermal mixing DNP mechanisms [7][8][9][10][11][12][13].
Here we shed light on the non-equilibrium stationary states of a strongly interacting electronic ensemble under continuous microwave driving and subject to dissipation to the environment. We model the dynamics of this system in terms of a Markovian master equation and use a mean-field approximation to compute the steady state phase diagram. This reveals phase transitions between states of high and low electronic polarisation as well as the emergence of a critical point that displays Ising universality [14]. These features are controlled by the distribution of the disordered electronic spin-spin interactions. The uncovered mean-field transitions imply the emergence of metastable states and accompanying intermittent dynamics [15][16][17], which we confirm numerically through simulations of small systems. Our results suggest that under appropriate conditions collective phenomena such as metastability, phase transitions and critical behaviour should be observable in driven-dissipative, paramagnetic systems. These predictions complement those of conventional theoretical approaches, based, e.g., on the so-called spin-temperature which, due to their restriction to certain parameter regimes, would only predict a homogenous quasi-equilibrium state [10][11][12][18][19][20][21][22][23]. Model -We model the evolution of the electron system within the framework of a Markovian Lindblad master equation. The density matrix ρ of a system consisting of N microwave-driven electrons evolves according tȯ ρ = −i[H, ρ] + Dρ. The Hamiltonian H at high static magnetic field, in the rotating frame approximation, is given by Here ω 1 is the strength of the microwave field, ∆ k are the offsets of the electron Larmor frequencies (detunings) from the microwave carrier frequency, and D kk , D kk are coefficients that parametrize the strength of the anisotropic and isotropic parts of the spin-spin dipolar and exchange interactions [3]. Depending on the degree of order and symmetries within the sample structure, D kk and D kk can either be well defined (e.g., for crystals) or considered to be random (e.g., for glasses). In amorphous materials ∆ k are also distributed due to the anisotropic interaction of the electrons with the static field, leading to inhomogeneous broadening of the EPR line [3,13,24]. Dissipative processes within the electron system are  7) in the (a , b )-plane featuring regions of unique and multiple solutions similar to that in panel (b) and a critical point G belonging to the same universality class as G (see text for details). The dark gray region illustrates the contraction of the multi-stability region caused by inhomogeneous broadening (see text and Appendix C for details). modeled by the dissipator D which describes single-spin relaxation and takes the form where L(X)ρ ≡ XρX † − X † X, ρ /2 is the Lindblad form of a dissipation operator [25]. The dissipation rates depend on the longitudinal (R 1 ) and transversal (R 2 ) relaxation rates of the electron spins as well as the thermal polarization p ∈ [0, 1]. Here, p is a function of the average electron Larmor frequency ω S and the temperature T . For typical experimental conditions (W -band, ω S ∼ 100 GHz, sample temperature between T ∼ 0 K and T ∼ 100 K) the thermal spin polarization takes on values between p ∼ 1 and p ∼ 10 −2 . Mean-field in the absence of disorder -In order to obtain a basic understanding of the phase structure of the driven electron system, let us first disregard any dispersion in the frequency offsets and interactions, by setting ∆ k = ∆ and D kk = D/(N − 1). (Note, that this N -dependence takes into account the fact that in practice the coupling strengths decay rapidly with the interspin distance and keep the interaction energy an extensive quantity.) In the non-disordered case, the last term of Eq.
(1) commutes with the rest of the Hamiltonian and does not influence the bulk polarization dynamics. Therefore, we can neglect it, leading to the mean-field Hamiltonian We now compute the stationary average bulk polarization p z = −2 Tr (S kz ρ ss )/N which serves as an order parameter for classifying the steady state ρ ss . To obtain the mean-field equation we define H k = ω 1 S kx +∆ k S kz which is the projection ofH onto the subspace of spin k. Here∆ k = ∆ + 3D N −1 k =k S k z is the effective energy shift or offset term experienced by the spin. This takes discrete values, i.e., where q = 0, ..., N −1 is the number of spins k = k in the up-state. The steady-state polarization p z (q) of a single spin for given q is [see Appendix A] where δ 0 = R 2 2 + ηω 2 1 and η = R 2 /R 1 is the ratio of the electron spin relaxation rates. Since the system is homogeneous, the steady-state polarization of the individual spins is identical and given by p z , which can be regarded as a self-consistency condition. Hence, the probability of having q up spins and N − q − 1 down spins is given by . Averaging Eq. (5) over all values of q finally yields the equation for the relative steady-state polarizationp z = p z /p: Low and high temperature regime -The relative polarization is bounded (|p z | ≤ 1), thus f (p z ) defines a continuous map of the unit intervalp z ∈ [0, 1] to itself. Therefore, by virtue of the Brouwer fixed point theorem [26], Eq. (6) always has at least one solution. We find that the solution is unique for small values of p corresponding to high temperatures and small numbers of spins N (see Appendix B).
For small values of N we can compare the results of the mean-field treatment to the exact solution of the quantum master equation given by the dissipator (2) and Hamiltonian (3). To this end we show in FIG. 1(a) the steady-state polarization spectrum, i.e. the dependence of the bulk polarizationp z on the average microwave offset ∆, for three typical sets of parameters for N = 4. Generally a good agreement is obtained. The observed spectra have N Lorentzian peaks occurring at ∆ = 3D(1/2 − q/(N − 1)), q = 0, 1, . . . , N − 1, with a half-width of δ 0 . The centre ∆ = 0 of the spectrum corresponds to q ∼ q 0 ≡ (N − 1)/2. The mean of the binomial distribution P (q, pp z ) where the maximal saturation is given byq = (N − 1)(1 − pp z )/2. Hereq is close to q 0 for small p and tends to shift from q 0 with increasing p. Hence, the intensities of the peaks are symmetric with respect to the centre of the spectrum at high temperatures (p ∼ 0) and undergo a shift from the centre at low temperatures (p ∼ 1).
For large N and high temperature we find a single broad region around ∆ ∼ 0, in which the polarization is saturated due to the applied field. The width of this region increases with interaction strength D (see Appendix B for details).
Multi-stability and phase transitions -The situation qualitatively changes when entering the regime of low temperatures p ∼ 1 and high numbers of spins N 1. In this case (see Appendix B) Eq. (6) can feature more than one solution. In FIG. 1(b) we show the phase diagram given by the number of solutions of Eq. (6) in terms of the scaled offset and interaction parameters a = ∆/ω 1 √ η, b = 3D/ω 1 √ η. FIG. 1(b) features a multistability region where three solutions coexist (gray) separated from the regions with a unique solution (brown) by two spinodal lines that coalesce at a critical point G. Similar phase diagrams have recently been discussed theoretically in other contexts, e.g., for open driven gases of strongly interacting Rydberg atoms [14,[27][28][29], or certain classes of dissipative Ising models [15,16,30]. The behavior of the steady-state polarizationp z upon crossing the multi-stable region is shown in FIG. 1(c). Solutions with smallp z ∼ 0 correspond to non-thermal quasi-saturated equilibrium states. States with large valuesp z ∼ 1 are unsaturated quasi-thermal equilibria. On crossing the spinodal curve 1 from large negative values of a, the unique stable quasi-thermal steadystate continues to exist but two other steadystate solutions appear: a stable quasi-saturated and an unstable intermediate one as shown in FIG. 1(c). Conversely, on crossing curve 2 towards large negative values of a, the unique stable quasi-saturated steadystate continues to exist but two other steadystates emerge, a stable and an unstable one.
The occurrence of multiple steady state solutions is an artifact of the mean-field approximation. It can be interpreted as the emergence of metastable states [16] near first-order phase transitions. An experimental signature of this type of physics is for example hysteretic behavior as recently studied in the context of interacting atomic gases [27][28][29]. We will return to this point further below.
The nature of the critical point G in the phase diagram FIG. 1(b) can be characterized by analysing the scaling behaviour ofp z near it. We find two directions that are singled out (see Appendix C for details): one is given by the curve that is tangent to both spinodal lines [see FIG. 1(b)], where we find |p z −p crit | ∼ y 1/2 , wherē p crit is the value ofp z at the critical point. Along the perpendicular direction we find |p z −p crit | ∼ x 1/3 . These are Ising mean-field exponents [31]. In the context of a classical Ising model, the directions x and y would correspond to magnetic field and temperature respectively (see also Ref. [14]). Disordered spin-spin interactions and augmented mean-field -The results so far indicate possible phase transitions in the polarization of the electron system controlled by the frequency offset ∆ and the average interaction strength D. However, typical sample materials are not single crystals and electrons are arranged randomly, such that the average interaction experienced by an electron is close to zero [13]. In order to take this into account we need an augmented mean-field description which accounts for a distribution in the coupling strengths.
Note that when the disorder in either the offsets ∆ k or the interactions D kk is large enough, unitary dynamics with Hamiltonian (1) is expected to undergo many-body localisation (MBL) [32]. In this case spatial fluctuations in the long-time state can be significant and determined by the disorder and the initial state, which raises the question of the appropriateness of mean-field. However, in the presence of dissipation, cf. Eq. (2), MBL is unstable and the stationary state is delocalised [33][34][35], suggesting that the mean-field analysis is still appropriate. (For other possible connections between MBL and DNP see [23].) For the sake of simplicity we assume that the interactions D follow a Gaussian distribution, χ(D) = exp(−D 2 /D 2 0 )/( √ πD 0 ), with zero mean and standard deviation D 0 . The offset frequency ∆ may also be disordered (e.g., from the g-anisotropy and hyperfine interactions with nuclei [3,24]), but we neglect that effect for now. Eq. (6) generalizes tō Here we replaced the function f (p z ) by its average with respect to the distribution P (q, pp z ): f 0 (D,p z ) = p z (q)/p = 1 − ηω 2 1 /(δ 2 0 + δ 2 ) with δ = ∆ − 3Dpp z /2. This is justified by the properties of the distribution P and by the fact that the averaged function f 0 coincides with the classical mean-field approximation of the Ising model in the limit of N 1, so Eq. (7) no longer depends on N [14-16] (see also Appendix D for details). The meanfield phase diagram resulting from Eq. (7) is displayed in FIG. 1(d) as a function of the dimensionless parameters a = ∆ 0 /ω 1 √ η (∆ 0 is the average offset, equal to ∆ in the case considered here) and b = 3pD 0 /2ω 1 √ η.
We assume that the strength of the microwave field is large: ω 2 1 η R 2 2 meaning that the electron system is fully saturated in the absence of spin-spin coupling (in which case the phase transitions observed are most pronounced). The structure is similar to that of FIG. 1(b). We observe regions with one and three solutions as well as spinodal lines forming a cusp at a critical point G . The scaling properties at this critical point are, again, those of mean-field Ising universality. Note, however, that the phase transition is controlled by the width of the distribution of the disorder strengths (D 0 ∝ b ), rather than the average interaction strength, which is in fact zero.
Eq. (7) can be modified to take into account disorder in the frequency offsets ∆ k . To this end the probability density χ(D) in Eq. (7) is replaced by a joint probability density χ(D, ∆) accounting for both homogeneous and inhomogeneous broadening. The disorder in ∆ k causes a shift and contraction of the multi-stability region which is illustrated by the dark gray region in FIG.1(d) where the dimensionless parameter c characterizes inhomogeneous broadening (see Appendix E for details).
Fluctuations and numerical simulations -The mean-field treatment above is of course not exact. Whether the predicted qualitative phase structure survives away from mean-field depends on the effect of fluctuations [30,36]. As shown in [15][16][17], phase coexistence at the mean-field level can be an indication -away from the thermodynamic limit -of the existence of long-lived metastable (rather than stationary) phases. These competing phases come with an intermittent dynamics of slow switching between them. We now show that this is indeed the case by considering the dynamics of the exact system, Eqs. (1), (2), by means of numerical simulations in small systems.
We study the time dependence of the polarization p z (t) = −(2/N ) k Tr (S kz ρ(t)) for a variety of values of a and b . For the set of parameters we consider, multiple disorder realizations of the dipolar coupling {D kk }, with D kk = D kk are taken. These are independent and identically distributed, sampled from a Gaussian distribution with variance defined by b (see Appendix F for details). Fluctuations are quantified through the variance of the integrated polarization, P z = 1/t t 0 p z (t )dt . In our simulations t is chosen long enough, such that fluctuations due to the transient, short time dynamics average out. In FIG. 2(a) we show the disorder averaged variance of P z as a function of b for several values of a , cf. FIG. 1(d). All curves display a peak indicating enhanced fluctuations for intermediate values of b , which is the region where metastable states and enhanced fluctuations are expected. Similar behaviour is observed in FIG. 2(b-d) for the probability distribution of P z , shown both for individual disorders and averaged over disorder.
Since the system is small, we do not expect self-averaging, and P z for individual realisations of the disorder to vary. Nevertheless, all histograms broaden significantly for intermediate values of b , clearly displaying enhanced fluctuations as indicated by the multi-stable region identified by our mean-field analysis.
Conclusions -Our results demonstrate that cooperative behaviour in strongly interacting ensembles of microwave driven electrons -a situation of relevance to DNP in NMR -can give rise to a non-trivial phase structure in the stationary state of these systems. Mean-field analysis predicts the existence of phases of distinct polarisation, with phase transitions between them controlled by the detunings in the microwave driving and the distribution of the dipolar electronic couplings. While the calculated phase diagram is mean-field in origin, our simulations show that -even for finite systems -dynamics will be correlated and intermittent, related to the coexistence of metastable states. The experimental demonstration of these predicted phenomena would ideally require a paramagnetic sample with minimal inhomogeneous broadening, kept at cryogenic temperatures and high magnetic field. In the context of our work, the microwave-driven single-spin master equation has the forṁ

ACKNOWLEDGMENTS
In terms of the relative polarization components we come to the Bloch equations (for R 2 R 1 ) The steady-state solution where the right-hand sides are all zero is unique and calculated as in full agreement with Eq. (5). df /dp z : it is proportional to p, and thus for small values of p, corresponding to high temperatures, we have df /dp z < 1. Under this condition the graph of the function f (p z ) can intersect the diagonal g(p z ) =p z only once and hence Eq. (6) has only one solution. This high temperatures behaviour is independent of the number of spins N , which is illustrated in FIG. 3(a). Here we plot maxp z df /dp z as function of p for different values of N and fixed other parameters, showing that the maximum slope for small p is always negative. The shape of the steady-state polarization spectrump z (∆) is described and good agreement between the master equation and the meanfield Eq. (6) for small N is illustrated in the main text. In FIG. 3(b) we show the high-temperature steady-state polarization spectrum resulting from Eq. (6) for large N and different values of D. Broadening of the saturation region around ∆ = 0 with increasing D is evident.
Appendix C: Structure of the phase diagram Mathematically, the phase diagram of a (smooth) general two-parametric family of self-consistent relations of the form can be studied from the point of view of the singularities in geometry of the 2-dimensional surface defined by the relation (C1) in the 3-space (a, b, u). The relation (C1) can be rewritten as which defines a critical point u of a (smooth) scalar function F (u) depending on the parameters a, b. This makes a subject of the mathematical theory of singularities combined with the geometry of the surface (C1) known as the catastrophe theory [37].
Consider the Taylor expansion of Eq. (C1) near a given value u = u * If c 0 = u * then near the value u = u * Eq. (C1) does not have solutions. If c 0 = u * then u = u * is a solution, and we have If c 1 = 1 then the solution u = u * is locally unique. If c 1 = 1, c 2 = 0 then u = u * is a degeneracy point where two solutions merge, If c 2 = 0, c 3 = 0 then u = u * is a degeneracy point where three solutions merge, etc. Since relation (C1) depends on two parameters a, b and one variable u, in a generic situation no more than three conditions on the coefficients c 0 , c 1 , c 2 can be simultaneously satisfied, so not more than three solutions can merge at u = u * . The latter takes place at the socalled cusp point G of the phase diagram [37] which is defined by the critical values a = a * , b = b * , u = u * with Consider now the Taylor expansion of Eq. (C1) near the cusp point up to terms of the third order, taking into account Eq. (C2), where the derivatives of f are taken at u = u * , a = a * , b = b * . The asymptotic cubic equation (C3) has three solutions ifD < 0 and has one solution ifD > 0, where the discriminantD is given by the expression whereD n is the term of the nth order in α, β. The lowest order term is the quadratic term originated from ξ 2 0 . This term forms the full squarē Making the rotation on the (α, β)-plane and rewriting the cubic termD 3 in the new parameters x, y, we obtain up to the third order where the coefficients s 0−4 are expressed via the derivatives of the function f (a, b, u) at the cusp point. We have s 0 = (r 2 + t 2 )/4ξ 2 3 > 0, so we can writē In other words, the critical curveD = 0 is asymptotically represented by the equation The last term can be removed by a shift transformation x → x + O(y 2 ) and neglecting a term ∼ y 4 , so this curve is asymptotically written as s 0 x 2 − s 1 y 3 = 0 : y = s 0 s 1 Inside the cusp region s 0 x 2 − s 1 y 3 < 0, Eq. (C4) has three solutions, outside the cusp region s 0 x 2 − s 1 y 3 > 0 only one solution exists. On crossing the cusp point G along the y-axis, the unique solutionv = 0 forks into three solutionsv = 0 andv = ±ȳ 1/2 . On crossing G along the x-axis, the unique solution has a singularitȳ v = x 1/3 . The described asymptotics are universal, i.e., valid for any two-parametric relation (C1) as soon as it has a critical point where relations (C2) hold [37]. The critical point G of the phase diagram of Eq. (6) satisfying Eq. (C2) was found numerically to be a * ∼ −0.18, b * ∼ 3.23,p crit ∼ 0.27 with the characteristic directions in the (a, b)-plane x ∼ 0.99(a − a * ) − 0.14(b − b * ), y ∼ 0.99(b − b * ) + 0.14(a − a * ).
In FIG. 4(b), the structure of the solutionp z is shown on crossing the critical point G along the tangent direction y, in FIG. 4(c) -the same on crossing along the perpendicular direction x.
The critical point G of the phase diagram of Eq. (7) corresponds to a * ∼ 0.26, b * ∼ 3.83,p crit ∼ 0.20 with the characteristic directions (not plotted) To justify the proceeding from the whole set q = 0, 1, . . . , N −1 to the meanq, rescale the integer variable q by a new variable by the rule where q = 0, 1/(N − 1), . . . , 1 defines a uniform subdivision of the unit interval. The probability density of the variable q is the same binomial distribution P (q, pp z ) and the detuning δ becomes a function of , Due to rescaling (D2), the mean and the variance of the distribution q are the mean and the variance of the distribution P divided by (N − 1) and (N − 1) 2 respectively, so we obtain In the limit N 1, the variance σ 2 becomes zero, so the distribution q is reduced to a single mean value¯ taken with the probablity 1. The summation over q can be replaced by an integration over the unit interval with the probablity density represented by the Dirac deltafunctionδ( −¯ ), This justifies the classical meanfield theory (D1) as a thermodynamic N 1 limit of the meanfield theory developed in the main text.
The simulations for FIG. 2 were done using the Quantum Jump Monte Carlo algorithm [38] to calculate the stochastic evolution (trajectory) of the pure state of the system over time. While all trajectories are initialized in the same state, the all up configuration, data from a trajectory is only considered after sufficient time has elapsed that there is no memory of the initial state (we can be certain such a time scale exists for this finite system due to the results of [36]), i.e. after the relaxation time. The remainder of the trajectory is then cut up in to time periods T of O(10 −2 s), chosen such that short time fluctations are averaged out so that only long time fluctuations influence the variance of the time integrated observable (similar to the approach used in Sec. III E of [16]).
Different disorder realizations are handled as follows: we begin by taking a set of random numbers from a Gaussian distribution of unit variance, defining the realization. For a given value of b we then rescale all of these numbers by the associated value of the standard deviation D 0 . As it can be shown that the probability density satisfies p 1 (x)dx = p D0 (D 0 x)d(D 0 x) where the subscript represents the variance of the Gaussian, this rescaling provides us with an equivalent set of numbers that were effectively drawn from a distribution with standard deviation D 0 .