# Dissipative Phase Transition in Central Spin Systems

## Abstract

We investigate dissipative phase transitions in an open central spin system. In our model the central spin interacts coherently with the surrounding many-particle spin environment and is subject to coherent driving and dissipation. We develop analytical tools based on a self-consistent Holstein-Primakoff approximation that enable us to determine the complete phase diagram associated with the steady states of this system. It includes first and second-order phase transitions, as well as regions of bistability, spin squeezing and altered spin pumping dynamics. Prospects of observing these phenomena in systems such as electron spins in quantum dots or NV centers coupled to lattice nuclear spins are briefly discussed.

###### pacs:

## I Introduction

Statistical Mechanics classifies phases of a given system in thermal equilibrium according to its physical properties. It also explains how changes in the system parameters allow us to transform one phase into another, sometimes abruptly, which results in the phenomenon of phase transitions. A special kind of phase transitions occurs at zero temperature: such transitions are driven by quantum fluctuations instead of thermal ones and are responsible for the appearance of exotic quantum phases in many areas of physics. These quantum phase transitions have been a subject of intense research in the last thirty years, and are expected not only to explain interesting behavior of systems at low temperature, but also to lead to new states of matter with desired properties (e.g., superconductors, -fluids and -solids, topological insulators Vojta (2003); Sachdev (2003); Ginzburg (2007); Kim and Chan (2004); Belitz and Kirkpatrick (1994); Hasan and Kane (2010)).

Phase transitions can also occur in systems away from their thermal equilibrium. For example, this is the case when the system interacts with an environment and, at the same time, is driven by some external coherent source. Due to dissipation, the environment drives the system to a steady state, , which depends on the system and environment parameters, . As is changed, a sudden change in the system properties may occur, giving rise to a so called dissipative phase transition (DPT) Carmichael (1980); Werner et al. (2005); Capriotti et al. (2005); Morrison and Parkins (2008); Eisert and Prosen (2010); Bhaseen et al. (2012). DPT have been much less studied than traditional or quantum ones. With the advent of new experimental techniques that allow to observe them experimentally, they are starting to play an important role Baumann et al. (2010). Moreover they offer the intriguing possibility of observing critical effects non-destructively because of the constant intrinsic exchange between system and environment Öztop et al. (2010). In equilibrium statistical mechanics a large variety of toy models exist that describe different kind of transitions. Their study lead to a deep understanding of many of them. In contrast, in the case of DPT few models have been developed.

The textbook example of a DPT occurs in the Dicke model of resonance fluorescence Hepp and Lieb (1973); Carmichael (1980). There, a system of spins interacts with a thermal reservoir and is externally driven. Experimental Gibbs et al. (1976) and theoretical studies Lawande et al. (1981); Puri et al. (1980); Bonifacio and Lugiato (1978); Schneider and Milburn (2002) revealed interesting features such as optical multistability, first and second order phase transitions, and bipartite entanglement.

In this paper, we analyze another prototypical open system: The model is closely related to the central spin system (CSS) which has been thoroughly studied in thermal equilibrium Gaudin (1976); Bortz and Stolze (2007); Schliemann et al. (2003). In its simplest form, it consists of a set of spin-1/2 particles (in the following referred to as the nuclear spins), uniformly coupled to a single spin-1/2 (referred to as the electron spin). In the model we consider, the central spin is externally driven and decays through interaction with a Markovian environment. Recently, the CSS model has found application in the study of solid state systems such as electron and nuclear spins in a quantum dot Schliemann et al. (2003) or an Nitrogen-Vacancy-center.

In what follows we first provide a general framework for analyzing DPT in open systems. In analogy with the analysis of low energy excitations for closed systems, it is based on the study of the excitation gap of the system’s Liouville operator . We illustrate these considerations using the central spin model. For a fixed dissipation strength , there are two external parameters one can vary, the Rabi frequency of the external driving field, , and the Zeemann shift, . We present a complete phase diagram as a function of those parameters, characterize all the phases, and analyze the phase transitions occurring among them. To this end, we develop a series of analytical tools, based on a self-consistent Holstein-Primakoff approximation, which allows us to understand most of the phase diagram. In addition, we use numerical methods to investigate regions of the diagram where the theory yields incomplete results. Combining these techniques, we can identify two different types of phase transitions, and regions of bistability, spin squeezing, and enhanced spin polarization dynamics. We will also identify regions where anomalous behavior occurs in the approach to the steady state. Intriguingly, recent experiments with quantum dots, in which the central (electronic) spin is driven by a laser and undergoes spontaneous decay, realize a situation very close to the one we study here and show effects such as bistability, enhanced fluctuations, and abrupt changes in polarization in dependence of the system parameters Krebs et al. (2010); Russell et al. (2007).

This paper is organized as follows. Section II sets the general theoretical framework underlying our study of DPT. Section III introduces the model, and contains a structured summary of the main results. In Section IV we develop the theoretical techniques and use those techniques to analyze the various phases and classify the different transitions. Thereafter in Section V numerical techniques are employed to explain the features of the phase diagram which are not captured by the previous theory. Possible experimental realizations and a generalization of the model to inhomogeneous coupling are discussed in Section VI. Finally we summarize the results and discuss potential applications in Section VII.

## Ii General Theoretical Framework

The theory of quantum phase transitions in closed systems is a well established and extensively studied area in the field of statistical mechanics. The typical scenario is the following: a system is described by a Hamiltonian, , where denotes a set of systems parameters (like magnetic fields, interactions strengths, etc). At zero temperature and for a fixed set of parameters, , the system is described by a quantum state, , fulfilling , where is the ground state energy. As long as the Hamiltonian is gapped (i.e., the difference between and the first excitation energy is finite), any small change in will alter the physical properties related to the state smoothly and we remain in the same phase. However, if the first excitation gap closes at a given value of the parameters, , it may happen that the properties change abruptly, in which case a phase transition occurs.

TPT | QPT | DPT | |

System | Hamiltonian | Hamiltonian | Liouvillian |

operator | – Lindblad | ||

Relevant | Free energy | Energy eigenvalues | ”Complex energy” eigenvalues |

quantity | |||

Gibbs state | Ground state | Steady state | |

State | |||

Phase transition | Non-analyticity in | vanishes | vanishes |

In the following we adapt analogous notions to the case of DPT and introduce the concepts required for the subsequent study of a particular example of a generic DPT in a central spin model.

We consider a Markovian open system, whose evolution is governed by a time-independent master equation . The dynamics describing the system are contractive implying the existence of a steady state. This steady state is a zero eigenvector to the Liouville superoperator . This way of thinking parallels that of quantum phase transitions, if one replaces . Despite the fact that these mathematical objects are very different (the first is a Hermitian operator, and the second a hermiticity-preserving superoperator), one can draw certain similarities between them. For instance, for an abrupt change of (and thus of certain system observables) it is necessary that the gap in the (in general complex) excitation spectrum of the system’s Liouville operator closes. The relevant gap in this context is determined by the eigenvalue with largest real part different from zero (it can be shown that for all eigenvalues of Riv (2012)). The vanishing of the real part of this eigenvalue – from here on referred to as asymptotic decay rate (ADR) Horstmann (2011) – indicates the possibility of a non-analytical change in the steady state and thus is a necessary condition for a phase transition to occur.

In our model system, the Liouvillian low excitation spectrum, and the ADR in particular, can in large parts of the phase diagram be understood from the complex energies of a stable Gaussian mode of the nuclear field. We find first order transitions where the eigenvalue of this stable mode crosses the eigenvalue of a metastable mode at zero in the projection onto the real axis. The real part of the Liouvillian spectrum closes directly as the stable mode turns metastable and vice versa. A finite difference in the imaginary parts of the eigenvalues across the transition prevents a mixing of the two modes and the emergence of critical phenomena such as a change in the nature of the correlations in the steady state at the critical point. In contrast, we also find a second order phase transition where the ADR vanishes asymptotically as both mode energies become degenerate (at zero) in the thermodynamic limit. Mixing of the two modes at the critical point gives rise to diverging correlations in the nuclear system. This observation parallels the classification of quantum phase transitions in closed systems. There, a direct crossing of the ground and first excited state energy for finite systems (mostly arising from a symmetry in the system) typically gives rise to a first order phase transition. An asymptotical closing of the first excitation gap of the Hamiltonian in the thermodynamic limit represents the generic case of a second order transition Sachdev (1999).

Besides the analogies described so far [cf. Table 1], there are obvious differences, like the fact that in DTP may be pure or mixed, and that some of the characteristic behavior of a phase may also be reflected in how the steady state is approached. Non-analyticities in the higher excitation spectrum of the Liouvillian are associated to such dynamical phases.

## Iii Model and Phase Diagram

### iii.1 The Model

We investigate the steady state properties of a homogeneous central spin model. The central spin – also referred to as electronic spin in the following – is driven resonantly via suitable optical or magnetic fields. Dissipation causes electronic spin transitions from the spin-up to the spin-down state. It can be introduced via standard optical pumping techniques. Atatüre et al. (2006); Tamarat et al. (2008). Furthermore, the central spin is assumed to interact with an ensemble of ancilla spins – also referred to as nuclear spins in view of the mentioned implementations Schliemann et al. (2003)– by an isotropic and homogeneous Heisenberg interaction. In general this hyperfine interaction is assumed to be detuned. Weak nuclear magnetic dipole-dipole interactions are neglected.

After a suitable transformation which renders the Hamiltonian time-independent, the system under consideration is governed by the master equation

(1) | ||||

where denotes the anticommutator and

(2) | ||||

(3) | ||||

(4) |

and () denote electron and collective nuclear spin operators, respectively.
is the Rabi frequency of the resonant external driving of the electron (in rotating wave approximation), while is the difference of hyperfine detuning and half the individual hyperfine coupling strength . for instance can be tuned via a static magnetic fields in direction.
Note that , describing the isotropic hyperfine interaction and its detuning.
The rescaling of the electron driving and dissipation in terms of the total (nuclear) spin quantum number
^{1}

In the limit of strong dissipation the electron degrees of freedom can be eliminated and Eq. (1) reduces to

(5) | ||||

where , and is the reduced density matrix of the nuclear system. This is a generalization of the Dicke model of resonance fluorescence as discussed in Carmichael (1980); Schneider and Milburn (2002); Morrison and Parkins (2008).

Master Eq. (1) has been theoretically shown to display cooperative nuclear effects such as superradiance even for inhomogeneous electron nuclear coupling Kessler et al. (2010). In analogy to the field of cooperative resonance fluorescence, the system’s rich steady state behavior comprises various critical effects such as first and second order DPT and bistabilities. In the following we provide a qualitative summary of the phase diagram and of the techniques developed to study the various phases and transitions.

### iii.2 Phenomenological Description of the Phase Diagram

For a fixed dissipation rate
^{2}

First we consider the system along the line segment (), where define a critical driving strength and critical hyperfine detuning, respectively. Here vanishes and the steady state can be constructed analytically as a zero-entropy factorized state of the electron and nuclear system. The nuclear field builds up to compensate for the external driving – forcing the electron in its dark state – until the maximal polarization is reached at the critical value . Above this point the nuclear system cannot compensate for the driving anymore and a solution of different nature, featuring finite electron inversion and entropy is found. The point features diverging spin entanglement and will be identified as a second order phase transition.

For the separable density matrix , the only term in Master Eq. (1) which is not trivially zero is the Hamiltonian term . However, choosing as an approximate eigenstate of the lowering operator (up to second order in ) with ( is the individual hyperfine coupling constant) the corresponding term in Eq. (1) vanishes in the thermodynamic limit. In Appendix A.1 we demonstrate that approximate eigenstates can be constructed as squeezed and displaced vacua in a Holstein-Primakoff Holstein and Primakoff (1940) picture up to a correction of order 1/J. The squeezing of the nuclear state depends uniquely on the displacement such that these states represent a subclass of squeezed coherent atomic states Kitagawa and Ueda (1993). Remarkably, this solution – where along the whole segment the system settles in a separable pure state – exists for all values of the dissipation strength .

In the limit of vanishing driving the steady state trivially is given by the fully polarized state (being the zero eigenstate of the lowering operator), as the model realizes a standard optical spin pumping setting for dynamical nuclear polarization Urbaszek et al. (2012). With increasing , the collective nuclear spin is rotated around the -axis on the surface of the Bloch sphere such that the effective Overhauser field in -direction compensates exactly for the external driving field on the electron spin. As a consequence along the whole segment the dissipation forces the electron in its dark state , and all electron observables, but also the entropy and some nuclear observables, are independent of .

Furthermore, the steady state displays increased nuclear spin squeezing in -direction (orthogonal to the mean polarization vector) when approaching the critical point. A common measure of squeezing is defined via the spin fluctuations orthogonal to the mean polarization of the spin system. A state of a spin- system is called spin squeezed Kitagawa and Ueda (1993), if there exists a direction , orthogonal to the mean spin polarization , such that:

(6) |

In Korbicz et al. (2005) it was shown that every squeezed state also contains entanglement among the individual constituents.
Moreover, if then the spin squeezed state contains
-particle entanglement Pezzé and Smerzi (2009); Hyllus et al. (2012).
In Appendix A.1 we show that the squeezing parameter in -direction for an approximate eigenstate is given as .
Note however, that this equation is valid only for . For higher squeezing the operator expectation values constituting the term of order can attain macroscopic values of order .
For we find that the nuclear spins are in a highly squeezed minimum uncertainty state, with -particle entanglement ^{3}

Since the lowering operator is bounded (), at where the nuclear field has reached its maximum value, the zero entropy solution constructed above ceases to exist. For large electron driving, where sets the dominant energy scale, the dissipation results in an undirected diffusion in the dressed state picture and in the limit the system’s steady state is fully mixed. In order to describe the system for driving strength , in Section IV.1 we develop a perturbative theory designed to efficiently describe a class of steady states where the electron and nuclear spins are largely decoupled and the nuclear system is found in a fully polarized and rotated state with, potentially squeezed, thermal Gaussian fluctuations (also referred to as rotated squeezed thermal spin states – RSTSS or the Gaussian mode) It is fully characterized by its mean polarization as well as the spin squeezing and effective temperature of the fluctuations (cf. Appendix B). Squeezed coherent atomic states, which constitute the solution along segment , appear as a limiting case of this class for zero temperature .

We conduct a systematic expansion of the system’s Liouville operator in orders of the system size , by approximating nuclear operators by their semiclassical values and incorporating bosonic fluctuations up to second order in an Holstein-Primakoff picture. The resulting separation of timescales between electron and nuclear dynamics is exploited in a formalized adiabatic elimination of the electron degrees of freedom. The semiclassical displacements (i.e., the electron and nuclear direction of polarization) are found self-consistently by imposing first order stability of the nuclear fluctuations. For a given set of semiclassical solutions we derive a second order reduced master equation for the nuclear fluctuations which, in the thermodynamic limit, contains all information on the nuclear state’s stability, its steady state quantum fluctuations and entanglement as well as the low excitation dynamics in the vicinity of the steady state and thus allows for a detailed classification of the different phases and transitions.

Using this formalism, we find that the system enters a new phase at the critical point , in which the nuclear field can no longer compensate for the external driving, leading to a finite electron inversion and a nuclear state of rising temperature for increasing driving strength. At the transition between the two phases, the properties of the steady state change non-analytically and in Section IV.2.2 we will find an asymptotic closing of the Liouvillian gap (cf. Section II) at the critical point, as the Liouvillian’s spectrum becomes continuous in the thermodynamic limit. We will characterize the critical point () as a second order phase transition.

Allowing for arbitrary hyperfine detunings , a phase boundary emerges from the second order critical point (line in Fig. 1), separating two distinct phases (blue) and (red) of the Gaussian mode. The subregion of indicates a region of bistability associated to the phase boundary and is discussed below.

At the semiclassical equations of motion feature two steady state solutions. Not only the trivial steady state of the spin pumping dynamics – the fully polarized state in direction – but also an inverted state where the nuclear system is fully polarized in direction is a (unstable) solution of the semiclassical system. Quantum fluctuations account for the decay of the latter solution of anomalous spin pumping behavior. The two semiclassical solutions (the corresponding quantum states are from here on referred to as the normal and anomalous spin pumping mode, respectively) persist for finite . As we show employing the formalism described above (Section IV.2.3), quantum fluctuations destabilize the mode of anomalous behavior in region of the phase diagram. The stable Gaussian solution in phase displays a behavior characterized by the competition of dissipation and the onsetting driving field . The nuclear state is highly polarized in the direction set by the decay, and the electron spin starts aligning with the increasing external driving field. Furthermore, the normal spin pumping mode of phase is characterized by a low effective spin temperature.

The analysis of the low excitation spectrum of the Liouvillian (Section IV.2.4) shows a direct vanishing of the ADR at the phase boundary between and , while the imaginary part of the spectrum is gapped at all times. At this boundary, the normal mode of phase destabilizes while at the same the metastable anomalous mode turns stable defining the second phase . The two mode energies are non-degenerate across the transition preventing a mixing of the two modes and the emergence of critical phenomena such as diverging entanglement in the system. Phase – anomalous spin pumping – is characterized by a large nuclear population inversion, as the nuclear field builds up in opposite direction of the dissipation. At the same time the electron spin counter aligns with the external driving field . In contrast to the normal mode of phase , phase features large fluctuations (i.e., high effective temperature) in the nuclear state, which increase for high , until at some point the perturbative description in terms of RSTSS breaks down and the system approaches the fully mixed state. Note that region also transforms continuously to via the lower two quadrants of the phase diagram (Fig. 1). In this supercritical region Clifford and Clifford (1999) no clear distinction of the two phases exist.

To complete the phase diagram, we employ numerical techniques in order to study steady state solutions that go beyond a RSTSS description in Section V. The subregion of labeled indicates a region of bistability where a second steady state solution (besides the normal spin pumping Gaussian solution described above) appears, featuring a non-Gaussian character with large fluctuations of order . Since this mode cannot be described by the perturbative formalism developed in Section IV (which by construction is only suited for low fluctuations ) we use numerical methods to study this mode in Section V for finite systems. We find that the non-Gaussian mode (in contrast to the Gaussian mode of region ) is polarized in direction and features large fluctuations of the order of . Additionally this solution displays large electron-nuclear connected correlations . It emerges from the anomalous spin pumping mode coming from region and the system shows hysteretic behavior in region closely related to the phenomenon of optical bistability Bowden and Sung (1979).

A fourth region is found in the lower half of the phase diagram (). In contrast to the previous regions, area has no effects on steady state properties. Instead the region is characterized by an anomalous behavior in the low excitation dynamics of the system. The elementary excitations in region are overdamped. Perturbing the system from its steady state, leads to a non-oscillating exponential return. This behavior is discussed at the end of Section IV.2.3, where we study the low excitation spectrum of the Liouvillian in this region within the perturbative approach.

In summary, all the phases and transitions of the system are displayed in Fig. 1. Across the whole phase diagram one solution can be described as a RSTSS – a largely factorized electron-nuclear state with rotated nuclear polarization and Gaussian fluctuations. Phase hereby represents a region of normal spin pumping behavior. The system is found in a cold Gaussian state, where the nuclear spins are highly polarized in the direction set by the electron dissipation and the electron spin aligns with the external driving for increasing field strength. In contrast, phase displays anomalous spin pumping behavior. The nuclear system displays population inversion (i.e., a polarization opposing the electron pumping direction) while the electron aligns in opposite direction of the driving field. Furthermore the state becomes increasingly noisy, quantified by a large effective temperature, which results in a fully mixed state in the limit of large driving strength . Along segment the state becomes pure and factorizes exactly with a nuclear field that cancels the external driving exactly. The nuclear state can be described using approximate eigenstates of the lowering operator which display diverging squeezing approaching the second order critical point . From this critical point a first order phase boundary emerges separating phases A and B. It is associated with a region of bistability (area ), where a second solution appears featuring a highly non-Gaussian character. The system shows hysteretic behavior in this region. Region is a phase characterized by its dynamical properties. The system shows an overdamping behavior approaching the steady state, which can be inferred from the excitation spectrum of the Liouvillian.

Let us now describe the phases and transitions involving the Gaussian mode in detail.

## Iv Perturbative Treatment of the Gaussian Mode

As seen in the previous section along the segment the system settles in a factorized electronic-nuclear state, where the nuclear system can be described as a lowering operator eigenstate up to second order in . Motivated by this result we develop a perturbative theory based on a self-consistent Holstein-Primakoff transformation that enables the description of a class of steady states, which generalizes the squeezed coherent atomic state solution along to finite thermal fluctuations (RSTSS, Appendix B). A solution of this nature can be found across the entire phase diagram and we show that this treatment becomes exact in the thermodynamic limit.

In Section IV.2 we discuss this Gaussian mode across the whole phase diagram. Steady state properties of the nuclear fluctuations derived from a reduced second order master equation provide deep insights in the nature of the various phases and transitions. Observed effects include criticality in both steady state and low excitation spectrum, spin squeezing and entanglement as well as altered spin pumping dynamics. Whenever feasible we compare the perturbative results with exact diagonalization techniques for finite systems and find excellent agreement even for systems of a few hundred spins only. First in Section IV.2.2 we apply the developed theory exemplarily along the segment , to obtain further insights in the associated transition at . In Section IV.2.3 we then give a detailed description of the different phases that emerge in the phase diagram due to the Gaussian mode. Thereafter in Section IV.2.4 we conduct a classification of the different transitions found in the phase diagram.

### iv.1 The Theory

In this section we develop the perturbative theory to derive an effective second order master equation for the nuclear system in the vicinity of the Gaussian steady state.

For realistic parameters, the Liouville operator of Eq. (1) does not feature an obvious hierarchy, that would allow for a perturbative treatment. In order to treat the electron-nuclear interaction as a perturbation, we first have to separate the macroscopic semiclassical part of the nuclear fields. To this end we conduct a self-consistent Holstein-Primakoff approximation describing nuclear fluctuations around the semiclassical state up to second order.

The (exact) Holstein-Primakoff transformation expresses the truncation of the collective nuclear spin operators to a total spin -subspace in terms of a bosonic mode ( denotes the respective annihilation operator):

(7) | ||||

In the following we introduce a macroscopic displacement () on this bosonic mode to account for a rotation of the mean polarization of the state, expand the operators of Eq. (7) and accordingly the Liouville operator of equation Eq. (1) in orders of . The resulting hierarchy in the Liouvillian allows for an perturbative treatment of the leading orders and adiabatic elimination of the electron degrees of freedom whose evolution is governed by the fastest timescale in the system. The displacement is self-consistently found by demanding first order stability of the solution. The second order of the new effective Liouvillian then provides complete information on second order stability, criticality and steady state properties in the thermodynamic limit.

The macroscopic displacement of the nuclear mode

(8) |

allows for an expansion of the nuclear operators [Eq. (7)] in orders of

(9) | ||||

where

(10) | ||||

(11) | ||||

(12) | ||||

and . Analogously, one finds

(13) | ||||

(14) | ||||

(15) | ||||

(16) |

This expansion is meaningful only if the fluctuations in the bosonic mode are smaller than . Under this condition, any nuclear state is thus fully determined by the state of the bosonic mode and its displacement .

The zeroth order superoperator acts only on the electron degrees of freedom. This separation of timescales between electron and nuclear degrees of freedom implies that for a given semiclassical nuclear field (defined by the displacement ) the electron settles to a quasi-steady state on a timescale shorter than the nuclear dynamics and can be eliminated adiabatically on a coarse grained timescale. In the following we determine the effective nuclear evolution in the submanifold of the electronic quasi-steady states of .

Let be the projector on the subspace of zero eigenvalues of , i.e., the zeroth order steady states, and . Since features a unique steady state, we find , where denotes the trace over the electronic subspace and . By definition it is . After a generalized Schrieffer-Wolff transformation Kessler (), we derive an effective Liouvillian within the zeroth order steady state subspace in orders of the perturbation

(20) | ||||

After tracing out the electron degrees of freedom the dynamics of the nuclear fluctuations are consequently governed by the reduced master equation

(21) |

The first order term in of Eq. (20) can be readily calculated

(22) |

where is an electronic operator

(23) | ||||

denotes the steady state expectation value according to , which depends on the system parameters and and on the semiclassical displacement via optical Bloch equations derived from as described below. Eq. (22) represents a driving of the nuclear fluctuations to leading order in the effective dynamics. Thus for the steady state to be stable to first order, we demand

(24) |

This equation defines self-consistently the semiclassical nuclear displacement in the steady state in dependence on the system parameters and .

The calculation of the second order term of Eq. (20) is more involved and presented in Appendix D.
We find the effective nuclear master equation to second order ^{4}

(25) | ||||

with

(26) | ||||

(27) | ||||

and

(28) | ||||

For a given set of system parameters the coefficients defining the nuclear dynamics [Eq. (26), Eq. (27) and Eq. (28)] depend only on the nuclear displacement . After choosing self-consistently to fulfill Eq. (24) in order to guarantee first order stability, Eq. (25) contains all information of the nuclear system within the Gaussian picture, such as second order stability as well as purity and squeezing of the nuclear steady state. Also it approximates the Liouville operator’s low excitation spectrum to leading order and thus contains information on criticality in the system. Eq. (25) therefore forms the basis for the subsequent discussion of the RSTSS mode and the corresponding phases and transitions in Section IV.

In order to calculate the coefficients of Eq. (28), we have to determine integrated electronic autocorrelation functions of the type and , where . The dynamics of single electron operator expectation values are governed by the optical Bloch equations derived from

(29) |

where and and

(30) |

where we defined and is given in Eq. (14). The steady state solutions can readily be evaluated

(31) | |||

(32) |

Defining the correlation matrix and , the Quantum Regression Theorem Lax (1963) yields the simple result:

(33) | ||||

(34) |

Finally the time integrated autocorrelation functions reduce to the simple expression

(35) | ||||

(36) |

These matrices straightforwardly define the coefficients of the effective master equation of the nuclear fluctuations Eq. (25). In Appendix D.1 we provide explicit formulas to calculate the relevant coefficients.

### iv.2 Phase Diagram of the Gaussian Mode

In this Section we use the theory developed above to study the RSTSS mode across the phase diagram. As outlined in the previous section we first determine self-consistently possible semi-classical displacements , which guarantee first order stability [Eq. (24)]. For each of these solutions we determine the effective master equation for the nuclear fluctuations Eq. (25), which in the thermodynamic limit contains all information on the steady state and the low excitation dynamics and we discuss properties like second order stability, criticality as well as purity and squeezing of the nuclear steady state. Using this information we provide a complete picture of the various phases and transitions involving the RSTSS solution.

#### Methods and General Features

In order to determine the semiclassical displacements which guarantee first order stability, we show in Appendix C that Eq. (24) is equivalent to the semiclassical steady state conditions. Due to a symmetry in the equation, the steady state displacements appear in pairs , . Any semiclassical displacement can be straightforwardly converted to the mean spin polarizations up to leading order in according to Eq. (10), Eq. (14), Eq. (31), and Eq. (32). In the thermodynamic limit the two sets of steady state expectation values extracted from and share the symmetry . In large parts of the phase diagram the solution () displays high nuclear polarization in the same (opposite) direction as the the electron spin pumping. We define the corresponding quantum states as the normal (anomalous) spin pumping mode.

The two solutions define two corresponding master equations of the nuclear fluctuations around the respective semiclassical expectation values according to Eq. (25). These master equations are subsequently used to determine second order stability of the nuclear fluctuations and, if the dynamics turn out to be stable, the steady state properties of the nuclear system. We emphasize that the effective Master Eq. (25) not only can be used to determine steady state properties, but also reproduces accurately the low excitation spectrum of the exact Liouvillian. It thus also describes the system dynamics in the vicinity of the steady state (increasingly accurate for large ).

From Eq. (25) one readily derives a dynamic equation for the first order bosonic moments

(37) |

with

(38) | ||||

(39) | ||||

(40) |

where all parameters are functions of the semiclassical displacements . This equation of motion – and thus the corresponding master equation itself – features a fixed point if the eigenvalues of the matrix have negative real part (). Due to the symmetry between and one finds that the eigenvalues of the two matrices corresponding to fulfill such that across the whole phase diagram only one solution is stable at a time and defines the corresponding phase in the phase diagram. Note however, that the unstable solution decays at a rate that is second order in . Preparing the system in this state consequently leads to slow dynamics, such that this solution exhibits metastability.

In the following we implicitly choose the stable for which the real parts of the eigenvalues of are negative and discard the unstable solution. Fig. 2 displays a selection of steady state expectation values in the thermodynamic limit across the phase diagram for the stable solution. Different expectation values illustrate the different nature of phase and and show distinct signatures of first and second order phase transitions which will be discussed in greater detail in Section IV.2.3 and IV.2.4.

The approximate steady state polarizations found in this way coincide with the exact values found via diagonalization techniques to an extraordinary degree ( relative deviation for J=150). Corrections to the perturbative solutions are of the order since the first order expectation values of the bosonic mode vanish by construction, since [Compare Eq. (9) and Eq. (13)]. In the thermodynamic limit the perturbative solution becomes exact.

The two eigenvalues of are typically of the form (except in region which will be discussed below) and define the complex energy of the mode.
In this case the matrix contains all information on the low excitation spectrum of the Liouvillian which is approximated by multiples of the mode energies within the perturbative treatment
^{5}

The ADR according to the perturbative descriptions based on Gaussian modes is displayed in Fig. 3. It is used to study the transitions involving the Gaussian mode in the thermodynamic limit. The ADR vanishes along a line indicating a phase boundary separating the normal and anomalous spin pumping phase, which is described in Section IV.2.4. Furthermore a non-analyticity of the ADR at a finite value defines region , which characterizes a dynamical phase and is explained in Section IV.2.3.

The dynamical matrix of the first order moments provides information on the stability of the semiclassical solutions, the criticality of the Liouvillian and the non-analyticities of region . In order to understand the character of the solutions in the different regions of the phase diagram we consider next the steady state covariance matrix (CM) of the bosonic system. For a quadratic evolution like the one of Eq. (25) the steady state covariance matrix contains all information on the state. We deduce the effective temperature and the squeezing of the nuclear spin system, which connects to criticality in the system.

For a one-mode system with vanishing displacements and [in the steady state of Eq. (25) this is always the case] the CM is defined as

(41) |

with the usual definitions and . Using Eq. (25) we straightforwardly calculate the steady state covariance matrix across the phase diagram. As , is symplectically diagonalizable, with

(42) |

where is orthogonal with . For a single mode, and are real numbers. While is a measure of the purity of the state [], the smallest eigenvalue of , determines the amount of squeezing in the system Kim et al. (2002). indicates squeezing in the bosonic mode. For , the covariance matrix Eq. (42) describes a thermal state of the bosonic mode and can be straightforwardly associated to a dimensionless effective temperature

(43) |

This definition is also meaningful for , since the squeezing operation is entropy-conserving. is also a measure for the entropy of the spin system, as to leading order it is connected to the bosonic mode via an unitary (i.e., entropy-conserving) transformation. The effective temperature of the different phases will be discussed below in Sections IV.2.2 & IV.2.3 [cf. Fig. 7].

We stress the point that all properties of the covariance matrix derived within the second order of the perturbative approach are independent of the system size . In particular, the amount of fluctuations (i.e., the purity) in the state does not depend on the particle number. In order to self-consistently justify the perturbative approach, has to be small with regard to . This implies that in the thermodynamic limit the perturbative results to second (i.e., leading) order become exact.

The inverse purity is displayed in Fig. 4 a). Except for for a small region around the Gaussian phase boundary the fluctuations are much smaller than , which justifies the validity of the perturbative approach and explains the excellent agreement with the exact diagonalization for this system size.

The squeezing in the auxiliary bosonic mode does not necessarily correspond to spin squeezing in the nuclear system. In order to deduce the spin squeezing in the nuclear system from the squeezing of the bosonic mode a transformation according to Eq. (11) and Eq. (15) is necessary. In Appendix A.1 we show that for Eq. (11) can be reformulated to connect the spin fluctuations to a squeezed and rescaled bosonic mode

(44) |

where is the squeezing operator and and .

Thus squeezing of the mode does in general not imply reduced spin fluctuations in a direction orthogonal to the mean spin polarization since the transformation between spin fluctuations and involves a squeezing operation itself and a scaling by a factor .

In general we thus have to apply a more involved squeezing criterion. In Korbicz et al. (2005) it was shown that for systems of spin-1/2 particles and for all directions the quantity

(45) |

signals entanglement if for some direction . Moreover,
indicates a generalized spin-squeezing of the state
^{6}

In the following we will use the quantity to investigate squeezing and bipartite entanglement in the nuclear system. In order to calculate we reconstruct the approximate nuclear operators according to Eq. (9) and Eq. (14) from the semiclassical displacement and evaluate the expectation values according to the steady state covariance matrix Eq. (41). Finally we maximize with regard to all possible directions to obtain . The results are discussed in Section IV.2.4. As discussed in more detail in the next Section, the fact that as on the line segment indicates a diverging entanglement length in the sense that -particle entanglement is present Hyllus et al. (2012).

#### A Second Order Phase Transition: The Segment

The segment at (Fig. 1) represents a very peculiar region in the phase diagram, where the solution below the critical point can be constructed analytically as seen in Section III.2. The electron and nuclear system decouple, resulting in a zero entropy product steady state. A nuclear polarization builds up to cancel the external driving up to the point of maximal Overhauser field (). At this point squeezing and entanglement in the system diverge, indicating a second order phase transition. In the following we exemplarily employ the formalism developed above along this line to obtain further insight about the criticality at . We calculate the analytical steady state solution as well as the effective master equation governing the nuclear fluctuation dynamics in its vicinity. We find that here the spectrum of the Liouvillian becomes continuous (implying a closing gap) and real. At the same time the creation operators of the elementary excitations from the steady state turn hermitian giving rise to diverging spin entanglement.

The first order stability condition Eq. (24) is fulfilled, if [compare Eq. (31) and Eq. (32)], which yields the possible semiclassical steady state displacements

(46) | ||||

corresponding to a normal (’’) and anomalous (’’) spin pumping mode, respectively.

Next, we explicitly calculate the second order corrective dynamics of the nuclear degrees of freedom for the normal mode. The vanishing of the effective driving forces the electron in its dark state – implying – and directly yields [Eq. (26) and Eq. (27)]. The remaining constants can be calculated as described above and introducing new bosonic operators (for the normal mode )

(47) |

with

(48a) | |||

(48b) |

one finds the effective evolution of the nuclear fluctuations given as

(49) | ||||

with

(50) | |||

(51) |

and fulfill boson commutation relations, since Eq. (47) defines a symplectic transformation (). The eigenvalues of the dynamical matrix associated to Eq. (49) are straightforwardly given as . The real part – representing the ADR of the system in thermodynamic limit (compare Fig. 5) – is always negative, indicating the stability of the normal spin pumping mode (). In an analogous calculation one shows that the semiclassical solution is not stable to second order since the eigenvalues of have a positive real part, i.e. the fluctuations diverge, violating the initial assumptions that the mode has to be lowly occupied.

Selected steady state expectation values derived from the stable displacement to leading order in (i.e., in the thermodynamic limit) are displayed in Fig. 6.

Already for we find excellent agreement between the perturbative and exact mean polarizations. The semiclassical nuclear field builds up to exactly cancel the external magnetic field forcing the electron in its dark state along and thus realizing the model of cooperative resonance fluorescence Carmichael (1980) even for weak dissipation [compare Eq. (5)]. This solution is only available if (defining segment ), i.e., up to the point where the nuclear field reaches its maximum. At this point the system enters a new phase of anomalous spin pumping (described below) and the steady state properties change abruptly.

Inserting solution in the coefficients of Master Eq. (49) yields

(52) | |||

(53) |

In the close vicinity below the critical point the real part of the gap in the Liouvillian’s spectrum closes as

(54) |

and the imaginary part as

(55) |

indicating criticality. Fig. 5 displays the ADR along in the thermodynamic limit (which is given on the segment by Eq. (52)) and for finite systems. It displays an avoided crossing at with a gap that vanishes in the thermodynamic limit. This closing of the gap coincides with diverging timescales in the system, which renders the model more susceptible to potential perturbing effects, a phenomenon well known in the context of criticality Bowden and Sung (1979).

In contrast to the general form Eq. (25), Eq. (49) contains only one Lindblad term and the dynamics drive the system into the vacuum of the squeezed mode . As the system approaches the critical value (i.e., ) the mode adopts more and more a like character and thus the squeezing of this mode’s vacuum increases. The (in general complicated) transformation between the squeezing of the bosonic mode and the spin operators (cf. Section IV.2.1) can readily be established along , since the operator is trivially related to the spin operators [cf. Eq. (11)]

(56) | ||||

The fluctuations in -direction, for example, are consequently given as