# Quantum sensing close to a dissipative phase transition: symmetry breaking and criticality as metrological resources

## Abstract

We study the performance of a single qubit-laser as a quantum sensor to measure the amplitude and phase of a driving field. By using parameter estimation theory we show that certain suitable field quadratures are optimal observables in the lasing phase. The quantum Fisher information scales linearly with the number of bosons and thus the precision can be enhanced by increasing the incoherent pumping acting on the qubit. If we restrict ourselves to measurements of the boson number observable, then the optimal operating point is the critical point of the lasing phase transition. Our results point out to an intimate connection between symmetry breaking, dissipative phase transitions and efficient parameter estimation.

###### pacs:
03.67.Ac, 37.10.Ty, 75.10.Jm, 64.70.Tg

## I Introduction

Quantum sensing and metrology is likely to be a key practical application of quantum technologies. It has been established both by theory and experiments that quantum effects can be exploited to increase the accuracy of measurement devices caves81prd; braunstein94prl; sanders95prl; bollinger96; giovannetti11natphys. Practical applications, however, face significant challenges. In an ideal scenario quantum metrology requires the preparation of many-particle entangled states by quantum operations that so far are only possible with a few degrees of freedom. Dissipation and noise pose severe limitations which often hinder the metrological advantages of entangled states huelga97prl; escher11; guta12natcomm; chin12prl; huelga16prl. Quantum setups such as superconducting circuits majer07; houck12natphys; devoret13sci and trapped ions Schneider12rpp; blatt12natphys offer us the opportunity to engineer quantum states of matter with a high degree of control over interactions and dissipation. It has been shown that dissipation may be actually exploited as an effective tool in quantum state engineering verstraete09; barreiro11. The question naturally arises, whether we can use dissipation to design metrological protocols and sensors dallatore13prl; gonzaleztudela13prl. We propose two working principles for such quantum sensors. First, one could exploit the sensitivity of a dissipative steady-state to an external field which explicitly breaks some suitable underlying symmetry. The second route could take advantage of the sensitivity at the critical point of a dissipative phase transitionZanardi08; wang14njp; salvatori14pra; guta16pra. Such a sensor would have the advantage that state preparation is not required and, furthermore, dissipation is a control parameter of the sensor dynamics, rather than an error source.

In this work, we show that a single qubit laser is a minimalist model where both ideas can be tested. A macroscopic laser with photons can be described by a coherent state of the light field with a mean value , which assumes the spontaneous breaking of the underlying lasing phase symmetry by choosing an arbitrary value of molmer97pra. This approach can be justified by assuming an infinitesimal field (e.g. an environmental fluctuation) that fixes the laser phase scully70pra. However, in a finite size system (e.g. a single qubit laser) an external field with finite amplitude, , is required to explicitly break the phase symmetry (see Fig.1). Here the thermodynamic limit is found when hwang15, at which the system undergoes a spontaneous symmetry breaking, i.e. . This relation implies that the order parameter must increase with the system size with an scaling yet to be determined, leading to a high sensitivity to .

Our article is structured as follows. Firstly, we present a semi-classical description in phase space of a single-qubit laser in the presence of a weak symmetry breaking driving field. This allows us to estimate analytically the quantum Fisher information related to the amplitude and phase of the driving, which further shows the connection between symmetry breaking and efficient parameter estimation. We identify the optimal observables that fully exploit the system metrological capacity. Non-equilibrium criticality is then examined as an alternative metrological resource with non-optimal protocols using the average number of bosons. We conclude with a discussion of possible error sources as well as applications.

## Ii Single qubit laser

In this work we consider a bosonic mode coupled by a Jaynes-Cummings interaction to a two-level system (qubit) with levels and . Additionally, we introduce a periodic driving which becomes the target weak field. Both the qubit and driving frequencies are assumed to be resonant with the bosonic mode. In an interaction picture rotating at the mode frequency, the coherent dynamics is described by the Hamiltonian

 H = HJC+Hd, HJC = g(σ+a+a†σ−),Hd=ϵ∗a+ϵa†, (1)

where , with and being the driving amplitude and phase, respectively. are the ladder operators of the two-level system, and . In addition to this coherent dynamics, the system is subjected to incoherent pumping of the qubit and losses of the bosonic mode with rates and , respectively. The resulting dissipative process is well captured by the following master equation for the system density matrix ,

 ˙ρ=−i[H,ρ]+L{σ+,γ}(ρ)+L{a,κ}(ρ), (2)

where Lindblad super-operators are defined as . In a mean field approximation to the case without driving () the steady-state is determined by the pump parameter, . This sets a dissipative phase transition into a lasing phase when BreuerPet, being the order parameter. To evaluate the response of the single-qubit laser to an external driving, we need to go beyond mean-field theory. Since we are only interested in the output laser field, we start by finding an effective Liouvillian able to describe the reduced dynamics of the bosonic mode. This can be accomplished in a strong pumping regime Sargent74; Mandel95, i.e , in which the qubit can be adiabatically eliminated, leading to an effective quartic master equation for the bosonic mode (see App.A for a detailed derivation),

 ˙ρf =−i[ϵ∗a+ϵa†,ρf]+L{a†,A}(ρf)+L{a,C}(ρf)+ +L{aa†,B}(ρf)−L{a†2,B}(ρf). (3)

We have defined the coefficients , , , and is the reduced density matrix of the bosonic field. Our expression is valid in a regime of strong incoherent pumping, such that the probability of occupation of the ground state can be neglected. This condition is justified both below the lasing phase transition, , or slightly above the threshold, (see App.A).

## Iii Semi-classical limit

The master equation obtained in (II) is still challenging to be tackled analytically. By using phase space methods we shall obtain a Fokker-Planck equation valid in a regime with high number of bosons Gardiner; Haken75. This will allow us to get analytical results that will be assessed below by comparing to exact numerical calculations. We start by introducing the coherent state or Glauber-Sudarshan P representation of the effective master equation Carmichael; Mandel95, defined as

 ρ(t)=∫d2αP(α,α∗,t)|α⟩⟨α|, (4)

where is the coherent state . The function plays a role analogous to that of a classical probability distribution over , with the normalization condition , and expectation values of normal ordered operators, . Note that is actually a quasi-probability distribution, since it is in general not a positive distribution function.

By substituting the representation (4) of into Eq. (II), one may convert the operator master equation into an equation of motion for . This can be accomplished by using the following equivalences,

 a|α⟩⟨α| =α|α⟩⟨α| (5) |α⟩⟨α|a† =α∗|α⟩⟨α| (6) a†|α⟩⟨α| =(∂∂α+α∗)|α⟩⟨α| (7) |α⟩⟨α|a =(∂∂α∗+α)|α⟩⟨α|. (8)

An integration by parts with the assumption of zero boundary conditions at infinity, which introduces an extra minus sign for each differential operator, converts the integrand of (4) into a product of and a c-number function of . This leads to a differential equation for . When the laser is operating near the steady state and above threshold, is a large number of the order of the average number of photons. Notice also that is a very small coefficient compared to , such that . Consequently, we shall retain only the most important terms in . This corresponds to dropping any contribution smaller than . In doing so, we end up with the following Fokker-Planck equation for ,

 ∂P∂t=−∂∂α[(A−C−B|α|2)α−ϵ′]P+c.c.+2A∂2P∂α∂α∗ (9)

where . Let us write this equation in cartesian coordinates, with and then

 ∂P∂t=−2∑i=1∂∂xi[(A−C−B→x2)xi−ϵ′i]P+A22∑i=1∂2P∂x2i (10)

where we introduce the two-dimensional vectors and . In the stationary state , Eq.(10) may be rewritten as , where the current is defined by

 Ji=[(A−C−B→x2)xi−ϵ′i]−A2∂P∂xi. (11)

When the drift vector satisfies the potential condition , as it does in our case, the solution to the Fokker-Planck equation is derived by imposing Mandel95. This leads to a differential equation for that can be directly integrated to give

 P(→x)=1Nexp{1A[(A−C−B2→x2)→x2−2→ϵ′⋅→x]}, (12)

where is a normalization constant. The steady-state solution (12) can be conveniently expressed in polar coordinates as follows,

 P(r,θ)=1Nexp(−λr4+μr2−2νrsin(θ−ϕ)), (13)

where we have introduced the parameters , , and . Note that the probability distribution (13) is positive, which indicates that the steady-state admits a classical description. Eq. (13) can be used to calculate expectation values in the steady state through parametric derivatives of the normalization constant, . In the App.B, it is detailed how to approximately calculate using the Laplace’s method together with explicit expressions of useful observables. In the absence of driving (), the laser phase is uniformly distributed in , implying that any average field quadrature vanishes. In contrast, when , the driving field explicitly breaks the phase symmetry and the state adopts a preferred phase with exponential sensitivity as illustrated in Fig.1. According to Eq.(13) we expect the output laser field to have a phase delay of with respect to the input driving.

Not any explicit symmetry breaking may lead to an advantageous sensing scheme. However, when this is associated to a spontaneous symmetry breaking (SSB) in the thermodynamic limit, the corresponding order parameter is expected to be very sensitive to such symmetry breaking field. Such subclass of non-trivial explicit symmetry breaking process is henceforth referred as induced symmetry breaking. In our case, this general symmetry argument is translated as a high sensitivity of the coherent component to . The average field quadrature will be shown to be particularly sensitive to the external driving, analytically given by

 ⟨^Pϕ⟩=2r0I1(2νr0)I0(2νr0)≈νr0≪12r20A|ϵ|, (14)

where are the modified Bessel functions of the first kind, and stands for steady average number of bosons with no driving (see App.B for detailed derivation),

 ⟨n⟩ϵ≈0=r20=(A−C)/B. (15)

The SSB of the lasing phase transition here implies . This entails a certain scaling of with the system size, here , now explicitly given by (14). Figure 2 shows the comparison of these results with numerical calculations of the exact and the adiabatic equation, Eqs.(2) and (II) respectively. From Fig.2 we differentiate two distinct regimes. First, a linear regime if is small enough, where scales linearly with the number of bosons. Essentially, the more pronounced the slope is, the higher the sensor sensitivity will be. Second, a saturation regime where and the laser admits a fully classical description Mandel95.

## Iv Quantum Fisher information and optimal measurements

Eq. (14) suggests that the induced symmetry breaking allows us to measure weak field amplitudes . The capability of this sensing scheme will be mainly determined by its resolution. The theory of quantum Fisher information Dowling15; Demkow14 provides us with an ultimate lower bound on the precision of parameter estimation that is possible in a quantum model, which will be used to assess the maximum metrological capacity of the single-qubit laser.

Assume that a target parameter is encoded in a certain density matrix . The quantum Cramer-Rao bound establishes a lower bound to the error in the estimation of ,

 Δ2φ≥1NexpFQ[ρφ] (16)

where is the quantum Fisher information (QFI) and the number of experiment repetitions. The QFI can be viewed as a quantitative measure of distinguishability of a state from its neighbors . Thus it can be used as a quantitative characterization of the maximal sensor resolution. A measurement scheme that saturates the bound Eq.(16) is called optimal. The symmetric logarithmic derivative operation (SLD) is known to be optimal for all quantum states Braunstein96. It is defined by the Hermitian operator satisfying the relation

 ∂φρφ=12(ρφLφ+Lφρφ). (17)

The QFI is then given by . In the eigenbasis of , the SLD is written as

 Lφ[ρφ]=∑i,jλi+λj≠02⟨ei(φ)|˙ρφ|ej(φ)⟩λi(φ)+λj(φ)|ei(φ)⟩⟨ej(φ)|, (18)

First we shall focus on the estimation of the field amplitude for a given known phase . By using the analytical result for the steady state (13), we aim for deriving theoretical results for the SLD as well as the QFI. To do so, it is necessary to solve the operator equation (17) for . In this context, a comprehensive solution of Eq.(17) is already known for Gaussian states in phase space, i.e quadratic in Monras13. Assuming the adiabatic elimination regime, i.e , the coefficients satisfy . Hence the function (13) can be well approximated by the following Gaussian-like approximation,

 P(r,θ)=N−1exp(−(r−r0)22σ2−νrsin(θ−ϕ)) (19)

where and . Even though this represents a simplification with respect to the original function (13), the state is still not Gaussian in the variables , for which exact solutions are known for the SLD and QFI Monras13. Even so, let us try to solve the equation (17) in the coherent state representation. Using (19), the l.h.s of the equation (17) gives

 ∂|ϵ|P(r,θ)=(−N−1∂|ϵ|N+iA(αe−iϕ−α∗eiϕ))P. (20)

It turns out that is equivalent to the average of the field quadrature . This result induces us to introduce the ansatz , with proper coefficients, which corresponds essentially to the measurement of a suitable field quadrature. Inserting this ansatz in the r.h.s of (17) and bearing in mind the equivalences (7,8), we have

 L|ϵ|ρ=∫∞0∫2π0rdθdr(S0+Sα+S∗(α∗−∂α))P (21)

with analogous result for . In a deep lasing regime (well above threshold but still within the validity regime of (II)) where , the derivative in Eq.(21) can be simplified assuming that , yielding

 ∂αP=e−iθ2(∂∂r−ir∂∂θ)P==(−α∗2σ2+r02σ2e−iθ+i|ϵ|Ae−iϕ)P≈i|ϵ|Ae−iϕP. (22)

Identifying now terms from both sides of the equation (17), the SLD reads

 Lϵ[ρ|ϵ|]=1A(−⟨^Pϕ⟩+|ϵ|A+^Pϕ). (23)

The contribution can be neglected in comparison with the contribution given by , leading to the SLD . Happily, this in turn implies that , a property that any SLD must fulfill according to its own definition (17). The QFI may now be calculated as in terms of a parametric derivative of the normalization constant introduced in (13), specifically as the fluctuations of (see App.B),

 FQ[ρ|ϵ|]=2r20A2⎛⎝1+I2(2νr0)I0(2νr0)−2(I1(2νr0)I0(2νr0))2⎞⎠≈νr0≪12r20A2. (24)

In Fig. 3 we show a comparison between the analytical result (24) and an exact numerical calculation of Eq.(2) by using (18). There are two important conclusions that are drawn from (24). Firstly, it shows that the metrological capacity for estimating is maximal when the induced symmetry breaking occurs, and decreases as the symmetry is already broken. This is intuitively natural since the parameter is directly associated with the symmetry breaking, and the gain of information is maximal at that point. This feature can be reasonably expected in any sensing scheme relying on spontaneous symmetry breaking as this one. Consequently, this type of sensing is advantageous when measuring extremely weak fields as the precision naturally increases in such domain. The parameters of the laser can be adjusted so that the amplitude remains in the first order approximation, where the precision remains constant for a fixed amplitude as Eq.(24) indicates. Secondly, the QFI scales linearly with the steady average number of bosons as . In the macroscopic limit, defined here as , diverges as a result of the sensitivity of the steady-state to an infinitesimal perturbation, giving rise to a spontaneous symmetry breaking. These results show a useful connection between symmetry breaking and efficient parameter estimation. The prior knowledge of in estimating may be eluded by performing an average of different quadratures over the range , decreasing the QFI by a factor but still conserving the same scaling.

In the light of these results, we examine whether a similar approach can be used for measuring the phase for a given amplitude. A completely analogous procedure can now be used to solve again the operator equation (17). Now the l.h.s of the equation (17) gives

 ∂ϕP(r,θ)=(−N−1∂ϕN+2ν(αe−iϕ+α∗eiϕ))P. (25)

The term is easily shown to be zero. Using then a linear ansatz , the r.h.s of (17) is analogous to Eq.(21). The comparison between both sides of the equation yields the SLD,

 Lϕ[ρϕ]=ν(^Xϕ). (26)

where is the field quadrature . The operator also satisfies as required by the definition (17) . The QFI is then , which turns out to be equivalent to

 FQ[ρϕ]=ν⟨^Pϕ⟩≈νr0≪12r20|ϵ|2A2, (27)

where we have used Eq.(14). This result predicts that the QFI scales linearly with and quadratically with the field amplitude as . Graphically, the behavior of the QFI in this case is indirectly given in Fig.2. In contrast to Eq.(24) for estimating the amplitude, here the QFI increases with since naturally a non-zero signal is required to have a localized phase. Note that the optimal observable depends itself on the target parameter, . To operate in the optimal measurement regime we need a first estimation of the observable, . If such estimation satisfies the condition , the quadrature leads to an optimal protocol for estimating , with a precision determined by Eq. (27). This requirement is analogous to the optimal free precession time in Ramsey spectroscopy Ludlow15.

In summary, our optimal scheme makes use of the coherent component to estimate within the linear regime of induced symmetry breaking, being and the optimal observables for estimating the amplitude and phase respectively. We stress the fact that the quantity appearing in Eqs.(24,27) refers to the number of bosons in the steady state, whose main contribution comes from the incoherent pumping but not the probe field, concretely in the lowest order. This implies that one can increase the precision in parameter estimation for a fixed driving intensity solely by increasing the laser pumping . Additionally, recall that none of the results presented in this work depend on the quantum state of the driving field, as the system steady-state is unique for all of them.

## V Criticality as a metrological resource

The results obtained in Eqs. (24, 27) constitute the maximal metrological capacity of the single qubit laser for estimating and respectively, as they saturate the Cramer-Rao bound (16). However, estimation by non-optimal observables, like the steady number of bosons , may be experimentally more convenient. Using the analytical results for and (see App.B), the expected relative error above threshold for estimating by means of is

 Δ|ϵ||ϵ|=1|ϵ|Δn∂n∂|ϵ|=Cpκgν2+g2γ(Cp−1Cp)+O(|ν|2). (28)

Eq.(28) indicates that the precision increases as we approach the critical point , at which the precision scales as 1. The ratio between the optimal and non-optimal protocols, suggests that both methods give comparable resolutions when . Figure 4 depicts exact numerical results for , confirming maximal precision around the critical point as the thermodynamic limit is approached. Such limit is reached when hwang15, or equivalently .

The maximal precision given by the critical point manifests a connection between non-equilibrium criticality in dissipative systems and efficient parameter estimation. An analogous result has been already explored for closed systems Zanardi08. Physically, it is intuitive to think that the system at the critical point becomes more sensitive to any perturbation, leading to a greater sensor resolution. The potential of criticality for sensing can be exploited in setups where the qubit-boson coupling, , can be controlled with the necessary accuracy to ensure that the system stays at the critical point. This is actually the case in, e.g., single trapped ion phonon lasers, where this coupling is implemented by a laser and its strength modulated by its intensity. Also, in superconducting qubits, qubit-photon coupling terms can be induced and controlled with periodic driving fields navarrete14prl.

A phase estimation by measuring the number of bosons is also possible if we extend the previous setup to arrange an adequate interferometric scheme. Concretely, we add a new reference field term to Eq. (II), where we assume that , are known parameters. Both the probe field and the reference field must be comparable to observe interference effects, so we shall assume for simplicity that they both have the same amplitude, . One may treat this new input field as we did in the previous sections, in which case the function for the steady state will be

 P(r,θ)=1Nexp(−λr4+μr2−2ν′rsin(θ−ϕ′)) (29)

with and . Comparing Eq.(29) with Eq.(13), we note that the addition of the reference field to the probe field leads to a total driving field with phase and amplitude . Interference has thus translated the information of into a new phase-dependent amplitude , which can be now estimated through measurements of the average boson number with the precision shown in Eq. (28). In the lowest order this leads to a precision , showing that the optimal operating condition is .

## Vi Possible sources of errors

One may wonder whether potential sources of error in real experiments could jeopardize our previous results. In the App.C we consider three possible sources of error: dephasing of the qubit, heating of the bosonic mode and detuning between the qubit and the mode. Our calculations show that the detuning is expected to be negligible as long as , while the dephasing and heating results in a renormalization of the constants and . In essence, our results are robust to any perturbation that respects the symmetry of the model and the universal scalings of the lasing phase transition.

## Vii Physical Implementations

Single-qubit photon lasers can be implemented with single atoms mckeever03 or superconducting qubits Astafiev07; Hauss08; navarrete14prl. Furthermore, our ideas can be also applied to single-qubit phonon lasers vahala09; Grudinin10prl. Here, the quantized excitations (phonons) of a trapped ion play the role of the photons in an optical laser, whereas internal electronic levels provide us with a qubit. Our scheme would lead to a the precise measurement of ultra-weak forces of the form resonant with the trapping frequency biercuk10natnano; schreppler14sci; ivanov15prapp; ozeri16arx. Phonon lasing has actually been already observed in a single trapped ion experiment vahala09. All the interactions and techniques required to implement this idea are routinely used in trapped ion experiments, see for example Leibfried03rmp for an excellent review on the topic.

To have full control of the parameters involved in our model we will consider a two-ion crystal in which one of the ions acts as a single-atom phonon laser, whereas a second auxiliary ion is used to provide us with a sympathetic cooling mechanism Barret03pra. To avoid the requirement of individual addressing of each of the two ions, different species could be used. We assume that ions are weakly coupled by the Coulomb interaction. We introduce phonon annihilation operators and associated to quantized vibrations of ions 1 and 2 respectively. The coupling term between the ions takes the form Porras04prl2; Haze12pra,

 Hc=tC(a†1a2+a1a†2). (30)

If we consider radial vibrations, then , where is the distance between ions, refers to the ion’s mass and is the trapping frequency.

Let us consider now the first ion’s quantum dynamics. To make the connection with trapped ion physics clearer, we will work in a spin basis where the role of states and is interchanged with respect to the discussion in the main text. In our trapped ion scheme, spin pumping will be induced by the radiative decay from an excited state to the ground state , whereas a spin-phonon coupling of the form will be induced. This is described by the following Liouvillian,

 L1(ρ)=−i[H1,ρ]+L{σ−1,γ}(ρ). (31)

The Hamiltonian acting on ion includes a blue-sideband coupling between the internal state of the ion and the local vibrational mode as well as the coupling to the external force that we aim to measure,

 H1=g(σ+1a†1+σ−1a1)+ϵ(a†1+a1). (32)

We have introduced ladder operators, , , associated to the internal state of ion 1. The blue side-band term can be induced by lasers with frequency , where is the frequency of the internal state transition Leibfried03rmp. Finally, the last term of Eq (31) is simply the radiative decay of the excited state Leibfried03rmp. To ensure that the dynamics of the ion is constrained to only two levels, one could simply choose and as the two levels of a cycling transition.

The only missing element is a cooling mechanism acting on ion 1. For this we will use ion 2 to provide us with a cooling medium by an effect known as sympathetic cooling. For this we assume that ion 2 is being continuously laser cooled with a rate ,

 L2=L{a2,κ}(ρ). (33)

If the Coulomb coupling is small relative to the cooling rate (), we can adiabatically eliminate ion 2 and obtain an effective cooling term for ion 1, with cooling rate . The reduced density matrix for ion 1, , is thus subjected to the following quantum dynamics,

 ˙ρ1=L1(ρ1)+L{κeff,a1}(ρ1). (34)

Our scheme is a phononic version of the single-qubit laser described in the main text. To assess the sensitivity of such a device in the measurement of external forces, we consider now some typical values for cooling rates and vibrational couplings. We focus on the optimal measurement protocol, which would imply measuring the quadrature, , defined above Eq. (6) of the main text. Quadratures of vibrational operators can be efficiently measured by coupling phonon observables to the ion’s internal state and detecting the emitted fluorescence (see for example Leibfried03rmp). By using our calculation of the error as estimated from the QFI we get

 Δϵ=1√FQ[ρ|ϵ|]=A√2r0. (35)

To estimate , we express it like , where we have assumed that we work in a regime with cooperativity parameter .

Our scheme can be applied to measure ultra-weak forces. The relation between the driving strength and the applied external force, , is , where

 x0=1√2mωT, (36)

is the size of the vibrational ground state. Our final expression for the force sensitivity reads (in standard units including ),

 ΔF≈ℏκeff√2nphx0, (37)

where we have used the fact that the number of phonons, . To get an estimate of the precision with which an ultra-weak force could be measured, we consider that ion 1 is Ca and 10 MHz, which yields 3.5 nm. Other typical values are kHz Haze12pra and kHz, leading to kHz. With those values we get

 ΔF≈53 yN/√nph, (38)

By increasing the number of phonons in the lasing regime to values such as , one could obtain precisions yN, well within the yocto-Newton regime and beyond the precision of results reported in experiments biercuk10natnano.

Large phonon numbers are in principle not difficult to get in a trapped ion phonon laser. For example, taking into account typical values of MHz our Eq.(15) yields the value with a side-band coupling kHz, well within the state-of-the-art Leibfried03rmp.

A limiting factor could be the presence of motional heating, . However, heating rates in linear Paul traps can be as low as 0.1 vibrational quanta per ms, which translates into kHz Turchette00pra. Under those conditions, and the effect of heating could be neglected or incorporated into minor corrections to the trapped ion sensor (see section VI).

## Viii Acknowledgments

Funded by the People Programme (Marie Curie Actions) of the EUâs Seventh Framework Programme under REA Grant Agreement no: PCIG14-GA-2013-630955. We thank Jacob Dunningham and Pedro Nevado for fruitful discussions.

Here we shall derive the effective quartic master equation claimed in (II), as a result of the adiabatic elimination of the fast spin variable. Firstly, we shall trace over the spin degree of freedom from the master equation for the single-qubit laser,

 ˙ρ=−i[H,ρ]+L{σ+,γ}(ρ)+L{a,κ}(ρ), (39)

thereby obtaining an equation for the reduced density matrix of the bosonic field . Namely, this equation reads

 ˙ρf=−ig(aρge+a†ρeg−ρgea−ρega†)−−i(ϵa†ρf+ϵ∗aρf−ϵρfa†−ϵ∗ρfa†)++κ(2aρfa†−a†aρf−ρfa†a), (40)

where we introduced the notation and . To obtain a closed equation for the reduced density matrix , we have to eliminate the operators from Eq. (40). We obtain the corresponding equations of motion for these operators using the original master equation,

 ˙ρge=−ig(a†ρee−ρgga†)−γρge, (41)

where we have neglected the contributions from and in comparison with . In the limit , we can adiabatically eliminate the operators and from (40) by taking in Eq.(41) and substituting their steady-state solutions,

 ρge=−igγ(a†ρee−ρgga†). (42)

As the resulting equation still depends on the operators and , we make use of the single-qubit master equation to obtain the equations of motions of these operators,

 ˙ρee =−ig(aρge−ρega†)+2γρgg (43) ˙ρgg =−ig(a†ρeg−ρgea)−2γρgg (44)

where we again neglect terms with and . One may now obtain a perturbative solution to the steady-states of Eqs.(43)(44) in terms of the field density matrix . To do so, let us adiabatically eliminate by taking in Eq.(44), yielding

 ρgg=−ig2γ(a†ρeg−ρgea)=g22γ2(2a†ρeea−a†aρgg−ρgga†a). (45)

In a first order approximation, the ground state population is negligible due to the fast pumping of the atoms (). Therefore, we expect to find and in first order. A second order correction is achieved by inserting this first order approximation into Eq.(45), hence

 ρgg =g2γ2a†ρfa (46) ρee =ρf−ρgg=ρf−g2γ2a†ρfa. (47)

One can finally insert Eqs.(46)(47) into Eq.(40) to arrive at the desired closed equation for ,

 ˙ρf=−i(ϵa†ρf+ϵ∗aρf−ϵρfa†−ϵ∗ρfa)++g2γ(2a†ρfa−aa†ρf−ρfaa†)++2g4γ3(aa†ρfaa†−a†2ρfa2)++κ(2aρfa†−a†aρf−ρfa†a). (48)

The second term in the r.h.s. of Eq.(48) accounts for the single photon emission by the excited qubit (linear gain), while the third represents the contribution of two cycles of emission and re-excitation (gain saturation). Eq.(48) can be cast in Lindblad form as presented Eq.(II). A few brief remarks are worth mentioning about the single-qubit laser physics. Using Eq.(48) and setting , we can easily derive an equation for the diagonal elements , namely

 ˙ρnn=−(2A−B(n+1))(n+1)ρnn+2Anρn−1,n−1−Bn(n−1)ρn−2,n−2−2Cnρnn+2C(n+1)ρn+1,n+1, (49)

where we defined the coefficients , and . In contrast to the classic Scully-Lamb treatment of the four-level laserSargent74, no detail balance solution can be found to Eq.(49). The rate equation for the average photon number can also be derived from (49),

 ⟨˙n⟩=2(A−C)⟨n⟩+2A−B(2⟨n2⟩+5⟨n⟩+5). (50)

According to Eq. (50), there will be an initial exponential increase in the mean photon number if , hence is the threshold condition for the laser phase. This agrees with the prediction of a mean field treatment to this problem, in which the lasing phase is found when the pump parameter satisfies BreuerPet.
We now address the conditions of validity of the adiabatic elimination. Using Eq. (46), the condition is translated into

 ⟨σ+σ−⟩=Tr{ρgg}=(gγ)2(1+n)≪1. (51)

Below threshold (), this is satisfied as long as . Above threshold (), we can estimate (see App. B), which leads to

 ⟨σ+σ−⟩≈12Cp−1Cp≪1. (52)

Consequently, for the adiabatic elimination to be self-consistent above threshold we must require .

## Appendix B Laplace’s method

We shall calculate different observables associated to the function obtained in (13), corresponding to the steady-state solution of the Fokker-Planck equation (9). To do so, it is first necessary to compute the normalization constant given in (13). This can be approximately integrated using the Laplace’s method, which is helpful for integrals of the form

 I(s)=∫∞−∞f(x)esg(x)dx≈√2πsg′′(x0)f(x0)esg(x0), (53)

in which stands for the global maximum of , represents its second derivative evaluated at , varies slowly around and is independent of the parameter . In our case, has the form

 N=∫∞0∫2π0rdθdre(−λr4+μr2−2νrsin(θ−ϕ)). (54)

Integrating over gives

 N=2π∫∞0rdrI0(2νr)e(−λr4+μr2). (55)

where are the modified Bessel functions of the first kind. Above threshold, where , the normalization constant is approximated by Eq.(53) as

 N≈√π3λI0(2νr0)exp(μ24λ), (56)

where . The laser field quadrature may now be computed by taking the parametric derivative , which gives

 ⟨^Pϕ⟩=2r0I1(2νr0)I0(2νr0). (57)

Assuming that , one may expand (57) in Taylor series as

 ⟨^Pϕ⟩≈2νr20−ν3r40+O(ν4), (58)

which in first order indicates a linear dependence in as claimed Eq.(14). To compute the uncertainty of a second derivative is required, specifically , the result of which reads

 Δ2^Pϕ=2r20⎛⎝1+I2(2νr0)I0(2νr0)−2(I1(2νr0)I0(2νr0))2⎞⎠. (59)

When , a Taylor expansion of (59) gives

 Δ2^Pϕ=2r20(1−3ν2r202+O(ν4)). (60)

On the other hand, the parametric derivatives with respect to can be related to the average number of bosons and its uncertainty. First, the average number of bosons is given by