{'fraction_non_alphanumeric': 0.03992974238875878, 'fraction_numerical': 0.009821428571428571, 'mean_word_length': 4.466992078098744, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Keldysh meets Lindblad: Correlated Gain and Loss in Higher-Order Perturbation Theory

Motivated by correlated decay processes producing gain, loss and lasing in driven semiconductor quantum-dots, we develop a theoretical technique using Keldysh diagrammatic perturbation theory to derive a Lindblad master equation that goes beyond the usual second order perturbation theory. We demonstrate the method on the driven dissipative Rabi model, including terms up to fourth order in the interaction between the qubit and both the resonator and environment. This results in a large class of Lindblad dissipators and associated rates which go beyond the terms that have previously been proposed to describe similar systems. All of the additional terms contribute to the system behaviour at the same order of perturbation theory. We then apply these results to analyse the phonon-assisted steady-state gain of a microwave field driving a double quantum-dot in a resonator. We show that resonator gain and loss are substantially affected by dephasing-assisted dissipative processes in the quantum-dot system. These additional processes, which go beyond recently proposed polaronic theories, are in good quantitative agreement with experimental observations.

A notable experimental example of this, which motivates our work, is an electronically open, DQD system coupled to a driven resonator, as illustrated in Fig. 1(a),(b) [1]. Substantial gain in the resonator field was observed when the DQD is blue-detuned with respect to the resonator, and capacitively biased to induce substantial population inversion. The observed gain is attributed to correlated emission of a resonator photon and a phonon into the semiconductor medium in which the DQD system is defined. This process ensures conservation of energy, since the phonon carries the energy difference, (ω q − ω r ), between the energy of the qubit and the energy of the resonator phonon.
In some experimental regimes of [1], the observed gain is well described by a theory based on a canonical transformation to a polaron frame [3]. In this frame, conventional quantum optics techniques and approximations (Born-Markov, secular etc) are used to derive dissipative Lindblad superoperators that are quadratic in both the qubit-phonon bath coupling strength, β j , and in the * c.muller2@uq.edu.au † stace@physics.uq.edu.au qubit-resonator coupling strength, g. However the same theory fails to describe substantial loss (sub-unity gain) in other experimental regimes, which strongly suggests that there are additional dissipative processes that are not captured in the polaron frame. One problem with relying on a canonical transformation as the basis for a perturbative expansion is that it is tailored to a specific frame, which emphasises some processes over others. It is therefore not guaranteed to find all dissipative processes that occur at a given order in perturbation theory.
In this paper we present a systematic approach to derive a Lindblad-type master equation at higher perturbative orders, which contains all relevant correlated dissipative processes at a given order. We use the Keldysh real-time diagrammatic technique to calculate the perturbative series, and explicitly evaluate all fourth order diagrams that generate correlated dissipation. In doing so, we make explicit the link between the Keldysh selfenergy and the Lindblad superoperators. We then apply this technique to the experimental situation described above, showing that the additional terms we find quantitatively explain anomalous gain and loss profiles in driven DQD-resonator systems.
This paper is structured as follows: Section II introduces the model system that motivates this work, namely a single artificial two-level system (qubit) coupled to a resonator as well as a dissipative environment [1,2]. Section III derives the Keldysh self-energy superoperator in terms of irreducible Keldysh diagrams, and decomposing it into a sum of Lindblad superoperators to form the Keldysh-Lindblad master equation. The complete set of fourth-order Keldysh-Lindblad dissipators, presented in Eq. (20), is the central theory result of this paper. Section IV introduces resonator driving, and solves the steady-state Keldysh-Lindblad master equation for a driven, DQD-resonator system in a mean-field approximation.
Section V shows that the additional Lindblad dissipators that arise in the Keldysh analysis of the Dyson series quantitatively account for the of gain and loss in a driven DQD-resonator system, across all experimental regimes of Ref. 1. We finish in Section VI with a discussion of our results, specifically the correlated decay rates we obtain, that describe the dynamics of qubit and resonator.
The Keldysh and Lindblad formalisms are each well used in their respective research communities, however the link between them is usually not made explicit. Equations of Lindblad type are the generators of completely-positive, trace-preserving dynamics [17] and a wide variety of techniques have been developed to analyse the dynamics of open quantum systems under the actions of Lindblad master equations, including the inputoutput [18] and cascaded 'SLH' formalisms [19], quantum trajectories and measurement [20][21][22][23][24], and phasespace methods [25]. There are widely used, but somewhat heuristic techniques in quantum optics to derive Lindblad master equations at second order in the bath coupling [25]. Going beyond this approximation requires a well formulated perturbative approach.
The Keldysh diagrammatic technique [26] is widely used in condensed matter physics as a tool to calculate Greens functions of non-equilibrium systems [27]. It has been previously applied with great success for calculating the properties of mesoscopic quantum systems [14,[28][29][30][31]. However in these previous cases the resulting master equations are not typically expressed explicitly in the Lindblad form. Furthermore, these derivation are either limited to second order or include only a small class of diagrams for which closed form expressions are available.
In contrast here we develop a systematic approach to account for all Keldysh diagrams at a particular order, with which to derive a master equation that is explicitly in Lindblad form.
This paper therefore has the subsidiary objective of making explicit the connection between the Keldysh and Lindblad approaches to open quantum systems. As such, we have provided extensive appendices giving the technical details of our calculations. As we show, Keldysh perturbation theory is a powerful approach to deriving consistent Lindblad dissipators at arbitrary order. It is especially well suited for the treatment of open quantum systems in which correlated dissipative process, which go beyond the usual second order in perturbation theory, are significant. As such, in addition to the specific problem of correlated gain and loss in a driven DQD-resonator, we foresee immediate applications of this technique in numerous context, like co-tunneling in open quantumdots [32], charge transport through chains of Josephson junctions [33] or the calculation of correlated decay rates in multi-level circuit-QED systems [34]. An electron (black dot) tunnels from the source lead to the DQD state |L , which couples to the DQD state |R with matrix element ∆q, and then tunnels to the drain, leaving the DQD in the empty state |Ø . Also shown are inter-dot bias, q , and coupling rates ΓL,R to metallic leads, with chemical potentials µL,R. (c) Dissipative processes due to phonon emission (dashed arrows) responsible for rates in Eq. (22). Each processes is correlated with resonator photon creation (downward wiggly arrows) or annihilation (upward wiggly arrows); γ (ωq ±ωr ) ↓± correspond to DQD relaxation from |e to |g , whilst γ (ωr ) ϕ− corresponds to DQD dephasing, leaving the populations of |e and |g unchanged.

II. DRIVEN DISSIPATIVE RABI MODEL
The Rabi Hamiltonian describe a single two-level system or qubit coupled to a harmonic mode. It is one of the main workhorses of modern quantum physics and describes such diverse situations as natural atoms interacting with visible light [18], superconducting artificial atoms coupled to microwave resonators [35] or semiconductor double quantum-dots coupled to superconducting resonators [36]. Here we additionally consider the situation where the two-level system is coupled to an environment, inducing dissipation in the dynamics. Here we adopt the language of a semiconducting open double quantum-dot, but the technique presented here is applicable to many other situations that are described by the Rabi-model.
The Hamiltonian of a DQD coupled to a superconducting resonator, expressed in the DQD position basis where σ (p) z = |L L|−|R R| is the dipole operator of the DQD (also referred to as a qubit), σ (p) x = |R L| + |L R| induces transitions between its two charge states and a is an annihilation operator for microwave photons in the cavity mode at frequency ω r . Here the position states indicate the presence of a single charge either on the right or on the left dot. The DQD has an asymmetry energy between its position states of q and they are tunnel-coupled with strength ∆ q . The coupling between DQD and resonator is described by the dipolar interaction, with coupling strength g. Here, as in the rest of the paper we use the convention = 1.
The DQD couples to a bath of environmental modes, b j , with eigenfrequencies ω j via the system-bath coupling term For descriptive purposes we will assume the environment to be a bath of phonons as this is generally the dominant source of dissipation for charge-based quantum-dots [5].
In general, the DQD environment could include other electronic fluctuations, such as charged intrinsic two-level defects [37,38] or charge traps on surfaces [39].
In addition to the coupling of the DQD to a bath, the microwave resonator will also be coupled to a dissipative environment leading to loss and spurious excitation of resonator photons. Also, in the experiments motivating this work, the quantum-dot was subject to an external bias, leading to charge transport through the DQD as illustrated in Fig. 1(b). Both these processes will be included in the master equation and we describe them in sections III B and III C.
We write the interaction Hamiltonian H I as a sum of the intra-system interaction H I,S (i.e. between the qubit and resonator), and the interaction between the system and phonon environment H I,B , where we have decomposed σ (p) z = cos θ σ z + sin θ σ x into the eigenoperators, σ x = |e g| + |g e| and σ z = |e e| − |g g| of the bare system-bath Hamiltonian where the bare qubit level splitting is ω q = ( 2 q + ∆ 2 q ) 1/2 and the qubit mixing angle tan θ = ∆ q / q , with the convention 0 ≤ θ < π. The qubit eigenstates are given by In the following we perform perturbation theory in the interaction Hamiltonian H I . We are thus assuming that the qubit resonator coupling is weak compared to their detuning, g < |ω q − ω r |. Additionally we make the usual weak coupling approximation between the system and the phonon environment. The standard situation in experiments is ω > g β j , where ω represents all relevant system and bath frequencies.
The objective of this paper is to derive a Lindblad master equatioṅ for the evolution of the (slowly varying) system density matrix, ρ, where the index m labels orders of perturbation theory. In particular, L 2 includes dispersive and dissipative terms at second order in H I , L 4 includes dissipative terms at fourth order in H I , L res accounts for the resonator decay and L leads ρ describes a transport bias across the DQD, as in Ref. 1. The superoperators L m are composed of linear combination of coherent evolution processes, expressed as commutators with some effective Hamiltonian, and dissipative evolution, expressed in the form of Lindblad dissipators For the purposes of this paper, we will refer to any term of the form in Eq. (9) as a Lindblad dissipator, though we note that at higher order, we find such terms with negative coefficients.
Keldysh perturbation theory provides a systematic way to evaluate self-energy terms in the Dyson equation, in order to derive the dynamical master equation for the system density matrix. To determine the system evolution, we evaluate the Keldysh self-energy superoperator Σ, which is written as a perturbation series in the interaction Hamiltonian. Terms in this series are expressed in the graphical language of Keldysh diagrams, in order to keep track of all relevant contributions and avoid double counting.
To begin the analysis, the Keldysh master equation for the density matrix in the Schrödinger picture, ρ (S) (t), is given bẏ where the superoperator Σ is the self-adjoint self-energy, which acts on the state from the right. A derivation of Eq. (10) can be found in Appendix A 1. In the interaction picture defined by H 0 , this becomeṡ Hereafter, we drop the interaction picture label (I) . In what follows, we seek to express the self-energy superoperator on the RHS of Eq. (11) in the form of dispersive terms (i.e. commutators with the Hamiltonian) and dissipative terms (i.e. Lindblad superoperators), as in Eq. (8).
A. Keldysh self-energy As described in Appendix A 1, the self-energy superoperator, Σ(t 1 , t), is given by the sum of irreducible superoperator-valued Keldysh diagrams, each of which evolves the state from an early time t 1 to the current time, t. We evaluate Σ(t 1 , t) to a given order of perturbation theory by truncating the sum at that order.
Keldysh superoperators are represented diagrammatically as two parallel lines representing the time-evolution of the forward ( ·|) and reverse (|· ) components of the density matrix. Vertices placed on the lines correspond to the action of the interaction Hamiltonian, so at m th order of perturbation theory, there are m interaction vertices. These act at specific interaction times, t j , over which we integrate. Tracing over bath degrees of freedom corresponds to contracting vertices together, indicated by dashed lines.
To evaluate the time integral appearing in the Keldysh master equation, Eq. (11), we note that each Keldysh diagram takes the form of a convolution kernel with oscillatory time-dependence, so that the Laplace transform in the time-domain becomes a product of terms in Laplace space, i.e. ρ(t) → ρ s . In Laplace space, Eq. (11) becomes where ω j labels a set of system frequencies.
One approach to the solution of these equations was presented in Ref. 7 starting from the ansatz In brief, the unknown residues in Eq. (13) are found by comparing the poles on both sides of Eq. (12), and imposing a self-consistency condition thereupon. This selfconsistent solution ultimately truncates the summation in Eq. (13), which corresponds to a secular or rotatingwave approximation. Further details are given in Appendix A 3.
In the experiments that motivate this work, the steadystate response of the system determines the phenomenology, so we focus here on the slow dynamics of the system. Retaining only the pole at s = 0 corresponds to Two examples of irreducible fourth-order Keldysh diagrams. All 32 irreducible diagrams at this order can be generated from the two shown by swapping sub-sets of interaction vertices between the lower and upper lines.
the conventional secular approximation made in deriving the quantum optical master equation [7], so we use the quasi-static ansatz so thatρ captures the quasi-static evolution associated with the weak coupling to the bath. We now describe the results of evaluating L m at each order.

First order
There are no dissipative terms surviving the perturbative treatment at first order in perturbation theory, since for any thermal state of the environment one finds b j th = b † j th = 0 [18,25]. This is generally true for any odd moments of bath operators.
In our case, the perturbative interaction Hamiltonian H I also contains coherent terms between the qubit and resonator, i.e. terms that couple two parts of the system to each other. In this case the above reasoning no longer applies and some terms might survive even in first order of perturbation theory. However these terms carry a rapidly rotating time dependence. Making a rotatingwave approximation (consistent with our approximate ansatz, Eq. (14)) eliminates these terms also. It follows that L 1ρ = 0.

Second order
In second order perturbation theory there are four classes of diagrams to consider, depicted in Fig. 2(a). Evaluating these diagrams generates a dissipative contribution arising from system-phonon interactions, as well as a coherent contribution from the dispersive, coherent terms arising from qubit-resonator coupling. We provide worked examples evaluating these diagrams in Appendix A 4. At second order we write, where we specify L 2,diss and H 2 below.
2nd order dissipative terms If the vertices in Fig. 2(a) represent the interaction with the underlying phonon bath, H I,B , the diagrams yield the usual dissipative terms in the master equation for a two-level system coupled to an environment [18,40]. They evaluate to (c.f. Appendix A 4 a) (16) where the second-order rates are γ ↓,2 = sin 2 θ C(ω q )/2, and in which the spectral function of the phonon environment is is the bath coupling operator in the interaction picture, J(ω) is the spectral density, and b † j b j = n th (ω j ) = (e βωj − 1) −1 is the Bose distribution with β = 1/k B T . θ(ω) is the stepfunction. See Appendix A 2 for more details.
2nd order coherent terms If the vertices in Fig. 2(a) represent dot-resonator interactions from H I,S , most diagrams carry fast oscillatory terms, so that they vanish in a rotating wave approximation. The residual, stationary terms contribute to the imaginary part of the self-energy, leading to a renormalisation of the Hamiltonian. We find (c.f. Appendix A 4 b) where the dispersive shift is χ = g 2 sin 2 θ/(4ω q − 4ω r ), similar to the standard treatment of the Jaynes-Cummings model in the dispersive regime [35]. The difference between our result Eq. (18) and the usual Jaynes-Cummings dispersive shift lies in the inclusion of counter-rotating terms in the system-bath coupling Eq. (5), leading to the additional ratio of system frequencies in Eq. (18). Our result does not make the rotatingwave approximation in this coupling and thus includes both the usual dispersive as well as the Bloch-Siegert shift [41,42] in one analytic expression.

Third order
Third order diagrams can again contribute due to the system-system coupling terms in our perturbative interaction Hamiltonian. Their contributions can be at order ∼ g 3 or ∼ gβ 2 , i.e. they are at even order in H I,B and odd order in H I,S . These terms contribute Hamiltonian corrections, i.e. small energy shifts and changes of the bare system parameters, so we ignore them here. Instead, we assume these contributions are implicitly included in renormalised experimental parameters. Thus we take L 3ρ = 0. 1. Correlated photon-phonon decay processes: Diagrams in which two vertices come from the intrasystem interaction Hamiltonian H I,S and the remaining two are from the system-bath interaction H I,B .
2. Coherent dispersive photon-photon processes: Diagrams in which all four vertices originate from the coherent intra-system interaction H I,S .
3. Two phonon decay processes: Diagrams in which all four vertices come from the system-bath interaction term H I,B We thus write where the order of the terms corresponds to the numbered list above.
In this paper we are only concerned with calculating the self-energy diagrams leading to correlated photonphonon decay processes, L 4,corr . We therefore evaluate neither H 4 nor L 4,2phonon here, but leave this to future work. In passing, we note that H 4 will renormalise the bare Hamiltonian terms, and L 4,2phonon will renormalise the rates in Eq. (17), and potentially include dispersive shifts as well.
4th order correlated photon-phonon decay terms Evaluating correlated photon-phonon decay selfenergy diagrams at fourth order results in a total of 21 individual dissipators, each associated with rates Γ j .
An example diagram at fourth order is evaluated in Appendix A 5 a for pedagogical purposes. After evaluating all the relevant fourth-order Keldysh diagrams, we are left with a large number of superoperators that contribute to L 4,corrρ . Collecting these into Lindblad dissipators requires careful analysis, which we describe in Appendix A 6.
The full set of correlated photon-phonon decay terms that appear at fourth order of Keldysh-Lindblad perturbation theory are (20) Each of the rates Γ j are sums of terms depending on the bath spectral function, C(ω), evaluated at frequencies ω ∈ {0, ±ω q , ±ω r , ±(ω q − ω r ), ±(ω q + ω r )}, c.f. Table B1. Additionally there are contributions to the Γ j that are proportional to the derivative of the bath spectral function, C (ω), evaluated at frequency ω = ±ω q . Those originate from diagrams where the time evolution involves degenerate frequencies, leading to double poles in Laplace space, as has been noticed before [43]. Their physical interpretation is uncertain. One possibility is that they are first order terms in a Taylor expansion of the bath spectral function and as such indicative of an implicit renormalisation of system frequencies. We have verified that the derivative terms do not contribute significantly in the experimental parameter regime we consider below and leave their detailed understanding to future work.
Analysing all 21 dissipators above in detail is not the main concern of this paper. However, the first six dissipators in Eq. (20) are easy to interpret, and have dominant contributions, γ j , to the corresponding rates Γ j , that can be expressed simply. These are where ↑ / ↓ denote qubit flipping processes, ϕ indicates qubit dephasing, and +/− photon excitation and loss. The rates are The remaining three rates in Eq. (21) are related by thermal occupation factors, In section V, we analyse in detail the contributions from Eqs. (22) to gain and loss in the parameter regime of recent experiments in semiconductor double quantum-dots coupled to a microwave resonator. There we show that the six dissipators in Eq. (22) quantitatively explain recent experimental results [1], whereas including the full set of rates from Eq. (20) does not change the picture qualitatively apart from a change in numerical parameters.

B. Resonator decay
The resonator couples weakly to the external electromagnetic environment [44], and we include this as an additional decay process in the usual way [25] through the Lindblad dissipators with the resonator relaxation and excitation rates κ ∓,r and the bare resonator linewidth κ = κ −,r − κ +,r .

C. External quantum-dot bias
In the experiments that motivate this work, the open DQD coupled to external leads, inducing a chargedischarge transport cycle, illustrated in Fig. 1(b). To model this process, we extend the DQD basis to include the empty state |Ø , in which the DQD is uncharged.
As electrons tunnel between the leads and the DQD, it passes transiently through the empty state. This process is described by the incoherent Lindblad superoperator [5,17,18,21,22] For simplicity, we will assume Γ R = Γ L = Γ in the following. Depending on the sign of q and the strength of Γ, the DQD population may become inverted in steady state.

IV. DRIVEN RESONATOR
The additional dissipators will tend to drive the system to its ground state. To see non-trivial behaviour arising from the master equation, the system needs to be driven out of equilibrium. Depending on the specific system, there are a variety of ways to achieve this. Here we consider resonator driving, in which an external microwave field is imposed on the resonator. To this end we introduce an additional resonator driving term in the system with the amplitude of the drive d and frequency ω d . In writing Eq. (25) we have assumed a rotating wave approximation, relevant to near-resonant driving ω d ≈ ω r . Under driving, we anticipate that the resonator will tend to a steady-state which is close to a coherent state, |α . We therefore apply a displacement transformation on the resonator modes [45,46], a →ã + α. This transformation forms the basis of a semi-classical, meanfield approximation for dissipative steady-state equations for qubit and resonator, from which we establish a selfconsistency condition that determines α.
The canonical unitary for the displacement transformation is so that D(α) † aD(α) =ã + α. The Hamiltonian in the displaced frame is transformed as where α may be time-dependent. When applying the displacement operation, we have to take care to include dissipative processes containing resonator operators in the transformation. For example, when applying the transformation to the dissipator describing direct photon decay ∼ D[a]ρ we find where the second term describes a coupling between the residual modeã and the coherent field amplitude α. This term, and similar ones from other dissipators, contributes to a self-consistent equation for the coherent field amplitude α. Correlated dissipators transform similarly. For example in the displaced frame, D[σ − a † ]ρ transforms as A. Separable and semi-classical approximations In principle, we could solve the steady-state (or dynamical behaviour) of the master equation, Eq. (8), accounting fully for correlations between the resonator and qubit. However, to make analytical progress we now make a simplifying mean-field approximation that enables us to calculate the qubit and resonator steadystates.
That is, we assume the resonator to be in the semiclassical regime, so that it is well approximated by a coherent state, |α . We will self-consistently determine the coherent amplitude α by finding a displaced resonator frame which is effectively undriven.
In this approximation, the resonator and qubit are separable,ρ =ρ r ⊗ρ q .
Under this approximation, Eq. (8) decomposes into two coupled mean-field equations, one for the resonator state, ρ r = |α α|, and one for the qubit state,ρ q , in which each sub-system experiences the mean field of the other. We findρ Evaluating these traces simplifies correlated dissipators. For example, tracing over the qubit in Eq. (29) contributes a term on the RHS of Eq. (31): where P i = Tr q {ρ|i i|}, so that σ + σ − = P e . The first term describes incoherent photon generation originating from the correlated decay process D[σ − a † ]ρ, while the second and third terms represent effective resonator driving, which depends on the qubit population and the resonator field expectation value α.
Finally, tracing Eq. (29) over the resonator degrees of freedom results in a contribution to the qubit master equation, Eq. (32)

B. Coupled steady-state equations
In the application to gain and loss in a driven DQDresonator, we are interested in steady-state properties of the system, so we solve the above equations subject tȱ ρ = 0, andα = 0. The mean-field approximation for the resonator field then yields coupled steady-state equations for the two sub-systems. In a frame rotating at the resonator drive frequency ω d , we find the steady-state of the resonator master equation, Eq. (31), becomes with the renormalized photon loss and generation rates and where the contributions from the fourth-order Keldysh perturbation theory are given by Here σ z = P g − P e is the expectation value of the qubit dipole operator and κ ±,r are the original resonator decay and excitation rates. P Ø is the population of the empty state, when no charge is on either island. For the parameter regime we consider below, P Ø 1. The full expressions for the rates κ (ωq) ϕ and κ (±ωq) are somewhat cumbersome and are given in appendix B 2 where we also give their relationship to the Γ j in Eq. (20).
The effective Hamiltonian for the resonator in the displaced frame is The resonator linewidth is renormalized by the interaction with the qubit to The term in Eq. (38) proportional toã † (and its hermitian conjugate) describes effective driving of the transformed resonator modeã, which would lead to a non-zero expectation value of the resonator operators ã . We now self-consistently require that the coefficient ofã † should vanish. Since the effective driving terms depend on α, this choice leads to the condition In this case, the effective driving terms in Eq. (38) are zero and Eq. (35) describes the residual quantum dynamics of an un-driven resonator subject only to relaxation at rate κ − and excitation at rate κ + In the displaced frame with the appropriately chosen value of α, the displaced resonator modeã will be in a thermal state, with its effective temperature T r depending on the ratio of its relaxation and excitation rates like In this case, we find ã = ã † = 0. In the following we additionally assume a small effective temperature of the resonator T r such that ã †ã = κ+ κ−−κ+ 1. The steady-state of the qubit master equation, Eq. (31), then becomes where the effective qubit Hamiltonian is and the correlated qubit dissipation rates are given by with the additional rates from correlated fourth-order processes The expressions for γ , and γ ϕ are somewhat lengthy and are given in full in Appendix B 2, where we also discuss their relationship to the Γ j 's in Eq. (20).
If in the previous calculation the effective resonator temperature T r is large, i.e., when κ + κ − , the expectation value of ã †ã acquires a non-zero value |α| 2 . If necessary, we can simply take this into account by noting that tracing over the resonator degrees of freedom now involves additional factors of |α| 2 . As an example Eq. (34) becomes whereα is determined from the effective thermal state at temperature T r .

V. APPLICATION TO GAIN AND LOSS IN MICROWAVE-DRIVEN DOUBLE QUANTUM-DOTS
We now apply this theory to a specific experimental situation, where an open semiconductor double quantumdot is coupled to a microwave resonator and subject to a transport bias. We are interested in calculating the steady-state resonator field depending on qubit parameters in an effort to replicate the experimental findings of Ref. 1. In particular, we calculate the gain of the DQDresonator system, by comparing the internal circulating microwave energy within the driven DQD-resonator with the same quantity for a driven resonator alone (i.e. without the DQD).
From Eq. (40) we find the resonator steady-state field amplitude is described by where δω r = δω r + 2χ σ z is the DQD-renormalised detuning. We calculate α by simultaneously solving for the qubit steady-state, Eq. (41) and for the resonator field, Eq. (45). The power gain is then given by is the steady-state resonator field that would be produced in the absence of the DQD coupling, (i.e. setting g = 0). We set δω r = 0, corresponding to resonant cavity driving, which is the experimental situation reported in Ref. 1.
In order to make a quantitative comparison between theory and experiment, we assume the bath spectral function is C(ω) = J(ω) (n th + θ(ω)) with J(ω) = J 1D (ω) + J 3D (ω), where the spectral densities for the first phonon mode in the quantum wire and the bulk substrate phonons are given by [3,5] with the speed of sound in the quantum wire c n = 4000 m/s and in the substrate c s = 5000 m/s [47], an inter-dot spacing d = 120 nm and nanowire radius a = 25 nm. We have taken k B T /ω r = 7.8, from [3]. We discuss details of fitting parameters F, P, and a gaussian smoothing parameter w below.

Theoretical and Experimental Gain-Loss Profiles
In Fig. 3(a) we plot the gain G in the microwave resonator due to the DQD coupling for three different theoretical treatments, along with experimental gain data extracted from Ref. 1, shown as points.  The three theoretical curves shown in Fig. 3(a) compare the dependence of the gain on various terms in the fourth-order correlated dissipators. The dashed blue curve includes terms in the first two lines of Eq. (21) but not the third line (i.e. the fourth-order theory restricted to γ (∓ωr) ϕ± = 0), which is equivalent to the polaronic theory in Ref. 3. The solid black curve further includes the six terms in Eq. (21), corresponding to DQD-mediated photon-to-phonon interconversion. The dotted red curve finally includes all 21 fourth-order rates that appear in the derivation of the full master equation, Eq. (20).
In Fig. 3(a), there is a clear qualitative difference between the polaronic theory (blue, dashed) and the experimental data in the regime q /ω r 0: the theory is unable to explain the depth of loss (sub-unity gain). In contrast, the additional dephasing-mediated processes in the last line of Eq. (21) give rise to enhanced losses beyond the polaronic terms, and are sufficient to quantitatively account for the entire range of gain and loss observed in the experimental data (black). We have also shown the results of the full fourth-order theory (red dotted) from Eq. (20), which also consistently captures the peak of the gain for q > 0, and the depth of the trough for q < 0. In the parameter regime described here, there are only small differences between the latter two curves, including a rescaling of F and P, indicating that the six terms in Eq. (21) quantitatively account for the extremes of gain and loss. Fig. 3(b) plots the rates in Eq. (22), and shows clearly that the dephasing-assisted loss rate (solid), γ (ωr) ϕ− , is significant compared to the other correlated decay processes, γ (ωq∓ωr) ↓± (dashed and dotted). This is the main reason for the qualitative difference between the theory curves in Fig. 3(a). Since the full fourth-order theory accounts for the experimental data, whilst the polaron-only theory does not, we conclude that dephasing-assisted loss is a substantial contribution to the system dynamics. Fig. 3(c) shows the steady-state DQD population imbalance, σ z ss (black), compared to the thermal equilibrium value, σ z th = Tr{σ z e −βH S }/Z (blue). In conventional gain/loss models, positive values of the difference σ z th − σ z ss (i.e. population inversion) drive gain, while negative differences drive loss. This is manifest in κ ± , where enhancement of P e leads to an increase of κ + , a corresponding decrease in κ , and thus an overall increase of α. In contrast, the dephasing-assisted loss contribution to Eq. (37) is independent of the state of the DQD in the limit that P Ø 1, which is the case here.

Note on Parameter Fitting
As in Ref. 3, we treat the dimensionless spectral strengths F and P as free parameters. We follow the same fitting approach, in which we choose parameter values so that the theory curves satisfactorily replicate the strong gain peak evident for q > 0. The specific values are shown as insets in Fig. 3(a).
Likewise, we also convolve the bare theory gain curves with a normalised Gaussian smoothing kernel ∝ e − 2 q /2w 2 (corresponding to a full-width-at-half-maximum of √ 8 ln 2 w) [3], to account for low-frequency voltage noise in the gates defining the inter-dot bias [36].
In effect, all theory curves have 3 fitting parameters, F, P, and w. Very roughly, F controls the gain peak height, P controls the tail behaviour of the gain at large q , and w controls the gain peak width. Thus, once we find a parameter set for a given theory curve that adequately accounts for the gain profile for q > 0, the dependence for q < 0 is determined. The overall scale for F and P is ultimately constrained by the value of the resonator linewidth κ: it is the relative contributions from the correlated decay rates compared to κ that set overall gain and loss through the effective resonator line-width κ , c.f. Eq. (39).
To make contact with the microscopic electron-phonon coupling, Appendix C gives estimates of the values of F and P calculated using material coupling constants, den-sities and speeds of sound, assuming simplified geometries for phonon modes of the 1D InAs wire and bulk SiN substrate [5]. We find F piezo = 0.85 and P piezo = 0.16, which are in reasonable order-of-magnitude agreement with values shown inset in Fig. 3(a).
The high phonon temperature assumed in Fig. 3(a) (T = 3K) may arise from ohmic heating from current flowing through the DQD, as estimated in Ref. 1. Increasing the phonon bath temperature in the theory leads to an increase in the loss rate, with commensurate improvement in the agreement with the experimental data for q /ω r < −3, as shown in Appendix C 1.

VI.
DISCUSSION AND CONCLUSIONS

Effective Rates
The large number of rates that appear in Eq. (20) is too great to analyse individually in detail, however the mean-field approximation embodied in Eq. (30) helps to distil significant combinations of the rates in Eq. (20) into effective decay rates on the resonator, Eq. (37), and the qubit, Eq. (43), separately. Figure 4 (and Figs. B1 and B2 in the appendix) show plots of the effective rates appearing in the qubit and resonator steady-state equations, Eq. (35) and Eq. (41), in the region ω q ∼ ω r . Many of the rates show divergences at resonance ω q = ω r or become negative in the vicinity of very small detuning. This behaviour is not unexpected, since our perturbative series relies on a small qubit-resonator coupling compared to the qubit and resonator detuning, g < |ω q − ω r |. Proximity to resonance between qubit and resonator will therefore make it necessary to take higher orders in perturbation series into account to accurately describe the system dynamics.
The effective resonator decay rates κ ±,4 are shown in Fig. 4 in various limits of P e and for an ohmic bath spectral function, C(ω) = ω(n th + θ(ω)), both for zero and non-zero bath temperature T . For a given choice of qubit parameters, and temperature, the 'true' value of the effective rates κ ±,4 is a convex combination of the limits plotted. For ω q > ω r , these effective resonator rates are all positive, so that the Lindblad operators generate well-defined completely positive (CP) maps on the resonator.
As discussed before, the displaced resonator modeã is un-driven such that its steady-state will be a thermal state at an effective temperature T r defined by the ratio of its relaxation and excitation rates κ + /κ − = e −ωr/(k B Tr) . Due to the effective coupling of resonator and qubit steady-states through the fourth-order dissipative rates, this effective temperature will be influenced by both the original resonator temperature as well as the temperature of the qubit environment. The effective resonator temperature depends specifically on the qubit steady-state population P e as seen from Eq. (37) and Fig. 4. In the case of population inversion in the qubit, κ−,4/g 2 @ Pe = 0 κ+,4/g 2 @ Pe = 0 κ−,4/g 2 @ Pe = 1 κ+,4/g 2 @ Pe = 1  the effective resonator excitation rate κ +,4 can dominate over κ −,4 .
In appendix B 3 we show plots of the correlated qubit relaxation, excitation and dephasing rates from Eq. (43). Their behaviour is similar to that of the effective resonator rates shown in Fig. 4.

Negative Lindblad superoperators
For ω q ω r , there are regimes in which the κ ±,4 can become negative. The situation is similar for the effective qubit rates γ ↓/↑,4 and γ ϕ,4 , as shown in appendix B 3. More generally, Eq. (20) has a number of Lindblad superoperators with explicitly negative coefficients. This can be seen e.g. in the prefactors of D[σ z a † a]ρ and D[σ z + σ z a † a]ρ, which have a relative minus sign relating them, so that one coefficient will always be negative. This mathematical phenomenon arises naturally at fourth-order perturbation theory, and has previously been connected to the description of explicitly non-Markovian dynamics [48][49][50][51] and its unravelling into a quantum jump description [52]. In contrast, the second order rates in Eq. (17) are always explicitly positive.
The appearance of negative coefficients to Lindblad superoperators is controversial, since there is no longer a guarantee that the resulting map is completely positive.
There are a number of possible resolutions to this mathematical problem. Firstly, the dominant fourthorder terms (i.e. those in Eq. (21)) are indeed explicitly positive, so that they do correspond to well-behaved Lindblad superoperators. That is, the most significant physics is represented by a CP map.
Secondly, the main objective of this paper is to find steady states of the propagator, so that we do not necessarily require complete positivity. Rather, we simply require that the steady state is a positive operator, as is indeed the case in the situation described here, illustrated by the red dotted curve in Fig. 3. However, this argument is unsatisfactory in general, since the dissipators can be interpreted in a dynamical equation, just as L 2ρ in Eq. (15) appears in the dissipative dynamics of the quantum optical master equation.
Thirdly, the explicitly negative signs in Eq. (20) are associated with somewhat anomalous dissipators such as D[σ z + σ z a † a]ρ, which on their own contain products of four resonator creation and annihilation operators (e.g. the expansion of this Lindblad dissipator includes the operator product σ z a † aρσ z a † a). These dissipators all arise from the diagonalisation procedure described in Appendix A 6 a. In fact, such dissipators always appear in combination with another dissipator (in this case D[σ z a † a]ρ), such that quartic products of resonator operators cancel. At higher orders of perturbation theory, we expect to see dissipators that are truly quartic in the resonator operators appearing, and these would renormalise the negative coefficients seen at the order we evaluate.
This suggests that the apparent negative signs appearing in Eq. (20) are a consequence of the truncation of an infinite power series; including higher-order terms should renormalise the dissipative rate at each order. As such, we speculate that resumming a large class of higher order Keldysh diagrams will lead to well defined Lindblad superoperators with positive prefactors.
This can be understood as an analogue to the seriesexpansion of a unitary operator, U = e −iH generated by a Hamiltonian H. While the zeroth and infinite order expansion are both explicitly unitary (corresponding to the physical requirement that probability is conserved), this is not guaranteed to hold for any finite order (e.g. the approximation U ≈ 1 − iH is not generically a unitary operator), so that conservation of probability is strictly violated. Nevertheless, linear-response theory explicitly makes these kinds of approximations, to good effect. By analogy, it is not surprising that the fourthorder Keldysh-Lindblad master equation we have derived includes terms that are non-CP, which are nevertheless useful for describing the system.

Keldysh-Lindblad synthesis
Equation (20) is the principal theory result of this paper. It represents the synthesis of the Keldysh and Lindblad formalisms, which both are widely used in mesoscopic physics, but with relatively little overlap in the community of practitioners who deploy them; the Lindblad formalism has been carried over from the quantum optics literature, whilst the Keldysh formalism has its roots in the non-equilibrium condensed matter literature.
Equation (20) contains a large class of Lindblad superoperators that go beyond the standard dissipators in the quantum-optical master equation. The first six terms in Eq. (20) are dominant in a typical experimental situation, and these have a natural interpretation as qubit mediated energy exchange between the resonator and the bath. The first four terms have previously been identified using a polaron transformation [3]. The terms γ ϕ+ D[σ z a † ]ρ + γ ϕ− D[σ z a]ρ, which correspond to a dephasing assisted energy exchange, represent the most significant new contribution presented here.

Conclusions
We conclude that the new dephasing-mediated gain and loss Lindblad superoperators in Eq. (21) account for substantial additional loss observed in recent experiments. These terms were derived using the Keldysh diagrammatic techniques, and arise at the same order of perturbation theory as other terms previously derived using a polaron transformation. Including all additional dissipators derived at the same order in perturbation theory does not qualitatively change the results apart from a rescaling of numerical parameters.
Synthesising Lindblad and Keldysh techniques to derive higher-order dissipative terms is thus a powerful approach to a consistent, quantitative understanding of quantum phenomena in mesoscopic systems. Lifting the simplifying mean-field approximation to study the effects of correlations between the DQD and resonator will be the subject of future work.
This paper opens a variety of avenues to explore. We have not evaluated the fourth-order dispersive terms, nor the fourth order pure bath (multi-phonon) terms. In some regimes, two-phonon emission may be significant, particularly if phononic engineering is used to suppress the spectral density at dominant single-phonon frequencies (e.g. by designing phononic band-gaps). In such a situation, two (or more) phonon emission rates can be calculated within the Keldysh formalism, and collected into Lindblad superoperators.
In future work we intend to analyse the dynamics and steady-state of the fully correlated equation, Eq. (20), without the simplifying mean-field approximations, Eq. (30). Also, the master equation Eq. (20) naturally lends itself to an input-output treatment allowing the calculation of qubit and resonator output correlation functions.

Derivation of Keldysh master equation
We describe a system interacting with an environment using Keldysh diagrams. The basic Hamiltonian has the form where H 0 = H S + H B describes the unperturbed system H S and the environment H B , and H I contains all the interactions terms and is considered small. In the following we will employ an interaction picture with respect to H I , where operators now evolve in time with respect to H 0 while states evolve according to H I . In the Schrödinger picture, the time evolution of the systems density matrix elements is where we have traced out the bath degrees of freedom, defined the projector unto system states P kk = |k k | with where U 0 ≡ U 0 (t, t 0 ) = e −iH0(t−t0) , and we have defined the forward and backwards time-evolution operators in the interaction picture by where T indicates time-ordering (later times left) and T anti-time-ordering (later times right). The interaction term H I (t) is itself in the interaction picture, i.e., evolves in time according to the unperturbed Hamiltonian H 0 as We have also assumed in Eq. (A2) that the system and bath factorise at the initial time, i.e.
We now expand the time-evolution operators U † I and U I on both sides of P kk (t) in Eq. (A2) in a perturbative expansion in powers of the coupling coefficients g and β j that appear in H I . At a given order, m, of perturbation theory in this expansion, the propagator Π acts on the density matrix by application of m instances of H I , which appear on the left and right of ρ in every possible configuration.
This expansion has a graphical representation in the form of Keldysh diagrams, shown in Fig. A1. The action of H I in Eq. (A2) acting on the left (which come from U † I ) are drawn as vertices on the lower branch of the diagram; those acting on the right (which come from U I ) are drawn as vertices on the upper branch. These diagrams are operator-valued, schematically acting on ρ(t 0 ) from the right. Diagrams up to fourth-order are shown in Fig. A1a. Here we have already made use of Wick's theorem to express multi-operator bath correlation functions as sums of two-time correlation functions, indicated by red dashed lines in the diagrams of Fig. A1. Fig. A1(b) and (c) show a diagrammatic refactoring of the expansion, to derive the Dyson equation, which relates the propagator, Π, to the self-energy Σ, and the free evolution Π 0 through Π(t 0 , t) = Π 0 (t 0 , t) + dt 2 dt 1 Π(t 0 , t 2 )Σ(t 2 , t 1 )Π 0 (t 1 , t).
(A7) Diagrammatically, the self-energy Σ is given by the sum over all irreducible diagrams, shown in Fig. A1(c).
When we substitute the the Dyson equation (shown in Fig. A1d) into the propagator in Eq. (A2), we find the time-dependent evolution of the density matrix is given by Figure A1. (Colour online) (a) Illustration of the density matrix propagator Π in diagrammatic perturbation theory for evolution from the initial time t0 to the final time t. Other times depicted are dummy integration variables, satisfying t0 < ... < t2 < t1 < t. Lines with arrows correspond to periods of free evolution. Dots denote the interaction Hamiltonian, HI (tj) acting at time tj. denotes swapping of upper and lower vertices. Background colours are used to help track terms in the following panels. Only terms at even orders in perturbation theory are shown, and dashed lines denote contractions of bath operators into two-time correlation functions. In our convention, Π acts on the density matrix from the right, so that ρ(t) = ρ(t0)Πt 0 →t, as in Eq. (A8). (b) Diagrammatic reordering of the Π. Here we factor each diagram into products of irreducible diagrams, which are those that cannot be divided by vertical lines without cutting any contractions. (c) Diagrammatic expression of the Dyson equation. The first term in the product (blue diagrams) is the propagator, Πt 0 →t 2 , and the second term (green and orange diagrams) is the self energy, Σt 2 →t 1 , composed of irreducible self-energy diagrams, followed by a period of free evolution. where in the second line the free propagator has been explicitly written as Taking the time-derivate of Eq. (A8) finally leads to the Keldysh master equation where E k is the energy of system state |k . From this equation we see that the self-energy superoperators Σk k →k k generate the residual dynamics arising from the interaction hamiltonian, and are responsible for both dispersive and dissipative effects at higher order. Note that the self-energy in Eq. (A10) is still formally exact, since it contains all orders of perturbation theory.
Noting that ρ(t) = k k ρ k k (t) |k k|, we write the master equation for the full density matriẋ where the self-energy superoperator Σ(t 1 , t) acts on the density matrix from the right. This convention ensures that the Keldysh diagrams we draw later are correctly understood as acting with earlier times on the left. In the interaction picture, we havė In what follows, we will work in the interaction picture, and drop the notation (I) .
The detailed relationship between the last term in each of Eq. (A10) and Eq. (A11) is established by defining Σk k →k k (t 1 , t) ≡ j k L j |k k |R j k for some set of time-dependent operatorsL j ,R j . Then We note that in order to preserve hermiticity of ρ, the self-energy operator is self-adjoint.
In what follows, we evaluate the self-energy operator up to fourth order, and collect terms in order to express the master equation in Lindblad form.

Note on environmental spectral functions
As in most treatments of open quantum systems, the bath spectral function plays a critical role in capturing the open-systems dynamics. In anticipation of later derivations, here we define the bath spectral function through the Fourier transformation of the two-time correlation function of the bath coupling operators as such that in all integral expression we can replace Here,X(t) is the bath coupling operator in the interaction picture, which in our case, for bosonic phonon modes, takes the form Defining the coupling operator spectral density J(ω) in the continuum limit of the bath modes with spectral density ν(ω) as and assuming the bath modes to be in thermal equilibrium with a bose distribution, b † j b j = n th (ω j ) = 1 e βω j −1 , we can explicitly evaluate the spectral function as which follows the detailed balance relation C(ω)/C(−ω) = e βω with the thermal factor β = 1/k B T .

Laplace space integration
To evaluate the self-energy, we transform the master equation into Laplace space. In Laplace space, the timeintegration is particularly simple, as all necessary timeintegrals take the form of convolutions. Additionally, we will focus our evaluation on finding the steady-state density matrix of the system. We make use of the technique developed in Ref. 7, which enables one in general to find dynamical steady-states, i.e. steady-states with a residual, periodic time-dependence. Here, we will only be interested in the time-independent steady-state. In the following we provide a short overview of the technique, more details can be found in Ref. 7.
We first transform the Keldysh master equation into Laplace space. Since the time-integrals can be written as multi-product convolutions, their solution in Laplace space is particularly simple. Defining where f s is the Laplace transformed function f (t), the Keldysh master-equation, Eq. (A12), in Laplace space becomes where ρ s is the Laplace transformed density matrix and the index j enumerates the different diagrams in the selfenergy. The operatorsô j depend on the specific diagram and the index k is dependent on the order N j of perturbation theory of the j-th term. The coefficients ξ j contain bath spectral functions C(ω) each of which comes with an integral over the Fourier frequency ω, c.f. App. A 2. Finally, the frequencies ω k are linear combinations of system frequencies from the time-evolution of the system operators, and Fourier frequencies from the definition of the bath spectral function. Examples of this expression for specific diagrams are evaluated in Eq. (A28) and Eq. (A35). To solve Eq. (A12) for the dynamical steady state of the system, the following ansatz has been shown to be useful [7] whereρ 0 is the time-independent part of the steady-state and theρ j components of the steady-state show residual oscillations at frequencies ω j . In Laplace space, this becomes where the limit lim s→0+ corresponds to calculating the residuals of the right hand-terms. Here we have assumed that the zero pole on the rhs of Eq. (A20) originates with the density matrix ρ s , and not from a term in the selfenergy Σ s , as is the case in all diagrams considered below. More details can be found in Ref. 7. When evaluating the residue of the zero frequency pole of the master equation, the following limit will become Figure A2. (Color online) One specific self-energy diagram in 2nd order, with each of the vertex operators specified. Here both vertices correspond to coupling to bath operatorsX, and their contraction is indicated by the blue solid line. This diagram contributes to the the standard qubit decay process ∼ D[σ−]ρ as detailed in the text.
where the δ-function allows us to restrict the frequency integral over the bath spectral functions C(ω) and the imaginary parts may contribute as frequency shifts to the effective Hamiltonian.

a. Second order dissipative Keldysh diagrams
To give a practical example of the general steps, we will show three specific examples. First we consider a simple dissipative diagram in second order shown in Fig. A2. Here we use the convention that when collecting the terms, we trace the lines in the direction of the arrows, starting from the bottom right. Also, remember that the density matrix at the earliest dummy integration time (here t 1 ) multiplies this diagram from the left. Evaluating this diagram leads to the the master equation term where for the sake of readability we neglected the prefactors ∼ sin θ.
We now assume separability of the bath and system, ρ(t ) = ρ q (t ) ⊗ ρ B over the interval t 0 < t < t 1 (n.b. the first vertex at the dummy time variable t 1 ), and that the bath is in its thermal steady state throughout this time interval. This gives Figure A3. (Color online) Another example of a self-energy diagram in 2nd order, with each of the vertex operators specified. Here both vertices correspond to coherent coupling between the qubit and the resonator, as indicated by the wavy green line.
which, after inserting the definition of the bath spectral function, Eq. (A15) becomes Here we have expressed the result in the form of temporal convolution, with f (t) = ρ q (t) and g(t) = exp (i(ω − ω q )t). Transforming into Laplace space allows us to evaluate the time integral and we get 2,diss ρ(t 1 ) with the Laplace transformed density matrix ρ s . Inserting the Ansatz, Eq. (A22), and taking the steady-state limit using Eq. (A24) leaves us with which is one of the terms appearing in the dissipator describing qubit relaxation γ ↓,0 D[σ − ]ρ, cf. Eq. (17). Here the imaginary term in limit s → 0 + will contribute only as a small renormalization of Hamiltonian parameters, which we neglect.

b. Second order dispersive Keldysh diagrams
We now show another example in second order, illustrated by the diagram in Fig. A3. Here, both vertices represent operators from the interaction between the qubit and the resonator. The contribution to the Keldysh master equation takes the form where prefactors are ignored. The integral can again be evaluated in Laplace space to arrive at Since there is no frequency integral now, the only relevant term in the limit s → 0 + is the imaginary contribution and we get the final result which forms part of the effective dispersive Hamiltonian H 2 , Eq. (18).

4th order diagrams
At 4th order in the interaction, a total of 32 irreducible diagrams contribute to the Keldysh self-energy. Half of these diagrams are schematically depicted in Figs. A4, from which the other 16 can be obtained by swapping all vertices between the upper and the lower line.

a. Fourth order dissipative Keldysh diagrams
We consider one particular dissipative diagram in fourth order, Fig. A5. This diagram acts on the state at t 3 , so that the corresponding term in the Keldysh master equation is where we again left off all prefactors for the vertex operators and already assumed the Born initial condition of separable system and bath density operators for t 0 < t < t 3 . Inserting the definition of the bath spectral functions gives This is a multi-product convolution over the integration variables. In Laplace space the time integrals can then be evaluated trivially and we get . (A35) Figure A5. (Color online) Particular realisation of a selfenergy diagram in 4th order, with each of the vertex operators specified. Here one pair of operators represents coupling to the dissipative environment and is connected by a blue solid line. The second pair of vertex operators describes coupling between the qubit and the resonator, and is connected by a wavy green line. This diagram contributes to the correlated decay process ∼ D[σ−a † ]ρ.
We now make use of a partial fraction decomposition of the integral kernel, with the added subtlety that one of the poles in Eq. (A35) is degenerate. We note that where we used shorthand ω 1 = ω − ω q + ω r . The second term deserves special attention since in the steady-state limit, s → 0 + , it leads to i.e., a degenerate pole leads to a derivative of the bath spectral function, as has been observed before [43]. Note that above we neglected the write the imaginary components appearing in the limiting procedure, since we will neglect those small Hamiltonian terms in the end. For the complete diagram Fig. A5, we thus get which contributes to the process of correlated qubit decay and photon generation. Specifically the first term is a relevant contribution to the primary rate expression γ (ωq−ωr) ↓+ , Eq. (22). We note in passing that the terms proportional to C (0) and C(0) as well as any imaginary parts in Eq. (A38) ultimately cancel with terms arising from other, similar diagrams.

Dissipators from master equation
This section describes how to obtain the Lindblad dissipators from the Keldysh steady-state equation. After evaluating the self-energy to a given order can write the Keldysh master equation in the forṁ whereσ is a row-vector of operators and the objects A are coefficients matrices in the space of these operators. For example, the scalar factor appearing in Eq. (A38) would contribute to A m . Hermiticity of the master equation requires that the matrices A be symmetric. Trace-conservation of ρ also constrains the A's, so that A s = A e . These conditions imply that the RHS Eq. (A39) can always be expressed as a sum of hermitian commutators, and lindblad dissipators D[ ], with real coefficients, γ j . They do not automatically constrain the γ j to be positive, which has implications for complete positivity of the evolution. We comment on this in the main text.
The vectorσ contains all operators that can be made up of combinations of terms in the interaction Hamiltonian H I at the appropriate order. Specifically in the case of the Rabi model and with the interaction Hamiltonian given by Eq. (5) up to fourth order in perturbation theory, we havê σ = 1, σ z , a † a, σ z a † a, a † , σ z a † , a, σ z a, σ − , σ + , σ − a † a, σ + a † a, σ − a † , σ + a, σ − a, σ + a † . (A40) As is obvious from our choice of parametrization of the matrices A, any diagonal term that appears in both A m as well as A s,e , directly corresponds to a dissipator in the master equation. The dissipative operator is then given by the eigenvector belonging to that particular diagonal element. As a simple example, collecting all terms in the RHS of Eq. (A39) that are proportional to the phonon spectral function at the difference frequency between the qubit and the resonator, C(ω q −ω r ), we find the matrices A only have a single non-zero entry: in the basis defined by Eq. (A40). The non-zero entry on the diagonal corresponds to the operator σ − a † in Eq. (A40), so we immediately identify this part as the Lindblad dissipator Similarly simple diagonal forms can be found for all terms probing the phonon spectral function at frequencies ±(ω q − ω r ), ±(ω q + ω r ) and ±ω r , leading to the six correlated dissipators in Eq. (21).

a. Non-diagonal forms for A
In general the matrices A in the RHS of Eq. (A39) are not purely diagonal, and some manipulation is required to bring them into a recognisable form in terms of sums of dissipators, D[ ]. Specifically in our case it is the cyclic property of the pauli matrices which makes this a nontrivial challenge.
To diagonalise the matrices A, we make use of a particular property of composite dissipators, namely that for any two operatorsô 1 andô 2 , we can write the dissipator of the sum of these operators as The last six terms in this expression correspond to symmetric off-diagonal terms in the matrices A in Eq. (A39), whereas the first two terms are purely diagonal. Thus, combinations of off-diagonal elements in the A's can be directly related to combinations of dissipators using Eq. (A43): Since the matrices A are symmetric, we can thus ascribe dissipators to all off-diagonal elements. All remaining terms are then already diagonal and correspond to dissipators of single operators of our basis vectorσ. In practise, we then perform this diagonalization procedure on the matrix A m according to the replacement rulê and verify that the resulting dissipator-like terms are sufficient to construct the matrices A e and A s . If additional terms are required to construct these matrices, those terms will appear as contributions to the Hamiltonian in the steady-state equation, ∼ −i [H,ρ], as opposed to dissipators.
b. Dissipators proportional to C(0) As an illustration we now consider an example where the matrices A are not diagonal to start with. In particular, collecting all terms proportional to the phonon spectral function at zero frequency, C(0), we find a contribution to the RHS of Eq. (A39) with the A m matrix taking the form with the 4 × 4 sub-matrices and where the matrix coefficients are listed in table A1. As which contributes to the dissipative superoperator at fourth order, L 4,corrρ .
As an aside, we note that some of the individual dissipators in Eq. (A49) contain four photon creation and annihilation operators, such as D[σ z a † a]ρ, the total expression including all terms only contains a maximum of two photon operators in each term, as is consistent with our choice of self-energy diagrams.
proportional to C(0). There are other terms proportional to C(−ω q ), C(ω q ), C (−ω q ) and C (ω q ). Generically, each of these terms take the form being block diagonal in the operator basis Eq. (A40), and we find A with the coefficients specified in Table A1. Performing the above specified diagonalisation procedure leads to the expression σA (−ωq) e. Dissipators proportional to C(ωq) The matrix representing terms proportional to the spectral function at the positive qubit frequency, C(ω q ), is represented by the blocks with diagonalisation resulting in a similar expression to Eq (A53), f. Dissipators proportional to C (ωq) The final remaining terms are proportional to derivatives of the spectral function. Terms ∝ C (ω q ) can be written as where the other two blocks are zero. Diagonalisation giveŝ The terms ∝ C (−ω q ) lead to with the other two sub-matrices again zero. Similar to Eq. (A57), diagonalisation leads to the Lindblad dissipatorŝ Appendix B: Fourth order master equation of the Rabi Model
Here we have additionally introduced the symmetrized and anti-symmetrized bath spectral functions C±(ω) = C(ω) ± C(−ω) and the unsymmetrized spectral function derivative

Complete coupled steady-state equations
Here we give the full expression for all rates and factors appearing in the coupled steady-state equations for resonator and qubit.
Appendix C: Phonon spectral functions from microscopic considerations Following the discussion in Ref. 5 we calculate expressions for the spectral function weights F and P in Eq. (46) assuming piezoelectric interactions between the phonons and the DQD, based on a simplified geometrical model of the semiconductor heterostructure. We find where µ 3 is the three-dimensional mass density of the material, and µ 1 = πa 2 µ 3 is the one-dimensional equivalent for a wire of radius a. Here we assumed that the interactions of substrate phonons with the electronic state in the DQD is mediated by the piezoelectricity of the wire material, while the phonon modes are defined by the bulk properties of the substrate. With the DQD length d = 120 nm, the wire radius a = 25 nm, the speed of sound in the wire c n = 4 × 10 3 m/s and in the substrate c p = 5 × 10 3 m/s [47], we find F piezo = 0.85, P piezo = 0.16.
Given that these estimates do not include any of the complexities of the real experimental system (e.g. nearby surface gates, shear coupling between the 1D wire and the bulk substrate, finite length of wire etc.) these values are in reasonable qualitative agreement with the values used in Fig. 3 in the main text, where F and P are treated as free parameters. The remaining microscopic parameter values used are summarized in Table C1, where we have corrected for angular averaging in 3D piezoelectric constants compared to 1D compounds.

Effect of higher temperature
To illustrate the effect of a raising the bath temperature, we compare in Fig. C1 the microwave power gain G at a bath temperature of T = 3 K (dashed red) and at T = 9 K (solid red). These calculations take into account all additional dissipators in Eq. (20). We see that the main effect is a broadening of the gain profile, leading to a better fit of the theory to experimental data (orange dots) in the region q /ω r < −5.  Figure C1. (Color online) Effect of raising the phonon bath temperature from 3 K to 9 K. We plot the power gain in the microwave resonator G as a function of DQD asymmetry q . The solid red line is for the elevated temperature, kBT /ωr = 23.4 (corresponding to T = 9 K), the dotted line is for comparison at the lower temperature shown in Fig. 3 (i.e. kBT /ωr = 7.8, corresponding to T = 3 K, from [3]). All other parameters are the same for both curves. Shown are calculations including all additional terms in the full master equation, Eq. (20). Parameters used in the calculations: ω d /ωr = 1, g/ωr = 0.0125, ∆q/ωr = 3, κ/ωr = 52 × 10 −6 , Γ/ωr = 0.34, F = 2.9, P = 0.25, and w = 1.7, as in Fig. 3.