Critical phenomena and nonlinear dynamics in a spin ensemble strongly coupled to a cavity. II. Semiclassical-to-quantum boundary

TL;DR

This paper investigates the quantum-to-classical transition in a driven-dissipative spin ensemble coupled to a cavity, revealing significant quantum effects near critical points even in large ensembles of up to 10,000 spins.

Contribution

It introduces a generalized cumulant expansion method to analyze finite-size quantum corrections and the validity of semiclassical equations in a driven spin-cavity system.

Findings

Quantum corrections remain significant near critical points.

Semiclassical equations are invalid close to phase transition thresholds.

Finite-size effects are crucial for understanding the quantum-classical boundary.

Abstract

We numerically study the dynamics and stationary states of a spin ensemble strongly coupled to a single-mode resonator subjected to loss and external driving. Employing a generalized cumulant expansion approach we analyze finite-size corrections to a semiclassical description of amplitude bistability, which is a paradigm example of a driven-dissipative phase transition. Our theoretical model allows us to include inhomogeneous broadening of the spin ensemble and to capture in which way the quantum corrections approach the semiclassical limit for increasing ensemble size $N$. We set up a criterion for the validity of the Maxwell-Bloch equations and show that close to the critical point of amplitude bistability even very large spin ensembles consisting of up to $10^4$ spins feature significant deviations from the semiclassical theory.

Figures

Equations

$$\begin{gathered} \mathcal{H}={\Delta }_c{a}^{†}a+\frac{1}{2}\sum _{j=1}^N{\Delta }_j{\sigma }_j^z+\sum _{j=1}^N[g_j{\sigma }_j^{-}{a}^{†}+g_j^{\ast }{\sigma }_j^{+}a] \\ +i[\eta (t)a^{†}-{\eta }^{\ast }(t)a], \end{gathered}$$

$${\mathcal{L}}_D(\rho )=$$

$$+{\gamma }_h\sum _{j=1}^N(2{\sigma }_j^{-}\rho {\sigma }_j^{+}-{\sigma }_j^{+}{\sigma }_j^{-}\rho -\rho {\sigma }_j^{+}{\sigma }_j^{-}),$$

$$\frac{d}{dt}\rho =\frac{1}{i}[\mathcal{H},\rho ]+{\mathcal{L}}_D(\rho ),$$

$$\frac{d}{dt}⟨a⟩=-(\kappa +i{\Delta }_c)⟨a⟩-i\sum _{j=1}^N{g}_j⟨{\sigma }_j^{-}⟩+\eta ,$$

$$\frac{d}{dt}⟨{\sigma }_j^{-}⟩=-({\gamma }_h+2{\gamma }_p+i{\Delta }_j)⟨{\sigma }_j^{-}⟩+i{g}_j⟨{\sigma }_j^{z}a⟩,$$

$$\frac{d}{dt}⟨{\sigma }_j^z⟩=-2{\gamma }_h(⟨{\sigma }_j^z⟩+1)+2i{g}_j(⟨{\sigma }_j^{-}{a}^{†}⟩-⟨{\sigma }_j^{+}a⟩).$$

$$⟨{\sigma }_j^{z}a⟩$$

$$⟨{\sigma }_j^{-}{a}^{†}⟩$$

$$\frac{d}{dt}⟨a⟩=-(\kappa +i{\Delta }_c)⟨a⟩-i\sum _{k=1}^N{g}_j⟨{\sigma }_j^{-}⟩+\eta ,$$

$$\frac{d}{dt}⟨{\sigma }_j^{-}⟩=-({\gamma }_h+2{\gamma }_p+i{\Delta }_j)⟨{\sigma }_j^{-}⟩+i{g}_j⟨{\sigma }_j^z⟩⟨a⟩,$$

$$\frac{d}{dt}⟨{\sigma }_j^z⟩=-2{\gamma }_h(⟨{\sigma }_j^z⟩+1)+2i{g}_j(⟨{\sigma }_j^{-}⟩⟨a^{†}⟩-\text{c.c.}),$$

$$⟨a⟩\to ⟨a⟩/\sqrt{N},g_j\to g_j/\sqrt{N},\eta \to \eta \sqrt{N},$$

$${⟨\hat{A}\hat{B}⟩}_c=⟨\hat{A}\hat{B}⟩-⟨\hat{A}⟩⟨\hat{B}⟩.$$

$${⟨\hat{A}\hat{B}\hat{C}⟩}_c=$$

$$+2⟨\hat{A}⟩⟨\hat{B}⟩⟨\hat{C}⟩,$$

$${⟨\hat{A}\hat{B}\hat{C}\hat{D}⟩}_c$$

$$+⟨\hat{C}⟩⟨\hat{A}\hat{B}\hat{D}⟩+⟨\hat{D}⟩⟨\hat{A}\hat{B}\hat{C}⟩+⟨\hat{A}\hat{B}⟩⟨\hat{C}\hat{D}⟩$$

$$\displaystyle+\braket{\hat{A}\hat{C}}\braket{\hat{B}\hat{D}}+\braket{\hat{A}\hat{D}}\braket{\hat{B}\hat{C}}\Bigr{)}$$

$$\displaystyle+2\Bigl{(}\braket{\hat{A}\hat{B}}\braket{\hat{C}}\braket{\hat{D}}+\braket{\hat{A}\hat{C}}\braket{\hat{B}}\braket{\hat{D}}$$

$$+⟨\hat{A}\hat{D}⟩⟨\hat{B}⟩⟨\hat{C}⟩+⟨\hat{B}\hat{C}⟩⟨\hat{A}⟩⟨\hat{D}⟩$$

$$\displaystyle+\braket{\hat{B}\hat{D}}\braket{\hat{A}}\braket{\hat{C}}+\braket{\hat{C}\hat{D}}\braket{\hat{A}}\braket{\hat{B}}\Bigr{)}$$

$$-6⟨\hat{A}⟩⟨\hat{B}⟩⟨\hat{C}⟩⟨\hat{D}⟩.$$

$$⟨\hat{A}\hat{B}\hat{C}⟩\approx$$

$$-2⟨\hat{A}⟩⟨\hat{B}⟩⟨\hat{C}⟩,$$

$$⟨\hat{A}\hat{B}\hat{C}\hat{D}⟩$$

$$+⟨\hat{D}⟩⟨\hat{A}\hat{B}\hat{C}⟩+⟨\hat{A}\hat{B}⟩⟨\hat{C}\hat{D}⟩+⟨\hat{A}\hat{C}⟩⟨\hat{B}\hat{D}⟩$$

$$\displaystyle+\braket{\hat{A}\hat{D}}\braket{\hat{B}\hat{C}}-2\Bigl{(}\braket{\hat{A}\hat{B}}\braket{\hat{C}}\braket{\hat{D}}$$

$$+⟨\hat{A}\hat{C}⟩⟨\hat{B}⟩⟨\hat{D}⟩+⟨\hat{A}\hat{D}⟩⟨\hat{B}⟩⟨\hat{C}⟩$$

$$+⟨\hat{B}\hat{C}⟩⟨\hat{A}⟩⟨\hat{D}⟩+⟨\hat{B}\hat{D}⟩⟨\hat{A}⟩⟨\hat{C}⟩$$

$$\displaystyle+\braket{\hat{C}\hat{D}}\braket{\hat{A}}\braket{\hat{B}}\Bigr{)}+6\braket{\hat{A}}\braket{\hat{B}}\braket{\hat{C}}\braket{\hat{D}}\,,$$

$$|⟨a_{st}⟩{|}^2{\left(1+C\frac{1}{1+|⟨a_{st}⟩{|}^2/n_0}\right)}^2={\eta }^2/{\kappa }^2$$

$$|⟨\widetilde{a_{st}}⟩{|}_{{}_{\text{(CE2,3)}}}^2\equiv \frac{|⟨a_{st}⟩{|}_{{}_{\text{(CE2,3)}}}^2}{|⟨a_{st}⟩{|}_{{}_{\text{(CE1)}}}^2},$$

$$\Delta^{\text{CE}}_{n-m}\equiv\frac{\Big{|}|\langle a_{st}\rangle|^{2}_{{}_{(n)}}\!-|\langle a_{st}\rangle|^{2}_{{}_{(m)}}\Big{|}}{|\langle a_{st}\rangle|^{2}_{{}_{(m)}}}\;,$$

$${\Delta }_{1-2}^{\text{CE}},{\Delta }_{2-3}^{\text{CE}},{\Delta }_{1-3}^{\text{CE}}<\delta ϵ,$$

$$M_{\mu }=\frac{N}{K}{e}^{-4\ln (2)\frac{{\Delta }_{\mu }^2}{{\Gamma }^2}},$$

$$C=\frac{g^2}{\kappa {\gamma }_h}\sum _{\mu =1}^L{M}_{\mu }\frac{1}{1+{\Delta }_{\mu }^2/{\gamma }_h^2},$$

$$\frac{d}{dt}⟨a⟩$$

$$\frac{d}{dt}⟨{\sigma }_k^{-}⟩$$

$$\frac{d}{dt}⟨{\sigma }_k^z⟩$$

Full text

I Introduction

The dynamics of open many-body quantum systems is of fundamental importance for various branches of physics Daley (2014). In particular, so-called hybrid quantum systems, whose technological relevance requires them to be open, have shifted to the center of attention over the last decade Xiang et al. (2012); Kurizki et al. (2015). As a prominent example, spin ensembles coupled to a cavity mode emerged as a powerful platform for quantum computation Tordrup et al. (2008); Wesenberg et al. (2009); Ping et al. (2012), quantum memories Simon et al. (2010); Kubo et al. (2012); Grezes et al. (2016), and quantum communication Duan et al. (2001); Kimble (2008). Beside their technological significance, the driven and dissipative character of spin-cavity systems makes them also well-suited to study fundamental aspects of non-equilibrium many-body physics Houck et al. (2012); Tomadin and Fazio (2010); Buchhold et al. (2013). Corresponding theoretical descriptions clearly profit from the fact that the interactions among the individual spins or atoms are mediated via few common cavity modes only. This leads to extremely long-ranged interactions and suppressed fluctuations Buchhold et al. (2013); Foss-Feig et al. (2017), often enabling accurate semiclassical descriptions of these systems.

Among semiclassical approaches the seminal Maxwell-Bloch equations play a distinguished role in quantum optics. Being based on neglecting the correlations between spins and the electromagnetic field they successfully describe many effects of lasers Haken (1985), superradiance Dicke (1954); Svidzinsky et al. (2015); Rose et al. (2017), critical slowing down Bonifacio and Meystre (1979); Angerer et al. (2017); Krimer (2019), and amplitude bistability Gibbs et al. (1976); Bonifacio and Lugiato (1976, 1978); Lugiato (1984); Brennecke et al. (2008); Angerer et al. (2017); Krimer (2019). In a separate companion paper we use the Maxwell-Bloch equations to study the dynamics of macroscopic spin ensembles near the critical point of amplitude bistability and analyze the effect of critical slowing down in the presence of inhomogeneous broadening Krimer (2019). Here, we will investigate to what extent the Maxwell-Bloch equations themselves are well-justified to describe the actual quantum dynamics of such systems on experimentally relevant time scales. This question is of particular interest near the critical point of a driven-dissipative phase transition as in the case of optical bistability. Here the semiclassical Maxwell-Bloch equations give rise to two stable steady state solutions combined with a hysteresis effect, while a full quantum treatment predicts a unique steady state that can deviate dramatically from semiclassical results Lugiato (1984). One major source of this deviation is the bimodal character of the quantum-mechanical probability distribution, which leads to a switching of the system between the two semiclassically stable branches Dombi et al. (2013); Fink et al. (2017); Casteels et al. (2017). Since, however, the switching time between the stable solutions diverges with the system size Fink et al. (2017); Casteels et al. (2017), the Maxwell-Bloch equations describe well how large systems evolve dynamically for experimentally observable time scales. The natural question to ask in this specific context is thus how many spins are required to constitute an ensemble that is sufficiently ”large” to enter the thermodynamic limit where the semiclassical solutions are sufficiently accurate to describe the dynamics of this quantum system close to the bistability region Rempe et al. (1991); Dombi et al. (2013).

To carry out such a comparison between the semi-classical and the quantum dynamics, the time evolution of a finite number of spins inside a cavity needs to be evaluated Sarkar and Satchell (1987); Rempe et al. (1991); Dombi et al. (2013); Fink et al. (2017); Shirai et al. (2018); Carmichael (2008); Dhar et al. (2018); Kirton and Keeling (2017); Mori (2013). Due to computational constraints, a full quantum mechanical treatment of cavity-spin systems is, however, limited to small spin ensembles or to few excitations in the system Carmichael (2008). Results for spin ensembles of up to $8$ spins Dombi et al. (2013) showed a rapid convergence of the quantum case to the semiclassical limit, but failed to provide a quantitative analysis of the corresponding boundary. Also further extensive theoretical efforts to describe spin-cavity systems of increasing ensemble size on a full quantum level Kirton and Keeling (2017); Shirai et al. (2018); Dhar et al. (2018) left the transition between the semiclassical and the full quantum case mostly uncharted. Since these two realms occupy opposing limits with respect to the system size, special techniques are required to bridge the gap between the microscopic and the macroscopic domain.

For a closed system of completely symmetric spin ensembles, e.g., the scaling of quantum corrections as a function of the number of two-level systems was studied using a WKB approach analyzing the eigenvalue spectrum Keeling (2009). In the present paper we study the onset of quantum corrections in the complementary setting of an open system with external driving and dissipation. Using a generalized cumulant expansion approach Kubo (1962); Leymann et al. (2013); Henschel et al. (2010); Leymann et al. (2014); Krämer and Ritsch (2015); Vardi et al. (2001); Vardi and Anglin (2001); Casteels et al. (2016), we establish a criterion for the validity of the semiclassical Maxwell-Bloch equations for a wide range of parameters taking the effect of inhomogeneous broadening explicitly into account.

Our paper is organized as follows: In Sec. II we introduce the model and derive the hierarchic set of equations of motion for expectation values (Sec. II.1). We review the semiclassical approximation (Sec. II.2) and present a generalized truncation scheme based on higher orders of cumulants (Sec. II.3). The effect of amplitude bistability is investigated in Sec. III starting with an homogeneous spin ensemble (Sec. III.1) and examining the semiclassical-to-quantum boundary (Sec. III.2). The presented analysis is then extended to the case of an inhomogeneously broadened spin ensemble (Sec. III.3). In Sec. IV we draw our conclusions.

II Model

The system we consider consists of $N$ two-level atoms or spins with transition frequencies $\omega_{j}$ coupled to a single-mode cavity with coupling strength $g_{j}$. The two-level emitters may or may not exhibit inhomogeneous broadening or coupling, both of which can be treated with the model at hand. The cavity is coherently driven by an external field of strength $\eta$ and frequency $\omega_{p}$. The starting point for the theoretical model is the Tavis-Cummings Hamiltonian (Tavis and Cummings, 1968), which, in a rotating frame with driving frequency $\omega_{p}$, reads ($\hbar=1$)

[TABLE]

where $\Delta_{c}\equiv\omega_{c}-\omega_{p}$ and $\Delta_{j}\equiv\omega_{j}-\omega_{p}$ are the detunings of the cavity frequency $\omega_{c}$ and of the individual spin frequencies $\omega_{j}$ with respect to the external driving field of frequency $\omega_{p}$. Here $a^{\dagger}$ and $a$ are the creation and annihilation operators of the single cavity mode and $\sigma_{j}^{z}$, $\sigma_{j}^{+}$, and $\sigma_{j}^{-}$ are the Pauli operators corresponding to the individual spins. In the following we consider the driving amplitude to be constant ($\eta(t)=\eta$) and without loss of generality assume $\eta^{*}=\eta$ as well as $g^{*}_{j}=g_{j}$.

The dissipation of the system is described by the Lindblad superoperator

[TABLE]

where the first term corresponds to the cavity loss at rate $\kappa$; the second term gives non-radiative dephasing of the individual spins with rate $\gamma_{p}/2$, and the last term describes their radiative decay with rate $\gamma_{h}$.

The total dynamics of the driven dissipative system is then given by the master equation

[TABLE]

with $\mathcal{H}$ and $\mathcal{L}_{D}$ given by Eq. (1) and (2), respectively, and $\rho$ being the density operator of the total cavity-spin system.

II.1 Equations of motion for expectation values

Full quantum solutions for the density operator of the system or of chosen subsystems can be obtained via direct integration of the master equation Carmichael et al. (1991); Auffèves et al. (2011); Sáez-Blázquez et al. (2017), quantum trajectory methods Daley (2014); Gardiner et al. (1992); Vukics and Ritsch (2007); Dombi et al. (2013); Fink et al. (2017) or, as recently shown in (Dhar et al., 2018) also by variational renormalization group methods. All of these approaches are limited, however, in the number of spins or excitations in the system. Since we are interested in the relation between quantum mechanical and semiclassical solutions over a wide range of parameters, we take a different approach here and directly solve for the expectation values of the operators of interest Leymann et al. (2014). Multiplying the master equation (3) with the given operator, taking the trace operation and using the cyclic permutation of the trace, it is straight-forward to obtain the following equations of motion (EoM) for the expectation values $\langle a\rangle$, $\langle\sigma_{j}^{-}\rangle$, and $\langle\sigma_{j}^{z}\rangle$:

[TABLE]

In the following we denote expectation values involving a product of n operators as nth order expectation values. The interaction part of the Hamiltonian couples EoM for nth-order expectation values to EoM for (n$+1$)th order expectation values and thereby creates an infinite hierarchy of coupled equations. (The EoM up to third-order expectation values are shown in the Appendix.) To solve the dynamics of the system, the hierarchy of equations thus needs to be truncated at some level. In the following we use a truncation procedure based on a cumulant expansion Kubo (1962); Henschel et al. (2010); Leymann et al. (2013, 2014); Krämer and Ritsch (2015) to obtain a closed set of equations that can be solved numerically.

II.2 Semiclassical approximation

The most prominent approach for driven dissipative spin-cavity systems of the type described in the previous section is to solve Eqs. (4)-(6) in the semiclassical limit by applying a self-consistent field approximation (Bonifacio and Lugiato, 1978),

[TABLE]

Using this full factorization, the hierarchy of equations (4)-(6) truncates at the 1st order according to our previous notation and one derives the seminal Maxwell-Bloch equations:

[TABLE]

where c.c. stands for the complex conjugate of the previous term. The full factorization of expectation values employed above is well justified for an infinite number of spins, since fluctuations in the large $N$ limit decrease with $1/N$ Mori (2013); Carmichael (2008). An important property of the semiclassical limit is that with the transformations

[TABLE]

and fixed spin distribution $\rho(\omega_{j})$, the governing equations are invariant under variation of the number of spins $N$ inside the ensemble. Any deviations from the scaling given by Eq. (12) can therefore be attributed to quantum corrections to the semiclassical equations (Dombi et al., 2013).

II.3 Cumulant expansion approach

The basis of the semiclassical limit is the full factorization given by Eqs. (7)-(8), neglecting the second order cumulant, which for two arbitrary operators $\hat{A}$ and $\hat{B}$ is defined as

[TABLE]

To include correlations between operators in the model it is necessary to keep the EoM for expectation values of higher order explicitly and truncate the hierarchy of equations at some higher level. The explicit formulas for third- and fourth-order cumulants read Kubo (1962)

[TABLE]

[TABLE]

The expressions for cumulants of high order are quite cumbersome but they allow us to represent an expectation value of a given order by a cumulant of the same order and by expectation values of lower order only. Therefore a closed system of equations can be obtained by neglecting higher orders of cumulants. In the following we call a truncation scheme that keeps EoM of expectation values up to nth-order but neglects cumulants of the (n$+1$)th-order as the nth-order cumulant expansion (CEn). The semiclassical Maxwell-Bloch equations, hence, can be considered as a first-order cumulant expansion (CE1).

For the present paper we apply a cumulant expansion of second (CE2) and third order (CE3) to the coupled EoM of the driven dissipative spin system (Appendix). The CE2 then consists of 12 coupled EoM for the expectation values: $\braket{a}$, $\braket{\sigma_{k}^{-}}$, $\braket{\sigma_{k}^{z}}$, $\braket{\sigma_{k}^{z}a}$, $\braket{\sigma_{k}^{z}\sigma_{j}^{-}}$, $\braket{\sigma_{k}^{-}a^{\dagger}}$, $\braket{\sigma_{k}^{+}\sigma_{j}^{-}}$, $\braket{\sigma_{k}^{-}a}$, $\braket{a^{\dagger}a^{\dagger}}$, $\braket{a^{\dagger}a}$, $\braket{\sigma_{k}^{z}\sigma_{j}^{z}}$, and $\braket{\sigma_{k}^{-}\sigma_{j}^{-}}$. All third-order expectation values are expanded according to

[TABLE]

where third-order cumulants $\braket{\hat{A}\hat{B}\hat{C}}_{c}$ are neglected. For the CE3 we eliminate this approximation and extend the EoM by $13$ additional equations for the expectation values $\braket{\sigma_{k}^{z}a^{\dagger}a}$, $\braket{\sigma_{k}^{-}a^{\dagger}a}$, $\braket{\sigma_{k}^{-}a^{\dagger}a^{\dagger}}$, $\braket{\sigma_{k}^{z}aa}$, $\braket{\sigma_{k}^{-}aa}$, $\braket{a^{\dagger}aa}$, $\braket{aaa}$, $\braket{\sigma_{k}^{z}\sigma_{j}^{z}a}$, $\braket{\sigma_{k}^{-}\sigma_{j}^{-}a^{\dagger}}$, $\braket{\sigma_{k}^{+}\sigma_{j}^{-}a}$, $\braket{\sigma_{k}^{z}\sigma_{j}^{-}a^{\dagger}}$, $\braket{\sigma_{k}^{z}\sigma_{j}^{-}a}$, and $\braket{\sigma_{k}^{-}\sigma_{j}^{-}a}$. Note that third-order expectation values containing only spin operators are not included and are truncated on the level of Eq. (16), which is justified since correlations among three spins play only a minor role as compared to correlations between spins and the collective cavity mode. All fourth-order expectation values that show up in the EoM of the CE3 are expanded as

[TABLE]

with fourth-order cumulants $\braket{\hat{A}\hat{B}\hat{C}\hat{D}}_{c}$ being neglected. The closed set of equations resulting from the CE2 and CE3 is solved numerically to obtain the dynamics and stationary states of the driven-dissipative spin system.

III amplitude bistability

We now analyze the validity of the semiclassical Maxwell-Bloch equations, using the CE2 and CE3 introduced in the previous section. In particular we will focus on amplitude bistability as a paradigm effect for cooperative phenomena in an open system far from equilibrium Lugiato (1984). Over the last decades amplitude bistability served as a role model for a non-equilibrium phase transition with experimental realizations in various systems Gibbs et al. (1976); Rempe et al. (1991); Brennecke et al. (2008); Angerer et al. (2017). To observe amplitude bistability in our model, we study the stationary transmission through a cavity-spin system under constant driving $\eta$. The transmission is proportional to the cavity probability amplitude $|\langle a\rangle|^{2}$, whose stationary value $|\langle a_{st}\rangle|^{2}$ can be obtained either directly by setting all time derivatives in Eqs. (9)-(11) to zero or by a temporal evolution of the system for a sufficiently long time.

III.1 Homogenous broadening

For simplicity, we start with the case of homogeneous coupling ($g_{j}=g$) and radiative decay only ($\gamma_{p}=0$). Moreover, we assume that the driving field is on resonance with the cavity and all individual spins ($\Delta_{c}=0$, $\Delta_{j}=0$). In the semiclassical limit, the stationary cavity probability amplitude $|\langle a_{st}\rangle|^{2}$ is then given by

[TABLE]

and is therefore completely characterized by the collective cooperativity parameter $C\equiv Ng^{2}/\kappa\gamma_{h}$, the photon saturation number $n_{0}\equiv\gamma_{h}^{2}/2g^{2}$, and the scaled driving amplitude $\eta/\kappa$. The stationary state equation above gives the typical, nonlinear input-output relations presented in Fig. 1 for collective cooperativity parameters ranging from $C=4$ to $14$. Note that we increase $C$ here by increasing the individual coupling $g$ while keeping all other parameters constant. At low driving, the enhanced cooperative emission of the ensemble leads to a stationary state of low transmission, sometimes called the lower or cooperative branch. Increasing cooperativities lead to an increasing suppression of the transmission (Bonifacio and Lugiato, 1978). For strong driving the spins start to saturate ($\langle\sigma^{-}_{j}\rangle_{st}\to 0$) and thereby decouple from the cavity (see Eq. (9)) leading to a stationary state of high transmission, which is independent of the cooperativity parameter $C$. This upper branch is therefore also called the independent-atom branch Carmichael (2008). For $C>8$ the stationary state equation (18) exhibits regions where three solutions coexist (two stable and one unstable). This bistable region is bounded by two critical points $\eta^{-}_{\text{crit}}$ and $\eta^{+}_{\text{crit}}$, which are characterized by an infinite slope $d|\langle a_{st}\rangle|^{2}/d\eta$ in the transmission curve. For $C<8$ no bistability occurs and only one point of maximal but finite slope $d|\langle a_{st}\rangle|^{2}/d\eta$ exists. In this case there is a continuous transition from the lower to the upper branch for increasing driving. In the bistable case, by contrast, the system changes discontinuously from a state of low transmission to a state of high transmission in a first-order phase transition at the critical point $\eta^{+}_{\text{crit}}$.

In the following we are interested in the validity of the semiclassical solution in the vicinity of the critical point $\eta^{+}_{\text{crit}}$ for a finite number of spins. Note that, under the trivial scaling given by Eq. (12), the semiclassical stationary state equation (18) is independent of the number of spins, $N$. The correlations $\langle\sigma_{j}^{z}a\rangle_{c}$ and $\langle\sigma_{j}^{-}a^{\dagger}\rangle_{c}$, however, can lead to significant deviations from this trivial scaling Dombi et al. (2013), as displayed in Fig. 2. Here we present the impact of the quantum corrections on the stationary cavity probability amplitude $|\langle a_{st}\rangle|^{2}$ for a collective cooperativity parameter of $C=14$ (with the tendency described in the following being similar for all $C$). Typical numerical results using the CE2 and CE3 are demonstrated on the example of ten different driving strengths, which are chosen such that we probe the stationary states on the lower transmission branch of the bistable regime as well as on the upper transmission branch above the critical point $\eta^{+}_{\text{crit}}$.

Our results in Fig. 2 show that, for small ensembles, the stationary transmission calculated by means of the CE2 and CE3 deviates significantly from the semiclassical solution even outside the bistable region. As expected, increasing the number of spins restores the semiclassical results, since the quantum fluctuations decrease as $1/N$ Dombi et al. (2013). However, the actual number of spins needed for the CE to agree well with the semiclassical solution substantially increases for driving strengths close to the critical point $\eta^{+}_{\text{crit}}$. Whereas the results obtained from the CE3 for spin ensembles of moderate size ($N=250$) agree reasonably well with the semiclassical solution for most driving strengths, this is not the case for the driving strengths close to the critical point of the low transmission branch $0.97\,\eta^{+}_{\text{crit}}$ and $0.99\,\eta^{+}_{\text{crit}}$, respectively. Here much larger numbers of spins ($N=2500$) are needed for the CE3 to reasonably approach the semiclassical lower transmission branch. Note that the stationary states shown in Fig. 2 are extracted from a temporal evolution of the system for sufficiently long time starting initially from an unexcited spin ensemble ($\langle\sigma_{j}^{z}\rangle=-1$, $\langle\sigma_{j}^{-}\rangle=0$) and an empty cavity ($\langle a\rangle=0$), subjected to constant driving. For this initial state the semiclassical dynamics converges towards the lower transmission branch of the bistable region and reaches the upper transmission branch only for $\eta>\eta^{+}_{\text{crit}}$ as indicated in Fig. 1.

Figure 2(c) presents the transient dynamics of the cavity probability amplitude $|\langle a(t)\rangle|^{2}_{{}_{\text{\!(CE3)}}}$ calculated using the CE3 for $C=14$ and driving strength $\eta=1.05\,\eta^{+}_{\text{crit}}$. As depicted already in Fig. 2(b) for $N=250$ the CE3 agrees well with the semiclassical solution, whereas for small spin ensembles ($N=10$ or $50$) the CE3 tends towards a stationary state of much lower transmission than that predicted by the semiclassical equations. Interestingly, our calculations for $N=86$ and $N=87$ spins indicate, that there is an abrupt crossover from spin ensembles with large deviations from the semiclassical limit towards spin ensembles where such deviations are small. Whereas for $N\leq 86$ the CE3 tends towards a stationary state of relatively low transmission, ensembles of $N\geq 87$ spins start to approach a state of high transmission, in accordance with the semiclassical solution.

To measure the validity of the semiclassical approximation for different driving strengths and cooperativity parameters, we normalize the stationary state obtained in the framework of the CE2 and CE3 by the corresponding semiclassical result (CE1),

[TABLE]

such that a value close to unity corresponds to the semiclassical regime. Figure 3 shows the normalized cavity probability amplitude $|\langle\widetilde{a_{st}}\rangle|^{2}_{{}_{\text{\!(CE2,3)}}}$ for the driving strength $1.05\,\eta^{+}_{\text{crit}}$ and cooperativity parameters ranging from $C=2$ to $C=20$. Focusing at first on the CE2 solutions, we can see that the normalized cavity probability amplitude approaches the semiclassical result $|\langle\widetilde{a_{st}}\rangle|^{2}_{{}_{\text{\!(CE2)}}}=1$ for increasing numbers of spins as expected from a linearized theory of fluctuations Carmichael (2008). As indicated already above, our results show that for increasing cooperativity parameters the transition towards the semiclassical solution becomes more abrupt and for $C>8$ resembles a first-order phase transition at $N\approx 45$.

Coming now to the CE3 solutions, a similar tendency is observed with the only difference that the transition is shifted to larger values of $N$ for increasing cooperativities $C$. Although these results suggest that for amplitude bistability the crossover from systems of large fluctuations towards systems of small fluctuations has a discontinuous nature, care must be taken, since in the crossover region the cumulant expansion has not yet converged, i.e. the CE2 and CE3 give quite different values for $|\langle a_{st}\rangle|^{2}$. Including higher orders of cumulants or a full quantum mechanical treatment of the problem is therefore required to ensure the convergence to a true quantum solution in this parameter regime.

It turns out that at values of $\eta$ slightly smaller than $\eta_{\text{crit}}^{+}$ (at which the first-order transition obtained in the framework of the semiclassical approach occurs), other time-dependent solutions can simultaneously exist for certain numbers of spins, $N$ - a common scenario in the driven dissipative dynamics described by a set of nonlinear differential equations. Specifically, starting from the simple initial conditions mentioned above (empty cavity with unexcited spin ensemble), we end up with periodic solutions after some transient time or sometimes the overall approach even becomes numerically unstable giving rise to unphysical solutions (indicated, e.g., by unphysical values of $|\langle\sigma_{k}^{z}\rangle|>1$). To overcome this problem we vary the initial conditions for $\langle\sigma_{k}^{z}\rangle$ between $-1$ and $-0.5$ to finally find those which lie in the so-called basin of attraction for the stationary state. Originating from such initial conditions, the system eventually settles to the stationary state under study as time increases.

III.2 Semiclassical-to-quantum boundary

To avoid the difficulties in the crossover region of Fig. 3, we focus in the following on the results of the CE close to the semiclassical stationary states and define a criterion for the validity of the semiclassical Maxwell-Bloch equations based on the convergence of the cumulant expansion. For this purpose we define the relative deviations

[TABLE]

where $n$ and $m$ stand for the different orders of the cumulant expansion, CE1, CE2, and CE3. Figure 4 presents the relative deviations between the first three orders of the cumulant expansion, $\Delta^{\text{CE}}_{1-2}$, $\Delta^{\text{CE}}_{2-3}$, and $\Delta^{\text{CE}}_{1-3}$, for $C=14$ and driving strengths $\eta=1.05\,\eta^{+}_{\text{crit}}$ and $0.95\,\eta^{+}_{\text{crit}}$, respectively. It turns out that the discontinuous nature of the crossover region is reflected also in the relative deviations, resulting in a rather complicated dependence on the number of spins $N$. Above the crossover regions the relative deviations decrease linearly with 1\big{/}N as expected from a linearized theory of quantum fluctuations in the small noise limit (Carmichael, 2008). The size of the relative deviations, however, does not only strongly depend on the cooperativity $C$ and the number of spins $N$ but also on the driving strength $\eta$. Our results show that for the same $C$ and the same $N$ the relative deviations for $\eta=0.95\,\eta^{+}_{\text{crit}}$ are significantly larger than for $\eta=1.05\,\eta^{+}_{\text{crit}}$. This asymmetry will be explored in more detail below.

In the following we define a small threshold value $\delta\epsilon=10^{-2}$, which serves as criterion for the convergence of the cumulant expansion as well as for the validity of the semiclassical solutions. We estimate the minimal number of spins, $N_{\text{sc}}$, for which all relative deviations drop below the threshold value, i.e.

[TABLE]

and call it a semiclassical-to-quantum boundary. Hence, for spin ensembles with $N>N_{\text{sc}}$ the semiclassical Maxwell-Bloch equations provide trustworthy results for the cavity probability amplitude $|\langle a_{st}\rangle|^{2}$ and all higher correlations like $\langle\sigma_{j}^{z}a\rangle_{c}$, $\langle\sigma_{j}^{-}a^{\dagger}\rangle_{c}$, etc., give only negligible relative contributions (lower than 1%). The vanishing influence of these quantum correlations, which is characterized by Eq. (21), captures a subtle but important point of open quantum systems. While the true quantum solution of the steady state master equation $\dot{\rho}=0$ is unique, the bistability of the semiclassical regime translates into a bimodality of the corresponding quasi-probability function. The quantum system switches between these two quasistationary components giving a unique time averaged expectation value which can substantially deviate from the semiclassical solution within the bistability region Dombi et al. (2013); Casteels et al. (2017); Fink et al. (2017). The switching time, however, strongly depends on the system size being divergent for large $N$ Casteels et al. (2017); Fink et al. (2017). Hence, the convergence criterion provided by Eq. (21) implies that for $N>N_{\text{sc}}$, the cavity probability amplitude $|\langle a_{st}\rangle|^{2}$ calculated from a full quantum mechanical evolution of the initial state can not be distinguished from the solutions of the semiclassical Maxwell-Bloch equations on experimentally feasible timescales.

Figure 5 shows this semiclassical-to-quantum boundary value as a function of both the cooperativity parameters $C$ and the driving strengths $\eta$. This main result of our paper demonstrates how the value of $N_{\text{sc}}$ increases close to the critical point $\eta^{+}_{\text{crit}}$. Exactly at $\eta=\eta^{+}_{\text{crit}}$ the time the systems needs to reach its stationary state for $C\geq 8$ diverges due to the effect of critical slowing down Bonifacio and Meystre (1979) and data points are therefore omitted for these parameters. For parameter regions where no bistability occurs, i.e. for $C<8$, there is no sharp distinction between the lower and upper transmission branch. Here the value of the semiclassical-to-quantum boundary $N_{\text{sc}}$ has its maximum at driving strengths slightly above the critical driving strength $\eta^{+}_{\text{crit}}$ and starts to peak at $\eta=\eta^{+}_{\text{crit}}$ only as $C$ approaches the threshold value of bistability. As can be seen in Fig. 5(b) for the cooperativity parameter $C=5$ the value of $N_{\text{sc}}$ above $\eta^{+}_{\text{crit}}$ is significantly larger than below that driving strength. This asymmetry with respect to $\eta^{+}_{\text{crit}}$ becomes less pronounced but is still present for increasing cooperativities up to $C\leq 8$. The peak in $N_{\text{sc}}$ for $C=7.8$ at the critical driving is the precursor of the emergence of a first-order phase transition and the effect of bistability, which emerges for cooperativity values above $C=8$.

It is worth noting that for $C>8$ the semiclassical-to-quantum boundary as defined in Eq. (21) behaves qualitatively differently for the lower transmission branch as compared to the upper transmission branch. As we approach the critical point $\eta^{+}_{\text{crit}}$ from below by analyzing stationary states disposed on the lower transmission branch, the value $N_{\text{sc}}$ progressively increases and eventually diverges exactly at $\eta^{+}_{\text{crit}}$, where the saddle-node bifurcation occurs [curves from the left with respect to the gray bar of Fig. 5(c)]. In contrast, the only stable solutions above $\eta^{+}_{\text{crit}}$ are those which are located on the upper transmission branch which lie, however, far away in phase space from the saddle-node bifurcation at $\eta^{+}_{\text{crit}}$. Therefore when approaching $\eta^{+}_{\text{crit}}$ from above, the value of $N_{\text{sc}}$ exhibits no divergence and is significantly smaller than below $\eta^{+}_{\text{crit}}$.

Another interesting tendency seen in Fig. 5(c) is that for $\eta<\eta^{+}_{\text{crit}}$ (lower transmission branch) $N_{\text{sc}}$ increases for increasing cooperativity parameters $C$. In contrast, for $\eta>\eta^{+}_{\text{crit}}$ (upper transmission branch) increasing cooperativity parameters $C$ lead to a decrease of $N_{\text{sc}}$. Our findings suggest that close to the critical point at the lower transmission branch even very large ensembles of up to $\sim 10^{4}$ spins can show behavior that goes beyond the semiclassical description. Note that for large cooperativities the differences in $N_{\text{sc}}$ for driving strengths within and outside the bistable region become very large. Whereas for $C=18$ and $\eta=1.01\,\eta^{+}_{\text{crit}}$ the semiclassical result agrees well (1% deviation) with the CE2 and CE3 already for ensembles of $N_{\text{sc}}\approx 500$ spins, the corresponding value grows to $N_{\text{sc}}\approx 3\times 10^{4}$ for $\eta=0.99\,\eta^{+}_{\text{crit}}$. This can be explained by very large quantum fluctuations near the critical point, which destabilize one of the two semiclassical basins of attraction. Note that the asymmetry between the lower and upper transmission branches encountered in Fig. 5(c) results from our choice to probe the critical point $\eta^{+}_{\text{crit}}$. We anticipate a reversed role of the two transmission branches when probing the stationary states of the upper transmission branch close to the critical point $\eta^{-}_{\text{crit}}$.

III.3 Inhomogeneous broadening

In the following we extend our investigations to inhomogeneously broadened spin ensembles and study how the broadening of the spin ensemble alters the previously defined minimal number of spins $N_{\text{sc}}$ that is required for the validity of the semiclassical Maxwell-Bloch equations. We therefore relax the condition $\Delta_{j}=0$ and allow for a Gaussian frequency distribution of the individual spins. This scenario becomes computationally much more demanding, since now the hierarchy of equations (summarized in the Appendix) has to be solved for each spin frequency $\Delta_{j}$ individually. Numerically we split the spin ensemble into $L$ equidistantly spaced frequency clusters $\Delta_{\mu}$, where the index $\mu$ runs from $1$ to $L$ and each frequency cluster is filled up with $M_{\mu}$ spins. After this procedure, the equations (A.1)-(A.25), representing the CE3, in total are $13L^{2}+L(L+1)/2+23L+9$ first-order ordinary differential equations, which can be solved for moderate values of $L$. For our calculations we chose $L=51$ and distributed the spins following a Gaussian distribution

[TABLE]

with $K=\sum_{\nu=1}^{L}e^{-4\ln(2)\frac{\Delta_{\nu}^{2}}{\Gamma^{2}}}$ being a normalization constant such that $\sum_{\mu=1}^{L}M_{\mu}=N$.

In Fig. 6 we present results for Gaussian spin distributions with three different full widths at half maximum, $\Gamma=0.1$, $0.5$, and $1.0$ MHz. Note that an increase of the width $\Gamma$ leads to a decrease of the collective cooperativity

[TABLE]

with all other parameters kept constant. This drop of the collective cooperativity for increasing widths of the distribution can be observed in Fig. 6(a), where we show the semiclassical bistability curves of the inhomogeneously broadened spin ensembles in comparison with the unbroadened case (corresponding to a collective cooperativity of $C=18$). The spin distributions with a full width at half maximum of $\Gamma=0.1$, $0.5$, and $1.0$ MHz, then corresponds to a collective cooperativity parameters of $C\approx 17.9$, $15.8$, and $12.7$, respectively.

The minimal number of spins, $N_{\text{sc}}$, for which the cumulant expansion converges towards the semiclassical results is presented in Fig. 6(b). A comparison with Fig. 5(c) indicates that the change in the semiclassical-to-quantum boundary due to the inhomogeneous broadening can be well understood in the way the collective cooperativity parameter $C$ changes with broadening. Our results for the inhomogeneously broadened spin distributions therefore confirm our earlier findings that even very large spin ensembles of about $10^{4}$ spins can show non-semiclassical behavior close to the critical point of bistability.

IV Conclusions

We have studied in detail the route towards the semiclassical limit for a dissipative spin-cavity system driven close to the critical point of amplitude bistability. In particular we analyzed the validity of the semiclassical Maxwell-Bloch equations close to the critical point following the transition from the lower to the upper transmission branch of amplitude bistability for varying cooperativities $C$ and for different numbers of spins $N$. We numerically solved the nonlinear sets of equations resulting from a second- (CE2) and third-order cumulant expansion (CE3) and compared the results with the semiclassical stationary solution for the cavity probability amplitude. Based on the convergence of the cumulant expansion towards the semiclassical results, we defined a criterion for the reliability of the Maxwell-Bloch equations and determined the minimal number of spins $N_{\text{sc}}$ necessary to ensure the validity of the semiclassical approximation.

Our results reveal that not only the distance to but also the way of approaching the critical point is a crucial factor which strongly influences the validity of the semiclassical equations. More specifically, we disclose that the large quantum fluctuations inside the bistable region lead to very large values of $N_{\text{sc}}$ in the proximity of the critical point of the lower transmission branch. Remarkably, here even very large spin ensembles of up to $\approx 10^{4}$ spins can feature deviations from the semiclassical cavity probability amplitude. Our results therefore suggest that a spin ensemble of the same size can behave semiclassically or quantum mechanically depending not only on the system parameters but also on the proximity of critical points and the way of approaching them.

Acknowledgements

We would like to thank H. Dhar, U. Hohenester, and P. Rabl for helpful discussions and the European Commission under Project No. NHQWAVE (MSCA-RISE 691209) for support. S.R. and M.Z. acknowledge support by the Austrian Science Fund (FWF) through Project No. F49-P10 (SFB NextLite) and the Doctoral Program CoQuS (W1210). Some of the computational results presented have been achieved using the Vienna Scientific Cluster (VSC).

Appendix: Hierarchy of coupled equations

In order to employ the CE2 and CE3 as described in Sec. II.3 the equations of motions have to be generated also for second- and third-order expectation values. The derivation of these equations is straightforward but the arising expressions soon become unwieldy for higher orders of expectation values. Since the resulting equations, which enable an accurate description of spin-cavity systems including the effect of inhomogeneous broadening, constitute a critical part of our paper, we provide them explicitly below:

First-order expectation values:

[TABLE]

[TABLE]

[TABLE]

Second-order expectation values:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Third-order expectation values:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]