Dissipative phase transition in a central spin system

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 ﬁrst-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 nitrogen-vacancy centers coupled to lattice nuclear spins are brieﬂy discussed.


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 [1][2][3][4][5][6]).
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, ρ 0 (g), which depends on the system and environment parameters, g. As g is changed, a sudden change in the system properties may occur, giving rise to a so called dissipative phase transition (DPT) [7][8][9][10][11][12]. 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 [13]. Moreover they offer the intriguing possibility of observing critical effects non-destructively because of the constant intrinsic exchange between system and environment [14]. 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 [7,15]. There, a system of spins interacts with a thermal reservoir and is externally driven. Experimental [16] and theoretical studies [17][18][19][20] 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 [21][22][23]. 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 [23] 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 L. 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 [24,25]. 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, H(g), where g denotes a set of systems parameters (like magnetic fields, interactions strengths, etc). At zero temperature and for a fixed set of parameters, g, the system is described by a quantum state, ψ 0 (g), fulfilling [H(g) − E ψ0 (g)]|ψ 0 (g) = 0, where E ψ0 (g) is the ground state energy. As long as the Hamiltonian is gapped (i.e., the difference between E 0 (g) and the first excitation energy is finite), any small change in g will alter the physical properties related to the state |ψ 0 (g) smoothly and we remain in the same phase. However, if the first excitation gap ∆ = E ψ1 (g) − E ψ0 (g) closes at a given value of the parameters, g = g 0 , it may happen that the properties change abruptly, in which case a phase transition occurs.
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 equatioṅ ρ = L(g)ρ. The dynamics describing the system are contractive implying the existence of a steady state. This steady state ρ 0 (g) is a zero eigenvector to the Liouville superoperator L(g)ρ 0 (g) = 0. This way of thinking parallels that of quantum phase transitions, if one replaces [H(g) − E ψ0 (g)] → L(g). 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 ρ 0 (g) (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 L(g) closes. The relevant gap in this context is determined by the eigenvalue with largest real part different from zero (it can be shown that Re(λ) ≤ 0 for all eigenvalues of L [26]). The vanishing of the real part of this eigenvalue -from here on referred to as asymptotic decay rate (ADR) [27] -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 [28].
Besides the analogies described so far [cf. Table I], there are obvious differences, like the fact that in DTP ρ 0 (g) 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.

A. 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 spindown state. It can be introduced via standard optical pumping techniques. [29,30]. 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 [23]-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 equatioṅ where {·, ·} denotes the anticommutator and S α and I α = σ α i (α = +, −, z) denote electron and collective nuclear spin operators, respectively. JΩ is the Rabi frequency of the resonant external driving of the electron (in rotating wave approximation), while δω = ω − a/2 is the difference of hyperfine detuning ω and half the individual hyperfine coupling strength a. δω for instance can be tuned via a static magnetic fields in z direction. Note that H I + H SI = a S I + ωI z , 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 J [51] is introduced here for convenience and will be justified later. Potential detunings of the electron driving -corresponding to a term ∆S z in the Hamiltonian part of the master equation -can be neglected if ∆ Ja. In the limit of strong dissipation γ a the electron degrees of freedom can be eliminated and Eq. (1) reduces toσ where γ eff = a 2 γ , Ω eff = Ωa 2γ and σ is the reduced density matrix of the nuclear system. This is a generalization of the Dicke model of resonance fluorescence as discussed in [7,10,20].
Master Eq. (1) has been theoretically shown to display cooperative nuclear effects such as superradiance even for inhomogeneous electron nuclear coupling [31]. 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.

B. Phenomenological Description of the Phase Diagram
For a fixed dissipation rate γ = a [52] the different phases and transitions of the system are displayed schematically in Fig. 1 in dependence on the external driving Ω and the hyperfine detuning ω. We concentrate our studies on the quadrant Ω, ω > 0 in which all interesting features can be observed. In the following, we outline the key features of the phase diagram.
First we consider the system along the line segment x (ω = ω 0 , Ω ≤ Ω 0 ), where Ω 0 = ω 0 = a/2 define a critical driving strength and critical hyperfine detuning, respectively. Here H I 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 Ω 0 . 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 Ω 0 features diverging spin entanglement and will be iden- Figure 1: Schematic of the different phases and transitions of Master Eq. (1). In the two main phases of the system A (yellow) and B (blue) -which together cover the whole phase diagram -the system is found in a RSTSS (cf. text). While phase A is characterized by normal spin pumping behavior (large nuclear polarization in the direction of the dissipation) and a low effective temperature, phase B displays anomalous spin pumping behavior (large nuclear polarization in opposing direction to the dissipation) and high temperature. They are separated by the first order phase boundary b which is associated with a region of bistability C (framed by the boundary c). Here a second non-Gaussian solution appears, besides the normal spin pumping mode of A. The region of bistability C culminates in a second order phase transition at (ω0, Ω0). Below this critical point the system is supercritical and no clear distinction between phases A and B exists. In this region a dynamical phase D emerges, characterized by anomalous behavior in the approach to the steady state. For a detailed description of the different phases and transitions see Section III B.
tified as a second order phase transition.
For the separable density matrix ρ 0 = |ψ ψ|, |ψ = |↓ ⊗ |α the only term in Master Eq. (1) which is not trivially zero is the Hamiltonian term S + ( a 2 I − + JΩ). However, choosing |α as an approximate eigenstate of the lowering operator I − |α ≈ α |α (up to second order in = 1/ √ J) with α = −2JΩ/a ≡ −JΩ/Ω 0 (a 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 [32] 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 [33]. Remarkably, this solution -where along the whole segment x the system settles in a separable pure state -exists for all values of the dissipation strength γ.
In the limit of vanishing driving Ω = 0 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 [34]. With increasing Ω, the collective nuclear spin is rotated around the y-axis on the surface of the Bloch sphere such that the effective Overhauser field in x-direction compensates exactly for the external driving field on the electron spin. As a consequence along the whole segment x 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 y-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-J system is called spin squeezed [33], if there exists a direction n, orthogonal to the mean spin polarization I , such that: In [35] it was shown that every squeezed state also contains entanglement among the individual constituents. Moreover, if ξ 2 n < 1 k then the spin squeezed state contains k-particle entanglement [36,37]. In Appendix A 1 we show that the squeezing parameter in y-direction for an approximate I − eigenstate |α is given as For Ω Ω 0 we find that the nuclear spins are in a highly squeezed minimum uncertainty state, with k-particle entanglement [53]. Close to the critical point k becomes of the order of √ J [ξ 2 ey = O(1/ √ J)] indicating diverging entanglement in the system.
Since the lowering operator is bounded (||I − || ≤ J), at Ω = Ω 0 where the nuclear field has reached its maximum value, the zero entropy solution constructed above ceases to exist. For large electron driving, where Ω Ω 0 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 Ω > Ω 0 , in Section IV A 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 T eff of the fluctuations (cf. Appendix B).
Squeezed coherent atomic states, which constitute the solution along segment x, appear as a limiting case of this class for zero temperature T eff = 0.
We conduct a systematic expansion of the system's Liouville operator in orders of the system size 1/ √ J, 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 Ω 0 , 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 B 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 (ω 0 , Ω 0 ) as a second order phase transition.
Allowing for arbitrary hyperfine detunings ω, a phase boundary emerges from the second order critical point (line b in Fig. 1), separating two distinct phases A (blue) and B (red) of the Gaussian mode. The subregion C of A indicates a region of bistability associated to the phase boundary b and is discussed below.
At Ω = 0 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 −z direction -but also an inverted state where the nuclear system is fully polarized in +z 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 B 3), quantum fluctuations destabilize the mode of anomalous behavior in region A of the phase diagram. The stable Gaussian solution in phase A 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 A is characterized by a low effective spin temperature.
The analysis of the low excitation spectrum of the Liouvillian (Section IV B 4) shows a direct vanishing of the ADR at the phase boundary b between A and B, while the imaginary part of the spectrum is gapped at all times. At this boundary, the normal mode of phase A destabilizes while at the same the metastable anomalous mode turns stable defining the second phase B. 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 B -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 A, phase B 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 A also transforms continuously to B via the lower two quadrants of the phase diagram ( Fig. 1). In this supercritical region [38] 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 A labeled C 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 J. Since this mode cannot be described by the perturbative formalism developed in Section IV (which by construction is only suited for low fluctuations J) 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 A) is polarized in +z direction and features large fluctuations of the order of J. Additionally this solution displays large electron-nuclear connected correlations S i I j − S i I j . It emerges from the anomalous spin pumping mode coming from region B and the system shows hysteretic behavior in region C closely related to the phenomenon of optical bistability [39].
A fourth region is found in the lower half of the phase diagram (D). In contrast to the previous regions, area D 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 D 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 B 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 A 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 B 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 x 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 I − which display diverging squeezing approaching the second order critical point Ω 0 . From this critical point a first order phase boundary emerges separating phases A and B. It is associated with a region of bistability (area C), where a second solution appears featuring a highly non-Gaussian character. The system shows hysteretic behavior in this region. Region D 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 x 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 = J −1/2 . 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 x 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 B 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 B 2 we apply the developed theory exemplarily along the segment x, to obtain further insights in the associated transition at Ω 0 . In Section IV B 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 B 4 we conduct a classification of the different transitions found in the phase diagram.

A. 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 L 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 selfconsistent 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 J-subspace in terms of a bosonic mode (b denotes the respective annihilation operator): In the following we introduce a macroscopic displacement √ Jβ ∈ C (|β| ≤ 2) 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 = 1/ √ J. 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 allows for an expansion of the nuclear operators [Eq. (7)] in orders of . . . and k = 2 − |β| 2 . Analogously, one finds This expansion is meaningful only if the fluctuations in the bosonic mode b are smaller than √ J. Under this condition, any nuclear state is thus fully determined by the state of the bosonic mode b and its displacement β.
According to the above expansions Master Eq. (1) can be written aṡ where The zeroth order superoperator L 0 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 quasisteady states of L 0 .
Let P be the projector on the subspace of zero eigenvalues of L 0 , i.e., the zeroth order steady states, and Q = 1 − P . Since L 0 features a unique steady state, we find P ρ = Tr S (ρ) ⊗ ρ ss , where Tr S denotes the trace over the electronic subspace and L 0 ρ ss = 0. By definition it is P L 0 = L 0 P = 0. After a generalized Schrieffer-Wolff transformation [40], we derive an effective Liouvillian within the zeroth order steady state subspace in orders of the perturbation After tracing out the electron degrees of freedom the dynamics of the nuclear fluctuations b are consequently governed by the reduced master equatioṅ The first order term in of Eq. (20) can be readily calculated where A is an electronic operator A ss denotes the steady state expectation value according to L 0 , which depends on the system parameters γ and Ω and on the semiclassical displacement β via optical Bloch equations derived from L 0 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 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 [54] 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 The dynamics of single electron operator expectation values are governed by the optical Bloch equations derived from L 0 where ∆ S := S − S ss and S = (S + , S − , S z ) T and where we definedΩ = Ω + a 2 √ kβ and L z 0 is given in Eq. (14). The steady state solutions can readily be evaluated Defining the correlation matrix S = ∆ S∆ S † ss and S t = ∆ S t ∆ S † ss , the Quantum Regression Theorem [41] yields the simple result: Finally the time integrated autocorrelation functions reduce to the simple expression 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.

B. 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 (± S x , S y , S z , I x , ± I y , ± I z ). 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 J).
From Eq. (25) one readily derives a dynamic equation for the first order bosonic momentṡ 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 (Re[λ 1,2 ] < 0). Due to the symmetry between β + and β − one finds that the eigenvalues of the two Σ matrices corresponding to β ± fulfill Re[λ 1,2 (β + )] = −Re[λ 1,2 (β − )] 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  The two eigenvalues of Σ are typically of the form λ 1,2 = a ± ib (except in region D 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 [55]. The low excitation spectrum contains information about criticality of the system and the dynamics in the vicinity of the steady state, and will be used to discuss and classify the different transitions in the phase diagram. In particular the eigenvalue of Σ with largest real part approximates the ADR in the thermodynamic limit in those regions of the phase diagram where the Gaussian mode is responsible for the lowest excitations in the Liouvillian spectrum (only in the region of bistability C this is not the case).
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 b indicating a phase boundary separating the normal and anomalous spin pumping phase, which is described in Section IV B 4. Furthermore a non-analyticity of the ADR at a finite value defines region D, which character-izes a dynamical phase and is explained in Section IV B 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 D. 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 x and p [in the steady state of Eq. (25) this is always the case] the CM is defined as with the usual definitions x = 1 Using Eq. (25) we straightforwardly calculate the steady state covariance matrix Γ ss across the phase diagram. As Γ = Γ T > 0, Γ is symplectically diagonalizable, with where O is orthogonal with det(O) = 1. For a single mode, D ≥ 1 and M ≥ 1 are real numbers. While D is a measure of the purity of the state [T r(ρ 2 ) = 1/ |Γ| = 1/D], the smallest eigenvalue of Γ, λ min ≡ DM −2 determines the amount of squeezing in the system [42]. λ min < 1 indicates squeezing in the bosonic mode. For M = 1, the covariance matrix Eq. (42) describes a thermal state of the bosonic mode and D can be straightforwardly associated to a dimensionless effective temperature This definition is also meaningful for M > 1, since the squeezing operation is entropy-conserving. T eff 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 B 2 & IV B 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 J. 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, D has to be small with regard to J. This implies that in the thermodynamic limit J → ∞ the perturbative results to second (i.e., leading) order become exact. The inverse purity D is displayed in Fig. 4 a). Except for for a small region around the Gaussian phase boundary b the fluctuations are much smaller than J = 150, which justifies the validity of the perturbative approach and explains the excellent agreement with the exact diagonalization for this system size.
The squeezing λ min 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 |β| < 1 Eq. (11) can be reformulated to connect the spin fluctuations to a squeezed and rescaled bosonic mode where S(r) = e (r * b 2 −rb †2 )/2 is the squeezing operator and cosh(r) = µ = (2k−|β| 2 )/[2 2k(1 − |β| 2 )] and sinh(r) = Thus squeezing λ min of the mode b does in general not imply reduced spin fluctuations in a direction orthogonal to the mean spin polarization since the transformation between spin fluctuations and b involves a squeezing operation itself and a scaling by a factor 0 < 2(1 − |β| 2 ) ≤ √ 2.
In general we thus have to apply a more involved squeezing criterion. In [35] it was shown that for systems of N spin-1/2 particles and for all directions n the

quantity
signals entanglement if C n > 0 for some direction n. Moreover, ∆I 2 n < J/2 indicates a generalized spinsqueezing of the state [56].
In the following we will use the quantity C = max{0, C n | n ∈ R 3 } to investigate squeezing and bipartite entanglement in the nuclear system. In order to calculate C n 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 C n with regard to all possible directions n to obtain C. The results are discussed in Section IV B 4. As discussed in more detail in the next Section, the fact that C → 1 as Ω → Ω 0 on the line segment x indicates a diverging entanglement length in the sense that O(1/(1 − C)) = O( √ J)-particle entanglement is present [37].

A Second Order Phase Transition: The Segment x
The segment x at ω = ω 0 (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 B. 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 (Ω 0 ). 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 (ω 0 , Ω 0 ). 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.
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 Ω = 0 forces the electron in its dark state -implying S + ss = S − ss = S + S − ss = 0 -and directly yields B = F = 0 [Eq. (26) and Eq. (27)]. The remaining constants can be calculated as described above and introducing new bosonic operators (for the normal mode β = β − ≤ 1) one finds the effective evolution of the nuclear fluctuations given asσ d and d † fulfill boson commutation relations, since Eq. (47) defines a symplectic transformation (|µ| 2 −|ν| 2 = 1). The eigenvalues of the dynamical matrix Σ associated to Eq. (49) are straightforwardly given as λ 1,2 = −Γ eff /2 ± iΘ eff . 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 β + > 1 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 b has to be lowly occupied. Selected steady state expectation values derived from the stable displacement β − to leading order in J (i.e., in the thermodynamic limit) are displayed in Fig. 6. Already for J = 150 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 x and thus realizing the model of cooperative resonance fluorescence [7] even for weak dissipation γ ≤ a [compare Eq. (5)]. This solution is only available if Ω ≤ Ω 0 (defining segment x), 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
In the close vicinity below the critical point Ω 0 the real part of the gap in the Liouvillian's spectrum closes as and the imaginary part as indicating criticality. Fig. 5 displays the ADR along ω = ω 0 in the thermodynamic limit (which is given on the segment x by Eq. (52)) and for finite systems. It displays an avoided crossing at Ω 0 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 [39]. In contrast to the general form Eq. (25), Eq. (49) contains only one Lindblad term and the dynamics drive the system into the vacuum |0 d of the squeezed mode d. As the system approaches the critical value Ω = Ω 0 (i.e., β − = −1) the mode d adopts more and more â p = 1 √ 2i (b − b † ) like character and thus the squeezing of this mode's vacuum increases. The (in general complicated) transformation between the squeezing of the bosonic mode b and the spin operators (cf. Section IV B 1) can readily be established along x, since the operator d is trivially related to the spin operators [cf. Eq. (11)] The fluctuations in y-direction, for example, are consequently given as . One readily shows that up to order O(1) and we used d = 0 in the steady state. In thep-vacuum |0 p it is p 2 d = 1/2, such that we evaluate ξ 2 ey = 2 ∆I 2 y /| I | (59) where we used | I | = J and inserted the semiclassical displacement β − . This is the same result we derived in Section III B and Appendix A 1 by constructing approximate eigenstates of the lowering operator I − and along x we find that C ≈ 1 − ξ 2 ey , as shown in Fig. 4 d). Note that hereê y is orthogonal to the direction of the mean spin I . This allows us to deduce that O( √ J) nuclear spins must be entangled close to the critical point, which establishes a 'diverging entanglement length' in this system. To see this, we employ a variant of the criterion Eq. (6) as discussed in [36]. There, it was shown that ξ 2 ey < 1/k sets a lower bound of N ξ −2 ey on the quantum Fisher information F Q of the state. In [37] it was shown that for states containing at most k-particle entanglement, F Q is upper bounded by N k. Consequently, the values of ξ 2 ey obtained close to the critical point (cf. Eq. (59) and Appendix A 1) imply that at least O( √ J)-particle entanglement must be present. Note that the bosonic description does not allow to describe the range ξ 2 ey = O(1/J), i.e., k = O(J), where the fluctuations become larger than the expansion parameter.
The nuclear squeezing and entanglement in the system diverges approaching the critical point, as the Lindblad operator d (defining the steady state |0 d ) becomes more and morep-like. The fluctuations in y-direction tend to zero, while at the same time -due to the Heisenberg uncertainty relation -the steady state is in a superposition of an increasing number of I z eigenstates. Since in a system with infinite range interactions (as the one we are considering) there is no obvious definition of a coherence length, the range of the involved I z eigenstates can be considered as an analogous concept.
At the critical value Ω = Ω 0 the symplectic transformation Eqs. (47) becomes ill defined (d becomes ap-like operator) while both the dissipation rate and the mode energy tend to zero. While the coefficients in Eqs. (48) diverge, the total master equation is well defined [due to the factors (1 − |β| 2 ) in Γ eff ] and straightforwardly can be written asσ The Liouville operator's spectrum is real and continuous with hermitian creation operators of the elementary excitations.
We stress the point that along segment x in the phase diagram highly dissipative dynamics drive the system in a pure and separable steady state with zero effective temperature T eff = 0 [cf. Fig. 7 b)]. At the critical point Ω 0 the steady state changes its nature abruptly as the system enters a high temperature phase.
Furthermore we remark that this steady state has no relation to the system's ground state. This is in contrast to the extensively studied Dicke phase transition [13,43,44] where the steady state is in close relation to the Hamiltonian's ground state (in fact in the normal phase it is identical). In the present model dissipation drives the system to a highly excited state of the Hamiltonian and the observed critical phenomena are disconnected from the Hamiltonian's low excitation spectrum.
We have seen that at the critical point (ω 0 , Ω 0 ) the gap of the Liouville operator's spectrum (in both real and imaginary part) closes in the thermodynamic limit [Eq. (54) and Eq. (55)]. Approaching the critical point the steady state fluctuations become more and more squeezed due to the increasingp-like character of the mode d. The spin squeezing close to the critical point [Eq. (59)] can be interpreted as a diverging coherence length in a system with infinite range interactions (the electron mediates interactions between remote spins). These are clear indications for a second order phase transition, which will be formalized in Section IV B 4.
In the present Section we study the different phases of the system, which involve the RSTSS solution (A, B and D) using the analytic tools developed above. By construction, the RSTSS solution describes steady states where the electron and nuclear state factorize to leading order in the system size and the nuclear system is found in a fully polarized and rotated state with Gaussian fluctuations, which are fully characterized by their effective temperature and squeezing. Fig. 2 displays different steady state observables of the Gaussian solution determined via the formalism described above to leading order in the thermodynamic limit. Temperatures T eff > 6 are cut off, as the temperature diverges along the phase boundary b. a) The first order phase boundary b separates the low temperature phase A from the high temperature phase B. b) T eff along ω = ω0: On segment x the system is in a zero entropy state (T eff =0). Above the second order critical point Ω > Ω0 the system enters a high temperature phase. Here the temperature rises with increasing driving strength.
In phase A the system is characterized by normal spin pumping behavior. Only the semiclassical displacement β − (normal mode) leads to a dynamical matrix Σ that has negative real parts of its eigenvalues, while for β + the eigenvalues have positive real parts, indicating the instability of that mode in second order. The nuclear system in the normal mode settles in a state highly polarized in -z direction following the direction of the electron spin pumping [ Fig. 2(a)]. Meanwhile, increasing the external driving Ω and approaching the phase boundary b, a nuclear field in x direction builds up, but only along x it can fully cancel the external driving [ Fig. 2(c)]. Therefore, in general the electron spin aligns more and more with the external field [ Fig. 2(b,d)]. Furthermore, the effective temperature (and thus the entropy) of the phase is low, as displayed in Fig. 7 a).
In region B in contrast, β + is the only stable solution, defining the phase of anomalous spin pumping behavior. The nuclear system now shows strong population inversion, i.e., the nuclear polarization is in direction opposite to the external pumping (z). In the same way the electron now aligns in opposite direction to the external driving field (x). Also, in contrast to phase A, the RSTSS now is in a high-temperature state. For larger electron driving the temperature increases until eventually the Gaussian description breaks down (as D ∝ J) and for Ω → ∞ the system is found in a completely mixed state [compare Fig. 4 b)].
In the upper half of the phase diagram (ω > ω 0 ) phase A changes abruptly into phase B at the boundary b and certain steady state spin observables [ I z , S x [Fig. 2 a) & b)] and I y (not displayed)] show distinct features of a first order phase transition, changing sign as the normal (anomalous) mode destabilizes (stabilizes). This transition is discussed in greater detail in the following Section IV B 4. Following this boundary towards the critical point (ω 0 , Ω 0 ) the two phases become progressively more similar. Below the critical point (ω < ω 0 ) there is no clear distinction between the normal and anomalous spin pumping mode anymore, a phenomenon known from thermodynamics as supercriticality. Phase A transforms continuously to phase B in this region. Close to the critical point, supercritical media typically respond very sensitive to the external control parameters of the phase diagram (e.g., temperature or pressure) [38]. In our system we observe that small changes in the parameter ω leads to large changes in electron spin observables.
Next, we consider the third region associated with the RSTSS solution, region D. We will find that this region differs from the previous ones by the fact that it cannot be detected in the system's steady state but rather in dynamical observables.
The eigenvalues of the dynamic matrix Σ can be calculated as λ 1,2 = −(R a − R b ) ± 2 4|ξ| 2 − χ 2 and provide information on the approximate low excitation spectrum of the Liouvillian. We can distinguish two cases for the low excitation spectrum, which differ only in the Hamiltonian properties of Eq. (25) (fully determined by χ and ξ [Eq. (39) & Eq. (40)]). In the first case the quadratic bosonic Hamiltonian can be symplectically transformed to be diagonal in a Fock basis (i.e., of the form ∝b †b ). This is the case if χ 2 > 4|ξ| 2 . As a consequence the two eigenvalues of Σ have an identical real part and imaginary parts ±2 χ 2 − 4|ξ| 2 . In the second case the Hamiltonian transforms symplectically into a squeezing Hamiltonian ∝ (b †2 +b 2 ). Here one finds χ 2 < 4|ξ| 2 , such that the eigenvalues become real and symmetrically distributed around −(R a − R b ). In region D in Fig. 1 we find the effective Hamiltonian for the nuclear fluctuations to be symplectically equivalent to a squeezing Hamiltonian. Fig. 8 shows the ADR exemplarily along the line ω = 0.5 ω 0 ( l II in Fig. 1) calculated according to the perturbative theory and via exact diagonalization, respectively. The perturbative theory approximates accurately the low excitation spectrum of the Liouvillian. We find that in region D the ADR splits up, when the coherent part of Eq. (25) changes to a squeezing Hamiltonian. As mentioned above this non-analyticity occurs at a non-zero value of the ADR and thus does not leave signatures in the steady state behavior. The steady state transforms smoothly along l II . However the nature of dynamical observables change within region D as the system displays anomalous behavior approaching the steady state. The splitting of the ADR coincides with the vanishing of the imaginary part of the lowest non-zero Liouvillian eigenvalues. Thus the system is overdamped in D. Perturbing the system from its steady state will not lead to a damped oscillatory behavior, but to an exponential, oscillation-free return to the steady state. The blue area in vicinity to region D in Fig. 3 does not represents a new phase but is another interesting feature of the system. Here, the ADR exceeds the value at Ω = 0 by a factor ∼ 3. For Ω = 0 the model describes the standard spin pumping setting. Large gaps in the low excitation spectrum indicate the possibility to improve the effective spin pumping rate (remember that also in this region the steady state is fully polarized, however not in −z direction as it is the case for the normal spin pumping configuration Ω = 0). Indeed simulations show that starting from a fully mixed state, the system reaches the steady state faster than in the standard setting (Ω = 0). This feature becomes more distinct in systems, where the electron pumping rate γ is limited. For γ = 0.1a the time to reach the fully polarized steady state from a fully mixed state is shortened by a factor ∼ 6.

Transitions
In this Section we consider the transitions involving the RSTSS solution in greater detail providing a classification in analogy to quantum phase transitions in closed systems (compare Section II).
As seen in the previous Section, certain steady state observables show clear signatures of a first order phase transition at b (Fig. 2). In order to understand this sharp transition we consider the ADR exemplarily along path l I in Fig. 9. The broken lines represent numeric results   Fig. 3. The vertical black lines indicate the asymptotic boundaries of the region of bistability. In the whole region the ADR tends to zero in the thermodynamic limit due to the appearance of a non-gaussian stable mode. Inset: The next higher excitations in the spectrum for J = 150 display equidistant splittings in regions far from the region of bistability. This is an indication for the bosonic character of the steady state, which is exploited in the perturbative approach.
of exact diagonalization of the Liouvillian for J = 50, 100 and 150, while the solid line indicates the result of the perturbative approach. As described in Section IV B 1 we implicitly choose the semiclassical displacement β − (for Ω < 1.5Ω 0 ) or β + (for Ω > 1.5Ω 0 ) for which the ADR is negative, indicating a stable solution. For increasing system size the ADR is increasingly well approximated by the perturbative solution.
We stress the point that the red line represents the first Gaussian excitation energy only. However, within the region of bistability (indicated by two vertical bars and discussed below in Section V), a non-Gaussian mode (discussed below) is responsible for additional excitations in the exact spectrum. The Gaussian mode eigenvalue (red line) in this region is reproduced approximately by higher excitations of the exact spectrum (not displayed) . The perturbative theory is still correct within the region of bistability but, as expected, it misses all non-Gaussian eigenstates of the exact Liouvillian.
At the boundary b ( Ω ≈ 1.5 Ω 0 ) the gap in the real part of the spectrum of the Liouvillian closes nonanalytically, indicating critical behavior. This observation is supported by the effective temperature (and thus the fluctuations in the system), which is increased in the vicinity of the boundary b, and diverges at the boundary [ Fig. 7 a) & Fig. 4 a)]. The vanishing of the ADR at b (i.e., the vanishing due to the RSTSS solution) can be observed at finite J (dashed lines in Fig. 9) and is not a feature appearing in the thermodynamic limit only. The position of this closing of the gap -which in the thermodynamic limit (solid line) is found at Ω ≈ 1.5 Ω 0 -is shifted for finite system sizes to lower drivings Ω.
The origin of this closing of the Liouvillian gap becomes more transparent if we take the mode energy of the respective metastable solution into account.
In Fig. 10 a)  Together with the discontinuous change in system observables such as mean polarizations we classify this Gaussian transition as of first order. Second, we consider the transition along ω = ω 0 (in-cluding the line segment x). In contrast to the situation before we find that the semiclassical displacements β + and β − merge approaching the critical point such that the two modes become asymptotically identical at Ω 0 [Eq. (46)]. Approaching the critical point, the eigenvalues of the two modes tend to zero (both the real and imaginary parts), causing the gap of the Liouvillian's spectrum to close [ Fig. 10 b), Eq. (54), Eq. (55)]. As we have seen in Section IV B 2 at (ω 0 , Ω 0 ) the spectrum becomes real and continuous signaling criticality. The perturbative treatment intrinsically is a description in the thermodynamic limit. If we consider the exact spectrum we indeed find an avoided crossing due to the mode mixing at the critical point with a gap that is closing for J → ∞ (cf. Fig. 5). As we discussed in Section IV B 2 the elementary excitations becomep-like, causing a diverging coherence length in the system [indicated by the diverging squeezing parameter C in Fig. 4(c,d)]. Together with the continuous but non-analytical change of the mean polarizations these properties classify the point (Ω 0 , ω 0 ) as a second order transition.

V. REGION OF BISTABILITY: NON-GAUSSIAN SOLUTION
As noted in Section III B along the Gaussian boundary b extends a region of bistability [C in (Fig. 1)] -culminating in the critical point (Ω 0 , ω 0 ) -in which a second stable solution appears. Within the perturbative framework from Section IV this highly non-Gaussian solution could not be detected because it features large fluctuations of the order of the system size J. In the following we use numerical techniques to construct and study this mode for finite systems. In the thermodynamic limit the ADR tends to zero within C, such that there exists a two dimensional subspace of steady states. Here we find two independent, physical solutions within the two dimensional kernel of the Liouvillian, one of which will turn out to be the Gaussian normal spin pumping mode described in Section IV. We analyze the nature and properties of the other, non-Gaussian solution, exemplarily along the line ω = 1.5 ω 0 ( l I in Fig. 1). Fig. 9 displays the ADR for different particle numbers. Within the indicated region of bistability (the black vertical lines represent the boundaries c and b, respectively) the ADR tends to zero with increasing particle number. Already for J = 150 one finds a small region, where the ADR is small enough (of the order of 10 −6 a) that one can construct two linearly independent (quasi) steady state solutions. Although we find the eigenmatrix ρ 1 associated with the ADR to be non-positive and traceless (the latter being a consequence of L being the generator of a trace preserving map) we can linearly combine it with the true steady state ρ 0 to obtain two linear independent, positive solutions with trace one, ρ lo (corresponding to the normal spin pumping mode) and ρ up . These solutions span the two dimensional space of steady states in that region.   11 illustrates the solutions ρ lo and ρ up around the bistable region in an equally weighted mixture. The density matrices are represented by their diagonal elements in the I z basis. In the plane the blue dots (red diamonds) represent the polarization in z-direction I z of the lower (upper) solution ρ lo (ρ up ). Coming from below the critical region (Ω < 1.15 Ω 0 ) the nuclear system is found in the Gaussian normal spin pumping mode, fully polarized, slightly rotated away from the −z direction and with fluctuations of the order of √ J. This Gaussian solution persists within the critical region where it becomes noisier until eventually -approaching the right boundary b at Ω = 1.5 Ω 0 -it destabilizes. In the thermodynamic limit the lower solution is stable up to the right boundary, where a first order transition occurs and the anomalous spin pumping mode appears. Approaching boundary b from above (Ω > 1.5 Ω 0 ) this mode transforms into a non-Gaussian solution, which -in contrast to the coexisting normal mode -features fluctuations of the order of J, is not fully polarized and shows large electron-nuclear and nuclear-nuclear connected correlations. Approaching the left boundary c at Ω = 1.15 Ω 0 this mode destabilizes eventually as the ADR becomes finite again and the normal mode is the only stable solution in the system.
The bistable behavior of the system in region C bears close resemblance to the phenomenon of optical bistability for saturable absorbers [45], where connections to phase transitions have been established [39]. In this region the system displays strong hysteretic behavior. Recent experiments in quantum dots, realizing a setting close to our model system display distinct signatures of hysteresis upon application of an external driving field on the electronic spin [24,25]. Our results suggest the observed optical bistability in central spin systems as a possible pathway to understand these experimental results, which will be a subject of further studies.

VI. IMPLEMENTATIONS AND EXTENSIONS OF THE MODEL
In the present section we discuss potential physical realizations of the Master Eq. (1) and address certain aspects of an extension of the model for inhomogeneous hyperfine couplings.
As mentioned above the model we study is a generic central spin model with various potential physical implementations. The most prominent ones represent singly charged semiconductor quantum dots, where the electron spin couples to the nuclear spins of the host material [23,34], and diamond nitrogen vacancy centers coupled to either nuclear ( 13 C spins of the host material) or electron (e.g., nearby nitrogen impurities) spin ensembles [46,47]. Recently diamond nano-crystals containing single NV centers coated with organic molecule spin labels, which are dipole coupled to the NV center spin have been manufactured [48].
NV centers represent a natural realization of the Master Eq. (1). Their ground state consists of three spin sublevels (of spin projection quantum number m = 0, ±1) featuring a zero field splitting due to anisotropic crystal fields of 2.88 GHz [46]. In a static magnetic field this zero field splitting can be compensated for and one of the transitions (e.g., m = 0 ↔ 1) is brought into near hyperfine resonance with the ancilla spin system, defining an effective two-level system. Since the m = 0 level does not carry a magnetic moment, the hyperfine interaction of the effective two level system and the ancilla system takes the anisotropic form of Eq. (4). Potential counterrotating terms of the dipole-dipole interaction are neglected in the static magnetic field in a rotating wave approximation. Optical pumping of the electron spin in the m = 0 spin state and resonant driving (either by optical Raman transitions or radio frequency fields) realizes Master Eq. (1) [30].
In general the hyperfine interaction in such a setting will not be homogeneous and the truncation to a symmetric subspace of total spin J is not justified. In the following we consider an extension of the model taking into account the inhomogeneous nature of the hyperfine coupling in a shell model. Along x we show that up to the critical point steady states can be constructed analytically as factorized product states involving nuclear eigenstates of the (inhomogeneous) lowering operator. In analogy to the homogeneous case, such solutions cease to exist after the critical point at which we find diverging nuclear squeezing. These results are supported by numerical simulations that confirm the analytical considerations and provide further indications that other features of the phase diagram aside from the second order transi-tion can be found in the inhomogeneous model.
In order to take into account inhomogeneities in the hyperfine coupling, we replace the homogeneous spin operators of Eq. (4) with inhomogeneous operators I α → A α (α = x, y, z). We approximate the actual distribution of coupling strengths by n shells of spins with identical coupling where A (i) α represent homogeneous spin operators within the ith shell. Each homogeneous shell is assumed to be in a symmetric subspace J i .
In analogy to the homogeneous case we can construct approximate eigenstates of the lowering operator A − |α = α |α . To this end we perform a Holstein-Primakoff transformation on the homogeneous spin operators within each shell and displace the respective bosonic mode b i by β i and expand the resulting operators in orders of 1/ √ J i . As we demonstrate in Appendix A 2 the choice of a particular displacement β i uniquely defines the squeezing of the respective mode b i if we demand that the corresponding state is an A − eigenstate to second order in the expansion parameters, i.e., of order O( i 1/J i ). The corresponding eigenvalue is then given as . As discussed in Section III B, |ψ = |↓ ⊗ |α is a steady state of the evolution to second order, if α = i g i √ k i β i = −JΩ/Ω 0 . In contrast to the homogeneous case (n = 1) the latter condition does not determine the steady state uniquely. Several sets of displacements within the different shells can fulfill the steady state condition. However, all these microscopic realizations lead to the same macroscopic behavior of the system such as the locking of the electron inversion S z = 0. Furthermore, at the critical point, the solution is unique again (β i = 1 for all shells) and the considerations on entanglement of Appendix A 1 can be straightforwardly generalized to the inhomogeneous case with the result that also here at the critical point the entanglement in the system diverges indicating a second order phase transition. Obviously, above the critical point no such solution can be constructed and the system observables change non-analytically. Fig. 6 shows numerical results which confirm the above considerations. We find numerically the exact steady state solution for a model of two inhomogeneously coupled shells (g 1 = 2g 2 ) of size J 1,2 = 8 (broken lines), as well as for a system of 5 nuclear spins with coupling strengths ({g i } i=1...5 = {0.67, 0.79, 0.94, 1.15, 1.4}, dotted lines). For low driving strengths Ω we find the Overhauser field building up linearly, as expected. The emergence of the thermodynamic phase transitions can be anticipated already for these low particle numbers.
These analytical and numerical arguments for the emergence of a second order phase transition in the inhomogeneous case, suggest the possibility to find other features of the homogeneous phase diagram also in inhomogeneous systems, such as NV centers in diamond.
Another attractive realization of a central spin system is provided by singly charged semiconductor quantum dots: up to several 10 4 nuclear spins are coupled to a central spin-1/2 electron, driving and spin pumping of the electronic state have been demonstrated experimentally with high efficiency [29,49]. In this setting, however, the inhomogeneity of the hyperfine coupling and the absence of an m = 0 central spin state lead to a situation in which the effective nuclear Zeeman term H I in Eq. (1) becomes inhomogeneous (it is composed of Knight field, nuclear Zeeman energy, and the (homogeneous) detuning) and does not vanish for any choice of parameters. Therefore the above argument for a persistence of the second order phase transition does not apply. However, critical phenomena similar to the ones described above were observed in optically driven quantum dots [24]. The adaptation of our model this and other more general settings is subject to future studies.

VII. CONCLUSIONS
In analogy to closed systems where critical phenomena arise from non-analyticities of the Hamiltonian low energy spectrum, in open systems critical phenomena are intimately related to the low excitation spectrum of the Liouville operator. We investigated a generic driven and damped central spin model and its rich steady state behavior, including critical effects such as bistabilities, first and second order phase transitions and altered spin pumping dynamics. We developed a two-step perturbative theory involving the expansion of nuclear fluctuations up to second order in a self-consistent Holstein-Primakoff transformation and the subsequent adiabatic elimination of the electron degrees of freedom in the vicinity to the steady state, which enabled us to provide a complete picture of the system's phase diagram. Linking common ideas from closed system phase transitions to the dissipative scenario, we were able to introduce a classification of the different transitions in the phase diagram.
The relevance of the considered model involves two aspects. On the one hand Eq. (1) describes a simple yet rich model, which displays a large variety of critical phenomena. The limitation to symmetric states allows for an efficient (and in the thermodynamic limit exact) perturbative treatment that gives deep insights into the nature of dissipative critical phenomena from a fundamental point of view. On the other hand the central spin model is general enough to have realizations in a large variety of physical systems (e.g., quantum dots, Nitrogen-Vacancy centers). Our understanding of the critical phenomena in this model could provide insight into recent observation of critical behavior in such systems [24,25]. Furthermore the main features of the phase diagram discussed above can also be found if the central (two-level) spin is replaced by a different physical system, e.g., a larger spin or a bosonic mode. The theory developed in Section IV can straightforwardly be adapted to different scenarios and opens the possibility to study dissipative critical effects in a variety of different physical systems [13].
Finally, we showed that in a more realistic adaptation of the model incorporating an inhomogeneous hyperfine coupling, the second order phase transition persists, indicating the possibility that the phase diagram remains qualitatively correct in this experimentally more realistic case. A more thorough analysis of the effects of inhomogeneities is subject to future work.
where in the last equation we used the expansion Eq. (14). Note that the Bloch vector is orthogonal (up to order O(1)) to the y-direction for all (real) α and of length | I | = I x 2 + I y 2 + I z 2 = J + O(1). Using Eq. (A6) and the angular momentum commutation relations one readily calculates where as usual ∆O 2 := O 2 − O 2 and we used the identity 1 − ( √ kβ) 2 = (1 − β 2 ) 2 . Thus, we find for the squeezing parameter in ydirection, The squeezing diverges for the state that realizes the maximal eigenvalue of the lowering operator ( √ kβ = 1). This corresponds to a state fully polarized in x direction. and in the usual way we define  The first term of the second order of Eq. (20) is of the same form as the first order and can readily be calculated: T r S (P L 2 P ρ) = −i[a/2( S + ss J − 2 + S − ss J + 2 ) (D1) + (a S + S − ss + δω)J z 2 , σ], dτ T r s (P L 1 P L 1 P ρ) (D6) where α, β = †, void and the Einstein sum convention is used.

Calculation of the coefficients
In order to determine the coefficients Eq. (28) we have to calculate terms of the kind ∞ 0 dτ ∆A α (τ )∆A β (0) ss