Quantum fluctuations and correlations in open quantum Dicke models

In the vicinity of ground-state phase transitions quantum correlations can display non-analytic behavior and critical scaling. This signature of emergent collective effects has been widely investigated within a broad range of equilibrium settings. However, under nonequilibrium conditions, as found in open quantum many-body systems, characterizing quantum correlations near phase transitions is challenging. Moreover, the impact of local and collective dissipative processes on quantum correlations is not broadly understood. This is, however, indispensable for the exploitation of quantum effects in technological applications, such as sensing and metrology. Here we consider as a paradigmatic setting the superradiant phase transition of the open quantum Dicke model and characterize quantum and classical correlations across the phase diagram. We develop an approach to quantum fluctuations which allows us to show that local dissipation, which cannot be treated within the commonly employed Holstein-Primakoff approximation, rather unexpectedly leads to an enhancement of collective quantum correlations, and to the emergence of a nonequilibrium superradiant phase in which the bosonic and spin degrees of freedom of the Dicke model are entangled.

Spin-boson models are paradigmatic theoretical models describing, for instance, the coupling of matter with electromagnetic fields or vibrational modes. A prominent example is the so-called Dicke model [1], which provides a simple framework for the study of the interaction between a large ensemble of atoms, described by N spin-1/2 (two-level) particles, and an electromagnetic cavity field, described by a bosonic mode [cf. Fig. 1(a)]. This model has been thoroughly investigated in equilibrium [2][3][4][5][6][7], where it displays a second-order ground-state transition from a normal to a superradiant phase. While the order-parameter behavior is captured by a mean-field treatment [8], studying quantum correlations requires the analysis of quantum fluctuations [9][10][11][12][13]. In equilibrium, this is typically done within the so-called Holstein-Primakoff approximation [14]. This exploits that the system Hamiltonian can be written in terms of collective (macroscopic) spin operators, which approximately behave as bosons when the system is close to its groundstate.
Nowadays, also due to a debate concerning a no-go theorem on the experimental realization of the equilibrium Dicke model [15][16][17], the focus is on the investigation of Dicke physics in an open quantum system setting [cf. Fig. 1(a)]. Open Dicke models feature a nonequilibrium superradiant phase transition, see Fig. 1(b), which is exactly captured by a mean-field approach [18,19]. However, in these settings, analyzing quantum fluctuations is challenging [8,[20][21][22]. As a consequence, very little is known about correlations in the nonequilibrium Dicke model phase transition, and even less in the presence of local dissipative processes, such as local spindecay [cf. Fig. 1(a)].
In this paper, we provide a complete characterization of quantum and classical correlations in open quantum FIG. 1. Open quantum Dicke model: superradiant phase transition and entanglement. (a) An ensemble of N spin-1/2 systems (with energy splitting ωz between upstate |↑ and down-state |↓ ) is coupled to a bosonic mode (the frequency ωm determines the energy-cost of creating one field excitation |n → |n + 1 ). The presence of an environment induces boson losses (at rate κ) and local spin-decay (at rate γ). (b) At a critical coupling strength λ = λc, the system undergoes a superradiant phase transition, characterized by a macroscopic occupation (∝ N ) of the bosonic mode, both in the presence and in the absence of local spin-decay. (c) The presence of local spin-decay leads to stronger quantum correlations and also "stabilizes" an entangled nonequilibrium superradiant phase.
Dicke models. We achieve this by developing an approach for treating quantum correlations, which is based on the theory of quantum fluctuation operators [23][24][25][26][27][28], that can be applied also in cases where the Holstein-Primakoff approximation cannot be exploited, e. g. in the presence of local dissipative terms which prevent a representation of the dynamics through collective spin operators. We fo-cus on various correlation measures, such as quantum discord and classical correlation, and show that they display non-analytic behavior at the critical coupling strength [see e. g. Fig. 1(c)]. Furthermore, we analyze bipartite entanglement between the spins and the bosonic mode. We find that the presence of local spin-decay -an unavoidable process in experiments which is usually considered detrimental for quantum effects -is unexpectedly beneficial for the build-up of quantum correlations. Our results indicate that this process leads to increased entanglement in the normal phase, and to the emergence of a nonequilibrium superradiant phase [cf. Fig. 1(c)] where entanglement is nonvanishing. Open quantum Dicke model.-The Dicke model consists of an ensemble of N spin-1/2 subsystems collectively interacting with a single bosonic mode, see Fig. 1(a). Spin operators for the kth particle are denoted as σ α k , with σ x = (|↑ ↓| + |↓ ↑|)/2, σ z = (|↑ ↑| − |↓ ↓|)/2, and σ y = −2iσ z σ x . Here, the states |↑ , |↓ are the single-particle spin states. The bosonic mode is described by creation and annihilation operators a † and a, respectively. For later convenience, we also introduce the spin operators σ ± = σ x ±iσ y and the bosonic quadrature operators q = i(a − a † )/ √ 2 and p = (a + a † )/ √ 2. The system is governed by a Markovian open quantum dynamics, under which the time-evolution of an operator O follows the Lindblad equationȮ(t) = L N [O(t)] [29][30][31] with generator The first term on the right-hand side of Eq. (1) gives the coherent contribution to the dynamics implemented by the Dicke Hamiltonian (setting = 1) Here, ω m > 0 is the bosonic mode frequency, ω z > 0 the energy splitting between spin states and λ > 0 the coupling parameter [cf. Fig. 1(a)]. The Dicke Hamiltonian (2) is written in terms of the collective operators where αβγ is the Levi-Civita symbol. The factor 1/ √ N , which rescales the collective spin-boson coupling in H D N , is necessary for a well-defined thermodynamic limit [8]. The last two terms in Eq. (1) account for irreversible dynamical effects. These are decay of bosonic excitations at rate κ as well as local (individual) spin-decay, |↑ → |↓ , at rate γ. As becomes clear from Eq. (1), the latter process is not described through collective, but rather local, spin (jump) operators σ − k .

Superradiant transition and mean-field results.-
The open quantum Dicke model undergoes a phase transition -as a function of the coupling strength λfrom a normal stationary phase, with subextensive (in N ) bosonic occupation, to a superradiant one where the bosonic mode becomes macroscopically occupied [8,18,19] [see sketch in Fig. 1(b)]. An order parameter for this transition is the stationary expectation of the renormalized number operator a † a/N in the thermodynamic limit of large number of spins, N → ∞.
The form of the order parameter suggests the definition of the so-called mean-field operators where the last term must be considered for α = x, y, z.
The first two operators are relevant as they provide the order parameter through (m q N ) 2 + (m p N ) 2 − 1/N = 2a † a/N . The mean-field operators of the spin ensemble [last terms in Eqs. (3)] are also important for studying the model. Indeed, by computing the action of the generator L N in Eq. (1) on the mean-field operators in Eqs. (3), one finds that these operators are all dynamically coupled [19]. In the thermodynamic limit, the time-evolved operators m α N (t) (with α = x, y, z, q, p) behave as multiples of the identity proportional to their expectation, [18,19]. Furthermore, they obey the differential equations (we drop the explicit time-dependence) These equations feature two different stationary regimes, separated by a critical value of the coupling strength For λ < λ c , there exists a unique stable stationary solution to the system in Eqs. (4), with the only nonzero value given by m z = −1/2. This is the normal phase. Quantum fluctuations.-The observables m α N provide, in the thermodynamic limit, N → ∞, a classical (mean-field) description of the Dicke model [19,28,41] which carries no information about correlations. In order to go beyond this, we introduce a new set of observables, so-called quantum fluctuation operators [23][24][25][26][27][28]. These will allow us to explore fluctuations and quantum correlations in the two stationary phases.
The fluctuation operators read For α = x, y, z, these are the usual spin fluctuation operators [27,28], and we have defined the bosonic ones (α = q, p) in full analogy. Roughly speaking, the operators in Eqs. (6) account for deviations of the operators m α N from the mean-field behavior. Remarkably, despite being collective, these retain a quantum character in the thermodynamic limit, in which the limiting operators F α = lim N →∞ F α N behave as bosons (for a rigorous discussion see e. g. Ref. [26]). This is straightforward to check for F q/p , since F q N = q − q and F p N = p − p . However, also collective spin fluctuations give rise to an emergent bosonic mode. This can be seen as follows. Looking at the commutator of fluctuation operators, one finds that [F α , F β ] = i γ αβγ m γ (α, β = x, y, z), which is a multiple of the identity. Now, we rotate the reference frame for the spin ensemble aligning the z-direction with the direction identified by mean-field variables, which iŝ n = m s /| m s | with m s = (m x , m y , m z ) T . In this rotated frame we havem x =m y = 0 andm z > 0, so that the only nonzero commutator is [F x ,F y ] = im z . A canonical bosonic mode is finally obtained by defining In what follows, we work with the set of fluctuations r = (Q, P,F z , F q , F p ) T . The first two elements represent an emergent bosonic mode describing collective properties of the spin ensemble; the last two are the fluctuations of the original bosonic mode, whileF z is a fluctuation operator which commutes with the others [26,28].
To analyze correlations in the Dicke model through fluctuation operators, we introduce the covariance ma-trixΣ αβ = {r α , r β } /2. For Gaussian states, this matrix contains the full information about fluctuations and can even be used to quantify collective correlations [42,43]. Before going to that, however, we briefly discuss the timeevolution ofΣ under the dynamics implemented by the generator in Eq. (1). For each parameter regime, we consider the dynamics of fluctuations emerging, in the thermodynamic limit, from an initial state which is stationary with respect to the mean-field observables and possesses Gaussian fluctuations. In this setting, the covariance matrix obeys the differential equation [28,44] where the matricesG andD, whose explicit form is given in [33], depend on the parameters of the model and on the stable stationary mean-field variables of Eqs. (4). The time-evolution in Eq. (7) has the structure of a bosonic Gaussian open quantum dynamics [45], suggesting that the Gaussianity of fluctuations is preserved at all times. Moreover, as long as λ = λ c , Eq. (7) has a unique stationary solutionΣ ∞ [33].
Since we are mainly interested in quantum correlations, we discard the information associated with the trivial fluctuationF z . This can be done by extracting from the stationary covariance matrixΣ ∞ the 4 × 4 minor obtained by neglecting its third row and its third column. The resulting covariance matrix contains the full information about the two bosonic modes Q, P and F q , F p . In particular, Γ s is the 2 × 2 matrix containing the second moments of the operators Q, P , Γ b contains those of F q , F p , and Γ c contains correlations between Q, P and F q , F p . Quantum and classical correlations.-In order to explore the correlation structure in the open quantum Dicke model, we focus on measures which can distinguish between correlations of different nature, e. g. quantum or classical, and that are fully determined by the covariance matrixΣ t−m ∞ [33]. Since the spin fluctuation operators involve all the spin degrees of freedom, the correlations that we discuss here are of collective type, i. e. reflecting how the spin ensemble as a whole is collectively correlated with the bosonic mode.
Firstly, we consider the classical correlation J [33, [46][47][48][49][50] between the spin ensemble and the bosonic mode. This quantity encodes the maximum information that can be extracted on one subsystem, by making generalized (Gaussian) measurements on the other one. In this sense, the classical correlation is "asymmetric" since it can be defined in two ways, i. e. either considering that measurements are performed on the spin ensemble or on the bosonic mode. Secondly, we study the so-called quantum discord D [33, [46][47][48][49][50], which is defined as the difference between the total correlation -quantified by the quantum mutual information -and the classical correlation J . This quantity measures the genuine quantum contribution to the total correlation between the two subsystems. According to its definition through the classical correlation, also the quantum discord is asymmetric under exchange of the role of the spin ensemble and of the bosonic mode. In the following, we consider both quantum discord and classical correlation assuming that the measurements are performed on the bosonic mode (see results in [33] for the other case).
In Fig. 2(a-b), we show the stationary behavior of classical correlation and quantum discord, as a function of the coupling strength λ and of the local spindecay rate γ. As shown, the classical correlation diverges at the nonequilibrium phase transition line, witnessing strong spin-boson correlations. Concerning the presence of quantum correlations, we observe that quantum discord D is different from zero almost everywhere in the phase diagram. It is maximal along the critical line, where it shows a non-analytic behavior even though it remains bounded.
We now consider the emergence of collective entanglement between the spins and the bosonic mode. This can be quantified from the covariance matrixΣ t−m ∞ , through the logarithmic negativity E N -a proper entanglement measure -defined as [33, [51][52][53] Here,ν − is the smallest (symplectic) eigenvalue [53] of the partially transposed covariance matrix obtained from 2Σ t−m ∞ by exchanging F p → −F p [39,40,54,55]. As we show in Fig. 2(c), the open quantum Dicke model displays collective spin-boson entanglement in a large parameter regime. We are particularly interested in understanding the impact of local spin-decay on entanglement. For small, yet nonvanishing, values of γ we identify a pronounced peak near the critical line λ c (γ). This suggests that a small rate of local spin-decay leads to larger entanglement. However, when γ vanishes, entanglement is dramatically reduced. This becomes evident when comparing the behavior of entanglement, as a function of λ, for γ = 0 and γ = 0. An example is shown in Fig. 2 In the absence of local spin-decay, entanglement vanishes at the critical point and is always zero in the superradiant phase [cf. Fig. 2(d)]. However, when local spin-decay is present, entanglement assumes larger values across the whole phase diagram and can also persist in the superradiant phase [cf. Fig. 2(e)]. Furthermore, for γ = 0, the logarithmic negativity shows a non-analytic behavior at the critical point and undergoes a "sudden death" well inside the superradiant phase, as shown in Fig. 2(e). These results show that local spin-decay has, rather surprisingly, an overall beneficial effect on quantum correlations, and on quantum entanglement in particular. Comparing Fig. 2(b) and Fig. 2(c), we also see that there exist parameter regions where the quantum discord assumes a finite value but the logarithmic negativity is zero. In this region, the quantum state of fluctuations is separable but nevertheless non-trivially quantum correlated.
Finally, we analyze quantum correlations within each subsystem separately. These are measured by the squeezing parameter ξ = 2 min(Θ 1 , Θ 2 ) [44,[56][57][58] where Θ 1 , Θ 2 denote the eigenvalues of Γ s for spin squeezing and of Γ b for boson squeezing. The parameter ξ quantifies the minimum variance among all possible quadrature operators. A state is called squeezed if ξ < 1, i. e. if the variance in one of the quadratures is smaller than the smallest possible simultaneous uncertainty of two canonically conjugated quadrature operators (also referred to as shot-noise limit [58]). Fig. 3(a) shows that there is no spin-squeezing in the stationary state of the model since ξ s is always larger than or equal to 1. In contrast, the system can feature squeezing in the bosonic mode below a threshold, i. e. γ 4. As shown in the inset of Fig. 3(b) the bosonic squeezing parameter ξ b takes its minimum values near λ = λ c (γ).
Discussion.-We explored the stationary structure of correlations in an open quantum Dicke model. We found that, in the absence of local spin-decay (γ = 0), the superradiant phase does not feature spin-boson entanglement. Even though this may appear somehow counterintuitive since superradiant phases arise in the strong coupling regime [59], disentangled nonequilibrium superradiant phases have also been observed in other settings [60]. However, as we have shown, the presence of local spin-decay (γ = 0) appears to be beneficial for the buildup of quantum correlations and can even be used to "stabilize" entanglement in superradiant stationary regimes. Furthermore, through other measures of correlations, we have shown that, even when there is no spin-boson entanglement, there are residual quantum correlations in the system which evidence non-classical properties across the whole phase diagram of the open quantum Dicke model. Also these correlations, and not only entanglement, could be exploited to achieve quantum-enhanced sensitivity in metrological applications [61].

MEAN-FIELD ANALYSIS
In this section, we show how the mean-field equations, Eqs. (4), can be derived from the Lindblad generator L N . We then obtain the stationary solutions for such equations and analyze their stability. Summation over repeated indices is implied here and in the following.

Derivation of the mean-field equations
Starting point for the derivation of Eqs. (4) is the superoperator in Eq. (1), which can be rewritten as , where we have used that for an arbitrary operator A We separate the problem of calculating the action of L N on the mean-field operators, by looking separately at the actions of the contributions i[H D N , X], D κ N [X] and D γ N [X]. For this, we note that From the latter we get for α = x, y, z, as well as For the dissipator D κ N , we immediately see that D κ N [m α N ] = 0, for α = x, y, z, and that For the dissipator D γ N , in contrast, we have that D γ N [m q N ] = 0 and D γ N [m p N ] = 0. Using the identities [σ α k , σ and σ ρ k σ ν k = δ ρν 1 k 4 + i ρνµ σ µ k 2 it also follows that , taking the expectation value and using that in the large N limit m α N → m α which are proportional to the identity, one obtains the mean-field equations reported in Eqs. (4).

Stationary state solution of the mean-field equations
Setting the time derivatives of the equations in (4) to zero, we get the system of equations which determines the stationary solution. Through several substitutions we can find Including also Eq. (S4), we further extract Thus with Eqs. (S2), (S4) and (S6) we see that for γ > 0 Considering the case γ = 0, the latter solution is still a valid solution whenever λ ≥ 0 (together with the solution having m z = 1/2) but the non-trivial solutions are now given by for λ ≥ λ c . Since the sign choice for m z is independent of the others, these are four solutions and the critical coupling λ c is the same as above, evaluated at γ = 0.

Stability analysis of the stationary mean-field solutions
For the solutions of the last section we perform a stability analysis using Lyapunov's indirect method [35]. For γ = 0 the constraint m z 2 = 1 4 − m y 2 − m x2 , occuring due to the conservation of S 2 , reduces Eqs. (4) to a system of four coupled first-order non-linear differential equationṡ which completely determines the dynamics. It can be written in the form˙ u = f ( u) with u = (m x , m y , m q , m p ) T and the Jacobian of f is The characteristic polynomial is calculated as Then, Hurwitz' theorem [34,36] states that all roots of this polynomial have negative real parts if and only if the inequalities hold. Employing this theorem we see that for our choice of parameters ω z = 4, ω m = 1, κ = 1 all roots of the characteristic polynomial have a negative real part if and only if m z = −1/2 and 0 < λ < λ c in the trivial stationary solution and m z = −λ 2 c /(2λ 2 ) and λ > λ c in the non-trivial solutions. Then the respective stationary solutions are asymptotically stable [35]. Moreover, numerical evidence shows that for λ = λ c small perturbations of the trivial solution (coinciding at this point with the non-trivial solutions) as initial conditions of the dynamics still drive the system to the trivial solution. We proceed with the case γ > 0. Here the Jacobian is All roots of this degree 5 polynomial have negative real parts (which is a necessary condition for stability of the solution) [34,36] if and only if the inequalities hold. For ω z = 4, ω m = 1, κ = 1 the trivial solution is asymptotically stable for 0 ≤ λ < λ c . On the other hand, for λ > λ c the non-trivial solutions are stable. Also in this case there is numerical evidence that for λ = λ c the stationary solution is approached eventually.

TIME-EVOLUTION OF THE COVARIANCE MATRIX
In this section, we give the derivation of the dynamics of the covariance matrix reported in the main text. We will first derive the dynamics of the covariance matrix in the original frame and then show how this is modified when considering the emergent normal mode for the collective spin fluctuations. Finally, we discuss how the asymptotic covariance matrix can be found from the dynamical equation.

Derivation of the dynamics for fluctuations
We start considering the original fluctuation operators that we collect in the following vector (S10) Defining the two-point functions C αβ N := F α N F β N , the entries of the covariance matrix in the original frame Σ are given by (S11) In general, the operators of the form F α N F β N possess an explicit time-dependence through expectation values over the state contained in the definition of fluctuations in Eq. (S10). Taking the total time-derivative of the two-point functions we thus get Here, we used that F α N = 0 by definition and that d/dtF α N is a scalar quantity. We thus have thaṫ Thus the problem is, like in the mean-field analysis, separated into three parts. First we want to evaluate the first two terms. To this end, we consider the commutator of the Dicke Hamiltonian and the fluctuation vector components and Eq. (S1) whence Exploiting the fact that F α N = 0, we can write In the above calculation we have used that m α N → m α , multiple of the identity, in the thermodynamic limit. Analogously, we can calculate For the remaining two parts ofĊ N we expand and Focusing on D κ N , we find for the last term on the right-hand side where the symplectic matrix s N is given by the commutation relations of the fluctuation operators s αβ The single-fluctuation κ-dissipator reads and therefore Collecting intermediately all the results concerning D κ N , we see Proceeding with D γ N , we have for the single-fluctuation dissipator and For the last term in Eq. (S13) we get and by means of σ ρ k σ ν k = δ ρν 1 k 4 + i ρνµ σ µ k 2 , Here . Now also collecting the results concerning D γ N gives and we conclude for Eq. (S12) in the thermodynamic limiṫ C = CA T + AC + CEs + sEC − sD s + CQ + QC + Z .
Considering then Eq. (S11), we finally get the differential equation for the covariance matriẋ We note that the differential equation for the covariance matrix involves, in general, time-dependent matrices G, s and Z. These may indeed be time-dependent through the time-dependence of the mean-field operators which appear in their matrix elements. However, for the purpose of this work, the matrices G, s and Z are time-independent since we investigate here the behavior of fluctuations when the state of the system is already stationary with respect to the mean-field observables.

Emergent normal mode
We transform the fluctuation vector as The J-matrix realizes a rescaling of the two remaining non-classical collective spin degrees of freedom, obtained after rotating, such that indeeds Similarly the covariance matrix transforms asΣ = JR θ,ϕ ΣR T θ,ϕ J and the time-evolution is now given by the differential equationΣ where for the sake of clarity we havẽ We stress here again that the considered initial state for the system is stationary with respect to the mean-field observables so that also R θ,ϕ and J are time-independent.

Stationary covariance matrix
In this section we show how the stationary covariance matrix can be obtained through a vectorization procedure. We will refer to an odd-dimensional square matrix as in "cross" form if its middle row and column consist only of zeros. An even-dimensional square matrix arising from an odd-dimensional one M by deleting the middle row and column is said to be in "reduced" form and we denote it by M red . We start with the differential equation for the covariance matrix from the last sectioṅ Σ(t) =Σ(t)G T +GΣ(t) +W (S14) where we definedW =Z −sDs. We focus on parameters chosen for Fig. 2 in the main text, i. e. ω z = 4, ω m = 1, κ = 1. Let I ⊆ R be an open interval. Here t ∈ I and t > t 0 ∈ I. At t 0 it is assumed that the quantum state is such that where m sub is the vector containing the stable solution of the mean-field equations in the normal phase, while m sup is the vector containing the solution of the mean-field equations in the superradiant phase. The cases γ = 0 and γ > 0 are treated separately and we first focus on γ > 0. The task is to find the stationary covariance matrixΣ ∞ which is such thatΣ ∞ = 0. The matrix equation to be solved, given bỹ is equivalent [37] to finding the 25 unknowns of the following linear system of 25 equations where ⊗ is the Kronecker product and the operation vec(·) arranges the entries of a matrix columnwise in a vector top down. Eq. (S15) has a unique solution if and only ifG ⊗ 1 5 + 1 5 ⊗G is invertible. The solution is, in vectorized form, given by vec(Σ ∞ ) = (G ⊗ 1 5 + 1 5 ⊗G) −1 vec(−W ).
Equivalently to invertibility, we want to prove that any eigenvalue ofG ⊗ 1 5 + 1 5 ⊗G is nonzero. If the spectrum ofG is σ(G) = {µ 1 , µ 2 , µ 3 , µ 4 , µ 5 } then the set of these eigenvalues is σ(G ⊗ 1 5 + 1 5 ⊗G) = {µ i + µ j |i = 1, ..., 5, j = 1, ..., 5}. Thus, any eigenvalue is nonzero if σ(G) ∩ σ(−G) = ∅, i.e. if no element of σ(−G) can be obtained by a point reflection of an element of σ(G) at the origin of the complex plane. Using again Hurwitz' theorem it can be proven that for λ = λ c all eigenvalues ofG lie in the open left half-plane. Consequently the matrixG In the γ = 0 case one cannot proceed the same way. In this setting, we focus on initial covariance matricesΣ(t 0 ) that are in cross form (see the definition at the beginning of this subsection). The matrixZ is the zero matrix and the differential equation (S14) reduces toΣ (t) = −sDs +Σ(t)G T +GΣ(t).
Known as the differential Sylvester equation [38], it has the unique solutioñ We note that This is in cross form sinceΣ(t 0 ) is in this form and thusG kΣ (t 0 )(G T ) l is in cross form, for all k, l ∈ N 0 . Similarly, sincesDs is in cross form, also eG (t−s)sDs eG T (t−s) is. It follows that the unique solution in Eq. (S16) has cross form for all t ∈ I with t > t 0 . Therefore it remains to solvė We want to find the stationary covariance matrixΣ t−m ∞ (t) = 0. Solving the matrix equatioñ is equivalent to finding the 16 unknowns of the linear system of 16 equations With the same steps as above, we establish the invertibility ofG red ⊗ 1 4 + 1 4 ⊗G red for λ ∈ {0, λ c } such that in this regime the unique stationary CM is given by vec(Σ t−m ∞ ) = (G red ⊗ 1 4 + 1 4 ⊗G red ) −1 vec(s redDredsred ).

QUANTUM AND CLASSICAL CORRELATIONS FOR TWO-MODE BOSONIC GAUSSIAN STATES
Here we discuss important tools of Gaussian quantum information theory and define the correlation measures exploited in the main text. In addition, we report supplementary results.

Measures of correlations for Gaussian bosonic systems
The total correlations contained in a quantum state ρ s,b (the notation reflects the bipartition into a spin subsystem and a boson subsystem) can be measured by the quantum mutual information , otherwise and f (x) = x + 1 2 log We used that, according to Williamson's theorem, ν + , ν − are the pairwise occuring symplectic eigenvalues of the covariance matrix, obtained as the diagonal elements of the symplectic diagonalized matrix 2Σ t−m ∞ . Based on Eq. (S17) one can also study entanglement in the system. In terms of the smallest symplectic eigenvalue ν − , a physically permissible (bonafide) covariance matrix has to fulfill ν − ≥ 1. The violation-degree of this condition under partial transposition of the underlying density matrix can be quantified by the logarithmic negativity E N = max(0, − log(ν − )).
Hereν − is the smallest symplectic eigenvalue of the partially transposed covariance matrix, obtained from 2Σ t−m ∞ by exchanging F p → −F p . Non-zero values of E N are necessary and sufficient for entanglement (non-separability) between the spin ensemble and the bosonic mode as a result of the Peres-Horodecki criterion.
Stationary quantum discord and classical correlation for measurements on the spin system As already mentioned in the main text, the quantum discord and the classical correlation allow for analog definitions with respect to measurements on the spins. Within the discussion of Gaussian states, in terms of covariance matrices, these definitions can be achieved by accordingly interchanging det(2Γ s ) ↔ det(2Γ b ) in the formulae given above.
The parameter dependence of J and D in this case is illustrated in Fig. S1. We can see from this figure that the discord has now maxima of approximately the same height, still distributed along the critical line λ c (γ). In contrast, the classical correlation shows essentially the same behavior as for the case where the measurements were performed on the boson system. It still diverges at the critical line (in the plots the asymptotic value is bounded by the chosen parameter resolution).
FIG. S1. Quantum discord and classical correlation for spin-measurements. (a) Classical correlation J and (b) quantum discord D as functions of γ and λ. Here J was maximized over all POVMs on the spin system. The insets visualize the λ-dependence of the respective quantities for γ = 2. The plots were produced assuming that ωm = 1 and ωz = 4. All parameters are in units of κ.