Phase transitions in electron spin resonance under continuous microwave driving

TL;DR

This paper investigates phase transitions in electron spin resonance under continuous microwave driving, revealing collective phenomena and critical behavior in strongly coupled electron systems interacting with an environment, using mean-field and numerical methods.

Contribution

It introduces a mean-field framework to analyze steady state phase transitions in driven, dissipative electron systems, supported by numerical simulations.

Findings

Identification of phase transitions between high and low polarization states.

Observation of metastable states and critical phenomena.

Good agreement between mean-field predictions and exact simulations.

Abstract

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.

Figures

Equations

$$H$$

$$\begin{array}[]{c}{\displaystyle\mathcal{D}=\sum_{k}\left[\gamma_{1+}{\cal L}(S_{k+})+\gamma_{1-}{\cal L}(S_{k-})+\gamma_{2}{\cal L}(S_{kz})\right],}\\ {\displaystyle\gamma_{1\pm}=\frac{R_{1}}{2}(1\mp p),\quad\gamma_{2}=2R_{2},\quad p=\tanh\frac{\hbar\omega_{S}}{2k_{\mathrm{B}}T}}\end{array}$$

$$\bar{H}=\sum _k\left({\omega }_1{S}_{kx}+\Delta S_{kz}\right)+\frac{3D}{N-1}\sum _{k<k^{\prime }}{S}_{kz}{S}_{k^{\prime }z}.$$

$${\bar{\Delta }}_k\in \delta (q)=\Delta +\frac{3D}{N-1}\left(q-\frac{N-1}{2}\right)$$

$$p_z^{\prime }(q)=p\left(1-\frac{\eta {\omega }_1^2}{{\delta }_0^2+{\delta }^2(q)}\right)$$

$${\bar{p}}_z=f({\bar{p}}_z)\equiv \sum _{q=0}^{N-1}P(q,p{\bar{p}}_z)p_z^{\prime }(q)/p.$$

$${\bar{p}}_z={\int }_{-\infty }^{+\infty }{f}_0(D,{\bar{p}}_z)\chi (D)dD.$$

$$\dot{\rho }=-i[H,\rho ]+\mathcal{D}\rho$$

$$H={\omega }_1{S}_x+\delta S_z,$$

$$\mathcal{D}=\frac{R_1}{2}\left[(1-p)\mathcal{L}(S_{+})+(1+p)\mathcal{L}(S_{-})\right]+2R_2\mathcal{L}(S_z).$$

$$\rho =1/2-p\left(X{S}_x+Y{S}_y+Z{S}_z\right),$$

$$\dot{X}=-\delta Y-R_{2}X,\dot{Y}=\delta X-{\omega }_{1}Z-R_{2}Y,$$

$$\dot{Z}={\omega }_{1}Y+R_1(1-Z).$$

$$X=\frac{{\omega }_1\delta }{R_2^2+{\delta }^2}Z,Y=-\frac{{\omega }_1{R}_2}{R_2^2+{\delta }^2}Z,$$

$$Z=1-\frac{{\omega }_1^2\eta }{{\delta }_0^2+{\delta }^2},{\delta }_0^2=R_2^2+{\omega }_1^2\eta ,\eta =\frac{R_2}{R_1}$$

$$u=f(a,b,u),$$

$$u-f(a,b,u)=\frac{\partial F}{\partial u}=0,F=\frac{u^2}{2}-\int f(a,b,u)du,$$

$$u^{\ast }+v=f(a,b,u^{\ast }+v)=f(a,b,u^{\ast })+\frac{\partial f}{\partial u}(a,b,u^{\ast })v+$$

$$\frac{1}{2}\frac{{\partial }^{2}f}{\partial u^2}(a,b,u^{\ast })v^2+\frac{1}{6}\frac{{\partial }^{3}f}{\partial u^3}(a,b,u^{\ast })v^3+\dots \equiv$$

$$c_0+c_{1}v+c_2{v}^2+c_3{v}^3+\dots$$

$$v=c_{1}v+c_2{v}^2+c_3{v}^3+\dots$$

$$0=c_2{v}^2+c_3{v}^3+\dots$$

$$0=c_3{v}^3+\dots ,$$

$$c_0=u^{\ast },c_1=1,c_2=0,c_3\ne 0$$

$$f(a^{\ast },b^{\ast },u^{\ast })=u^{\ast },\frac{\partial f}{\partial u}(a^{\ast },b^{\ast },u^{\ast })=1,$$

$$\frac{{\partial }^{2}f}{\partial u^2}(a^{\ast },b^{\ast },u^{\ast })=0,\frac{{\partial }^{3}f}{\partial u^3}(a^{\ast },b^{\ast },u^{\ast })\ne 0.$$

$$u^{\ast }+v=f(a^{\ast }+\alpha ,b^{\ast }+\beta ,u^{\ast }+v)\sim$$

$$u^{\ast }+{\xi }_0+(1+{\xi }_1)v+{\xi }_2{v}^2+{\xi }_3{v}^3+\dots$$

$$0\sim {\xi }_0+{\xi }_{1}v+{\xi }_2{v}^2+{\xi }_3{v}^3$$

$${\xi }_0=\frac{\partial f}{\partial a}\alpha +\frac{\partial f}{\partial b}\beta +\frac{1}{2}\frac{{\partial }^{2}f}{\partial a^2}{\alpha }^2+\frac{{\partial }^{2}f}{\partial a\partial b}\alpha \beta +\frac{1}{2}\frac{{\partial }^{2}f}{\partial b^2}{\beta }^2+$$

Full text

††preprint: APS/123-QED

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 Zavoisky (1945); Lancaster (1967); Schweiger and Jeschke (2001). 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. Griffin et al. (2010); Atsarkin and Köckenberger (2012) (ed.); Wenckebach (2016)) — in particular, this concerns the cross effect and thermal mixing DNP mechanisms Kessenikh et al. (1963); Hwang and Hill (1967); Hu et al. (2004); Borghini (1968); Atsarkin and Rodak (1972); Abragam and Goldman (1982); Karabanov et al. (2016).

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 Marcuzzi et al. (2014). 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 Ates et al. (2012); Rose et al. (2016); Foss-Feig et al. (2016), 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 Provotorov (1962); Borghini (1968); Atsarkin and Rodak (1972); Atsarkin (1978); Abragam and Goldman (1982); Jannin et al. (2012); Hovav et al. (2013); Serra et al. (2012); Luca and Rosso (2015).

Model — We model the evolution of the electron system within the framework of a Markovian Lindblad master equation. The density matrix $\rho$ of a system consisting of $N$ microwave-driven electrons evolves according to $\dot{\rho}=-i[H,\rho]+\mathcal{D}\rho$. The Hamiltonian $H$ at high static magnetic field, in the rotating frame approximation, is given by

[TABLE]

Here $\omega_{1}$ is the strength of the microwave field, $\Delta_{k}$ are the offsets of the electron Larmor frequencies (detunings) from the microwave carrier frequency, and $D_{kk^{\prime}}$, $D^{\prime}_{kk^{\prime}}$ are coefficients that parametrize the strength of the anisotropic and isotropic parts of the spin-spin dipolar and exchange interactions Schweiger and Jeschke (2001). Depending on the degree of order and symmetries within the sample structure, $D_{kk^{\prime}}$ and $D^{\prime}_{kk^{\prime}}$ can either be well defined (e.g., for crystals) or considered to be random (e.g., for glasses). In amorphous materials $\Delta_{k}$ are also distributed due to the anisotropic interaction of the electrons with the static field, leading to inhomogeneous broadening of the EPR line Schweiger and Jeschke (2001); Poole and Farach (1979); Karabanov et al. (2016).

Dissipative processes within the electron system are modeled by the dissipator $\mathcal{D}$ which describes single-spin relaxation and takes the form

[TABLE]

where ${\cal L}(X)\rho\equiv X\rho X^{\dagger}-\left\{X^{\dagger}X,\rho\right\}/2$ is the Lindblad form of a dissipation operator Karabanov et al. (2015). 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\in[0,1]$. Here, $p$ is a function of the average electron Larmor frequency $\omega_{S}$ and the temperature $T$. For typical experimental conditions ($W$-band, $\omega_{S}\sim 100~{}{\rm GHz}$, sample temperature between $T\sim 0~{}{\rm K}$ and $T\sim 100~{}{\rm K}$) the thermal spin polarization takes on values between $p\sim 1$ and $p\sim 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 $\Delta_{k}=\Delta$ and $D_{kk^{\prime}}=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

[TABLE]

We now compute the stationary average bulk polarization $p_{z}=-2\sum{\rm Tr}\,(S_{kz}\rho_{\rm ss})/N$ which serves as an order parameter for classifying the steady state $\rho_{\rm ss}$. To obtain the mean-field equation we define $H_{k}=\omega_{1}S_{kx}+\bar{\Delta}_{k}S_{kz}$ which is the projection of $\bar{H}$ onto the subspace of spin $k$. Here $\bar{\Delta}_{k}=\Delta+\frac{3D}{N-1}\sum_{k^{\prime}\not=k}S_{k^{\prime}z}$ is the effective energy shift or offset term experienced by the spin. This takes discrete values, i.e.,

[TABLE]

where $q=0,...,N-1$ is the number of spins $k^{\prime}\neq k$ in the up-state. The steady-state polarization $p^{\prime}_{z}(q)$ of a single spin for given $q$ is [see Appendix A]

[TABLE]

where $\delta_{0}=\sqrt{R_{2}^{2}+\eta\omega_{1}^{2}}$ and $\eta=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 $P(q,p_{z})=\binom{N-1}{q}\frac{(1-p_{z})^{q}(1+p_{z})^{N-1-q}}{2^{N-1}}$. Averaging Eq. (5) over all values of $q$ finally yields the equation for the relative steady-state polarization $\bar{p}_{z}=p_{z}/p$:

[TABLE]

Low and high temperature regime — The relative polarization is bounded ($|\bar{p}_{z}|\leq 1$), thus $f(\bar{p}_{z})$ defines a continuous map of the unit interval $\bar{p}_{z}\in[0,1]$ to itself. Therefore, by virtue of the Brouwer fixed point theorem Brouwer (1911), 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 polarization $\bar{p}_{z}$ on the average microwave offset $\Delta$, for three typical sets of parameters for $N=4$. Generally a good agreement is obtained. The observed spectra have $N$ Lorentzian peaks occurring at $\Delta=3D(1/2-q/(N-1))$, $q=0,1,\ldots,N-1$, with a half-width of $\delta_{0}$. The centre $\Delta=0$ of the spectrum corresponds to $q\sim q_{0}\equiv(N-1)/2$. The mean of the binomial distribution $P(q,p\,\bar{p}_{z})$ where the maximal saturation is given by $\bar{q}=(N-1)(1-p\,\bar{p}_{z})/2$. Here $\bar{q}$ 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\sim 0$) and undergo a shift from the centre at low temperatures ($p\sim 1$).

For large $N$ and high temperature we find a single broad region around $\Delta\sim 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\sim 1$ and high numbers of spins $N\gg 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=\Delta/\omega_{1}\sqrt{\eta}$, $b=3D/\omega_{1}\sqrt{\eta}$. FIG. 1(b) features a multi-stability 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 Marcuzzi et al. (2014); Carr et al. (2013); de Melo et al. (2016); Weller et al. (2016), or certain classes of dissipative Ising models Ates et al. (2012); Weimer (2015); Rose et al. (2016). The behavior of the steady-state polarization $\bar{p}_{z}$ upon crossing the multi-stable region is shown in FIG. 1(c).

Solutions with small $\bar{p}_{z}\sim 0$ correspond to non-thermal quasi-saturated equilibrium states. States with large values $\bar{p}_{z}\sim 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 Rose et al. (2016) 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 Carr et al. (2013); de Melo et al. (2016); Weller et al. (2016). 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 of $\bar{p}_{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 $|{\bar{p}_{z}}-{\bar{p}}_{\mathrm{crit}}|\sim y^{1/2}$, where ${\bar{p}}_{\mathrm{crit}}$ is the value of ${\bar{p}_{z}}$ at the critical point. Along the perpendicular direction we find $|{\bar{p}_{z}}-{\bar{p}}_{\mathrm{crit}}|\sim x^{1/3}$. These are Ising mean-field exponents Goldenfeld (1992). In the context of a classical Ising model, the directions $x$ and $y$ would correspond to magnetic field and temperature respectively (see also Ref. Marcuzzi et al. (2014)).

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 $\Delta$ 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 Karabanov et al. (2016). 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 $\Delta_{k}$ or the interactions $D_{kk^{\prime}}$ is large enough, unitary dynamics with Hamiltonian (1) is expected to undergo many-body localisation (MBL) Nandkishore and Huse (2015). 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 Levi et al. (2016); Medvedyeva et al. (2016); Fischer et al. (2016), suggesting that the mean-field analysis is still appropriate. (For other possible connections between MBL and DNP see Luca and Rosso (2015).)

For the sake of simplicity we assume that the interactions $D$ follow a Gaussian distribution, $\chi(D)=\exp(-D^{2}/D_{0}^{2})/(\sqrt{\pi}D_{0})$, with zero mean and standard deviation $D_{0}$. The offset frequency $\Delta$ may also be disordered (e.g., from the $g$-anisotropy and hyperfine interactions with nuclei Schweiger and Jeschke (2001); Poole and Farach (1979)), but we neglect that effect for now. Eq. (6) generalizes to

[TABLE]

Here we replaced the function $f(\bar{p}_{z})$ by its average with respect to the distribution $P(q,p\,\bar{p}_{z})$: $f_{0}(D,{\bar{p}_{z}})=p^{\prime}_{z}(\bar{q})/p=1-\eta\omega_{1}^{2}/(\delta_{0}^{2}+\delta^{2})$ with $\delta=\Delta-3Dp\,\bar{p}_{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\gg 1$, so Eq. (7) no longer depends on $N$ Rose et al. (2016); Ates et al. (2012); Marcuzzi et al. (2014) (see also Appendix D for details). The mean-field phase diagram resulting from Eq. (7) is displayed in FIG. 1(d) as a function of the dimensionless parameters $a^{\prime}=\Delta_{0}/\omega_{1}\sqrt{\eta}$ ($\Delta_{0}$ is the average offset, equal to $\Delta$ in the case considered here) and $b^{\prime}=3pD_{0}/2\omega_{1}\sqrt{\eta}$. We assume that the strength of the microwave field is large: $\omega_{1}^{2}\eta\gg 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^{\prime}$. 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}\propto b^{\prime}$), 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 $\Delta_{k}$. To this end the probability density $\chi(D)$ in Eq. (7) is replaced by a joint probability density $\chi(D,\Delta)$ accounting for both homogeneous and inhomogeneous broadening. The disorder in $\Delta_{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 Schirmer and Wang (2010); Weimer (2015). As shown in Rose et al. (2016); Ates et al. (2012); Foss-Feig et al. (2016), 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)\sum_{k}{\rm Tr}\,(S_{kz}\rho(t))$ for a variety of values of $a^{\prime}$ and $b^{\prime}$. For the set of parameters we consider, multiple disorder realizations of the dipolar coupling $\{D_{kk^{\prime}}\}$, with $D^{\prime}_{kk^{\prime}}=D_{kk^{\prime}}$ are taken. These are independent and identically distributed, sampled from a Gaussian distribution with variance defined by $b^{\prime}$ (see Appendix F for details). Fluctuations are quantified through the variance of the integrated polarization, $P_{z}=1/t\int_{0}^{t}p_{z}(t^{\prime})dt^{\prime}$. 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^{\prime}$ for several values of $a^{\prime}$, cf. FIG. 1(d). All curves display a peak indicating enhanced fluctuations for intermediate values of $b^{\prime}$, 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^{\prime}$, 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.

Appendix A Steady-state of single-spin microwave-driven dynamics

In the context of our work, the microwave-driven single-spin master equation has the form

[TABLE]

with

[TABLE]

[TABLE]

In terms of the relative polarization components

[TABLE]

we come to the Bloch equations (for $R_{2}\gg R_{1}$)

[TABLE]

[TABLE]

The steady-state solution where the right-hand sides are all zero is unique and calculated as

[TABLE]

[TABLE]

in full agreement with Eq. (5).

Appendix B Uniqueness of solution for

high temperatures and small

numbers of spins

To understand the structure of the solution space of Eq. (6) as a function of the thermal polarisation $p$ and the number of electrons $N$, we consider the derivative $df/d\bar{p}_{z}$: it is proportional to $p$, and thus for small values of $p$, corresponding to high temperatures, we have $df/d\bar{p}_{z}<1$. Under this condition the graph of the function $f(\bar{p}_{z})$ can intersect the diagonal $g(\bar{p}_{z})=\bar{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 $\max_{\bar{p}_{z}}df/d\bar{p}_{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 spectrum $\bar{p}_{z}(\Delta)$ 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 $\Delta=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

[TABLE]

can be studied from the point of view of the singularities in geometry of the 2-dimensional surface defined by the relation (8) in the 3-space $(a,b,u)$. The relation (8) can be rewritten as

[TABLE]

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 (8) known as the catastrophe theory Arnold (1984).

Consider the Taylor expansion of Eq. (8) near a given value $u=u^{*}$

[TABLE]

[TABLE]

[TABLE]

If $c_{0}\not=u^{*}$ then near the value $u=u^{*}$ Eq. (8) does not have solutions. If $c_{0}=u^{*}$ then $u=u^{*}$ is a solution, and we have

[TABLE]

If $c_{1}\not=1$ then the solution $u=u^{*}$ is locally unique. If $c_{1}=1$, $c_{2}\not=0$ then $u=u^{*}$ is a degeneracy point where two solutions merge,

[TABLE]

If $c_{2}=0$, $c_{3}\not=0$ then $u=u^{*}$ is a degeneracy point where three solutions merge,

[TABLE]

etc. Since relation (8) 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 so-called cusp point $G$ of the phase diagram Arnold (1984) which is defined by the critical values $a=a^{*}$, $b=b^{*}$, $u=u^{*}$ with

[TABLE]

which means

[TABLE]

[TABLE]

Consider now the Taylor expansion of Eq. (8) near the cusp point up to terms of the third order, taking into account Eq. (9),

[TABLE]

[TABLE]

which implies

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where the derivatives of $f$ are taken at $u=u^{*}$, $a=a^{*}$, $b=b^{*}$. The asymptotic cubic equation (10) has three solutions if $\bar{D}<0$ and has one solution if $\bar{D}>0$, where the discriminant $\bar{D}$ is given by the expression

[TABLE]

[TABLE]

where $\bar{D}_{n}$ is the term of the $n$th order in $\alpha$, $\beta$. The lowest order term is the quadratic term originated from $\xi_{0}^{2}$. This term forms the full square

[TABLE]

Making the rotation on the $(\alpha,\beta)$-plane

[TABLE]

and rewriting the cubic term $\bar{D}_{3}$ in the new parameters $x,\,y$, we obtain up to the third order

[TABLE]

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\xi_{3}^{2}>0$, so we can write

[TABLE]

[TABLE]

In other words, the critical curve $\bar{D}=0$ is asymptotically represented by the equation

[TABLE]

The last term can be removed by a shift transformation $x\to x+O(y^{2})$ and neglecting a term $\sim y^{4}$, so this curve is asymptotically written as

[TABLE]

This equation defines a cusp curve on the $(x,y)$-plane with two branches tangent to the $y$-axis at the cusp point $G$, see FIG.4(a) where the local geometry of the singular surface (11) is shown. In the rotated local coordinates, the cubic equation (10) representing the relation (8) takes the form

[TABLE]

Inside the cusp region $s_{0}x^{2}-s_{1}y^{3}<0$, Eq. (11) 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 solution $\bar{v}=0$ forks into three solutions $\bar{v}=0$ and $\bar{v}=\pm\,\bar{y}^{1/2}$. On crossing $G$ along the $x$-axis, the unique solution has a singularity $\bar{v}=x^{1/3}$. The described asymptotics are universal, i.e., valid for any two-parametric relation (8) as soon as it has a critical point where relations (9) hold Arnold (1984).

The critical point $G$ of the phase diagram of Eq. (6) satisfying Eq. (9) was found numerically to be

[TABLE]

with the characteristic directions in the $(a,b)$-plane

[TABLE]

[TABLE]

In FIG. 4(b), the structure of the solution $\bar{p}_{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^{\prime}$ of the phase diagram of Eq. (7) corresponds to

[TABLE]

with the characteristic directions (not plotted)

[TABLE]

[TABLE]

Appendix D Link to the classical meanfield

theory of the Ising model

As shown in the main text, the projection of the averaged Hamiltonian of Eq. (3) to the subspace of a randomly chosen spin $k$ is written as

[TABLE]

where

[TABLE]

The classical meanfield theory consists in replacing each operator $S_{k^{\prime}z}$ by its bulk steady-state observable (see, for example, Rose et al. (2016); Ates et al. (2012); Marcuzzi et al. (2014))

[TABLE]

This leads to the single-spin Hamiltonian

[TABLE]

Applying Eq. (5) justified in Appendix A, we obtain for the relative steady-state polarization

[TABLE]

Up to differences in notations, this is the classical self-consistent relation for the steady-state of the Ising model driven by a transversal field Rose et al. (2016); Ates et al. (2012); Marcuzzi et al. (2014).

The same result is obtained if we replace in Eq. (6) the summation over all $q$ by a single mean value of the binomial distribution $P(q,p\bar{p}_{z})$

[TABLE]

Indeed,

[TABLE]

To justify the proceeding from the whole set $q=0,\,1,\,\ldots,\,N-1$ to the mean $\bar{q}$, rescale the integer variable $q$ by a new variable $\epsilon$ by the rule

[TABLE]

where $\epsilon_{q}=0,\,1/(N-1),\,\ldots,\,1$ defines a uniform subdivision of the unit interval. The probability density of the variable $\epsilon_{q}$ is the same binomial distribution $P(q,p\bar{p}_{z})$ and the detuning $\delta$ becomes a function of $\epsilon$,

[TABLE]

Due to rescaling (13), the mean and the variance of the distribution $\epsilon_{q}$ are the mean and the variance of the distribution $P$ divided by $(N-1)$ and $(N-1)^{2}$ respectively, so we obtain

[TABLE]

[TABLE]

[TABLE]

In the limit $N\gg 1$, the variance $\sigma^{2}_{\epsilon}$ becomes zero, so the distribution $\epsilon_{q}$ is reduced to a single mean value $\bar{\epsilon}$ 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 delta-function $\tilde{\delta}(\epsilon-\bar{\epsilon})$,

[TABLE]

[TABLE]

[TABLE]

This justifies the classical meanfield theory (12) as a thermodynamic $N\gg 1$ limit of the meanfield theory developed in the main text.

Appendix E Effect of inhomogeneous

broadening

To estimate the effect of inhomogeneous broadening, we considered a system represented by two Gaussian spin packets of the same zero mean and standard deviation $D_{0}$ separated by a difference $2\Delta^{\prime}$ between the detunings. Here the Gaussian density $\chi(D)$ in Eq. (7) remains unchanged while the function $f_{0}(D,\bar{p}_{z})$ is modified as

[TABLE]

[TABLE]

The effect of $\Delta^{\prime}\not=0$ can be estimated varying the dimensionless parameter $\displaystyle c=\frac{\Delta^{\prime}}{\omega_{1}\sqrt{\eta}}$. For $c\not=0$, the phase diagram in the $(a^{\prime},b^{\prime})$-plane still features multi-stable regions but the latter are shifted and contracted with growing $c$. The contraction of the multistability region is explained by the fact that large differences between electron Larmor frequencies tend to quench the spin interactions and thus quench the multiplicity of the solution of the self-consistent relation Eq. (7).

Appendix F Quantum Jump Monte

Carlo simulations

The simulations for FIG. 2 were done using the Quantum Jump Monte Carlo algorithm Daley (2014) 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 Schirmer and Wang (2010)), i.e. after the relaxation time. The remainder of the trajectory is then cut up in to time periods $T$ of $\mathcal{O}({10}^{-2}\text{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 Rose et al. (2016)).

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^{\prime}$ 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_{D_{0}}(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}$.