Spectral functions and negative density of states of a driven-dissipative nonlinear quantum resonator

We study the spectral properties of Markovian driven-dissipative quantum systems, focusing on the nonlinear quantum van der Pol oscillator as a paradigmatic example. We discuss a generalized Lehmann representation, in which single-particle Green's functions are expressed in terms of the eigenstates and eigenvalues of the Liouvillian. Applying it to the quantum van der Pol oscillator, we find a wealth of phenomena that are not apparent in the steady-state density matrix alone. Unlike the steady state, the photonic spectral function has a strong dependence on interaction strength. Further, we find that the interplay of interaction and non-equilibrium effects can result in a surprising"negative density of states", associated with a negative temperature, and we identify two different mechanisms responsible for this phenomenon.


I. INTRODUCTION
Recent experimental progress in controllable quantum systems has renewed the interest in driven-dissipative quantum phenomena. Such systems typically have nontrivial, non-thermal steady states determined by the balancing of drive and dissipation. Examples include atomic and optical systems such as ultracold gases in optical lattices [1] or trapped ions [2], as well as solid state systems such as arrays of nonlinear superconducting microwave cavities [3][4][5] or microcavity polaritons realizing quantum fluids of light [6]. The simplest regime to consider is where dissipative effects are describable by a standard Markovian Lindblad master equation. Even here, considerable complexity can arise if there are interactions (nonlinearities). A vast amount of theoretical work has focused on finding (either exactly or approximately) the steady state of such systems, and the corresponding steady-state expectation values of observables [7][8][9][10][11].
While describing steady states is clearly of interest, many experimental probes involve studying how a system responds to a weak applied perturbation. One is then naturally interested in understanding the Green's functions that describe the linear response of the system to external perturbations. For Markovian systems, these correlations functions can be readily computed using the quantum regression theorem, and have been studied in a variety of different contexts, from the standard example of resonance fluorescence of a driven two-level atom [12][13][14], recently discussed in the case of arrays of coupled qubits [15], to the second-order correlations probing bunching/anti-bunching of time-delayed photons (see, e.g., [16]). The topic of correlation functions is also a standard topic in almost any quantum optics textbook (see, e.g., [17,18]).
Despite this existing work, methods for obtaining physical intuition from the behaviour of correlation functions (or their corresponding spectral functions) remain of interest. For closed, equilibrium quantum many body systems, the Lehmann representation [19,20] (see also, e.g., [21]) is a powerful tool. It expresses a single particle FIG. 1. Schematic plot of the setup considered in this manuscript. A cavity mode with Kerr nonlinearity is driven by an incoherent pump and subject to two-photons losses. We investigate its spectral features, encoded in the cavity mode spectral function. This quantity could be directly measured by considering the reflection of a weak probe tone. spectral function in terms of the system energy eigenstates, and allows one to interpret the spectral function in terms of Fermi Golden rule rates for the addition (or removal) of a particle. This directly connects to experimental probes (e.g. ARPES or tunneling spectroscopy), and is invaluable in constructing intuitive pictures.
In this work, we formulate the corresponding Lehmann representation of correlations functions of a driven dissipative system, and show that it also serves as a powerful interpretive tool. We focus on systems described by a Markovian Lindblad master equation, and discuss how the spectral functions are directly connected to the dynamical modes of the corresponding Liouvillian. As a concrete example, we analyze a simple but non-trivial model of driven-dissipative nonlinear quantum van der Pol oscillator, describing a single-mode bosonic cavity with a Kerr interaction subject to incoherent driving and nonlinear loss (see Ref. 22 for a comprehensive review). This model has recently received attention in the context of quantum synchronization [23][24][25]; it is also directly realizable in superconducting circuit architectures, where strong Kerr interactions and engineered two-photon losses have been experimentally achieved [26,27]. While the model has a relatively simple steady state, its spectral features are instead remarkably rich [22,28]. Unlike the steady state, the spectral function depends strongly on the size of the Kerr interaction, and reveals physics beyond that in the steady state density matrix. Specifically we show that the model features both population inversion in the density matrix and a negative density of states (NDoS), two aspects which are tightly connected in equilibrium but whose interplay in the drivendissipative case appears to be more complex. In particular we find a regime where NDoS emerges, even in absence of a population inversion in the stationary density matrix.
The paper is organized as follows. In Sec. II we define Greens' functions for an open system governed by a Lindblad master equation and discuss their decomposition in terms of eigenstates of the Liouvillian. In Sec. III we introduce the specific model of a driven-dissipative Kerr resonator, and in Sec. IV we discuss its spectral properties. We conclude in Sec. V.

A. Lindblad dynamics and Liouvillian spectrum
We consider an open quantum system described by the Lindblad master equation whereρ is the reduced density matrix of the system and L is the Lindblad/Liouvillian super-operator which takes the general form ( = 1) The first term on the RHS describes (unitary) Hamiltonian evolution, whereas the remaining terms describe incoherent driving and dissipative processes (each corresponding to a time-independent operatorL α ). Note that we will use throughout calligraphic letters to indicate super-operators. We assume the most common case where Eq. (1) has a unique time-independent stationary stateρ s . We are interested in calculating two-time correlation functions, which depend only on time differences due to the timetranslational invariance of the stationary state. Under the same assumptions one makes to derive Eq. (1), one can write dynamical correlators in terms of the Liouvillian, a result known as the quantum regression formulae [17,18]. For example, assuming t > 0, we have that whereÂ(s),B(s) are generic Heisenberg-picture operators of the system.
Recall that for closed systems in thermal equilibrium, it is extremely useful to relate Green's functions directly to the energy eigenvalues and eigenstates of the system; this is achieved by the Lehmann representation [29,30]. Our goal is to do something analogous in our driven dissipative system. Now however, the relevant spectrum is not that of the system HamiltonianĤ, but rather that of the Liouvillian. Such spectral decompositions for correlations functions for Lindblad open systems have been derived before (see e.g. [31,32] in the context of electron transport through correlated impurities or Ref. 22 in the contex of non-linear oscillators where a related decomposition as sum of partial spectra was introduced), but we include it here again for completeness and clarity.
To proceed, the first step is to properly enumerate the the eigenmodes of the Liouvillian. It is useful to treat the space of operators (so-called Liouville space) as a Hilbert space, with the inner product Â ,B ≡ tr Â †B [33]. This allows one to define the adjoint Liouvillian by Â ,L(B) = L † (Â),B . The left eigenvectors (l α ), right eigenvectors (r α ), and eigenvalues λ α of the Liouvillian are then defined viaL Note thatl α andr β are mutually orthogonal: tr l † αrβ = N δ α,β . For convenience, we fix the normalization constant N = 1, such that we have a simple completeness relation The eigenmodes of the Liouvillian directly determine how the system relaxes to the steady state. First, note that the unique (by assumption) steady state of our system corresponds to the unique right eigenstate ofL with zero eigenvalue. We label this eigenstate by the index α = 0, thusρ s =r 0 /tr [r 0 ] and λ 0 = 0. Suppose now that at t = 0 the system starts in some stateρ(0) that is not the stationary state. At later times, the system reduced density matrix will be given bŷ with c α = tr l † αρ (0) .
We can thus interpret each Liouvillian eigenmode α as a possible dynamical decay mode of some initial deviation from the steady state, with a decay rate given by −Re λ α . In general, a given decay mode will involve both diagonal elements of the density matrix in the energy-eigenstates basis (i.e. populations) as well as off-diagonal elements (i.e. coherences). However, in some cases the situation simplifies and one can cleanly separate the eigenmodes into processes only involving populations (T 1 processes) or only involving coherences (T 2 processes); we will see this explicitly in Sec. III A.

B. Spectral representation of correlation functions
To derive the Lehmann representation of the correlation function in Eq. (3), we note that the operatorB acting on steady state density matrixρ s causes the system to deviate from the steady state. Just as in Eq. (7), this deviation can be expressed as a linear combination of the Liouvillian decay modes (i.e. in terms of the righteigenstates of the Liouvillian). We thus obtain: At an intuitive level,B "excites" the various dynamical eigenmodes of the Liouvillian; these modes then oscillate and decay as a function of time. The factor involvinĝ A corresponds to the change in Â (compared to the steady state value) associated with exciting a particular dynamical eigenmode α. In what follows, we will focus on the connected average Â (t)B(0) − Â (t) B (0) ; equivalently, we will shiftÂ andB each by a constant so that they vanish in the steady state. In this case, the steady state (i.e. α = 0) does not contribute to the sum in Eq. (9). It is interesting to see how one recovers the standard closed-system Lehmann represetation by taking the zerodissipation limit of Eq. (3). This limit implies only keeping the Hamiltonian term in Eq. (1), i.e. replacingL with −i Ĥ , • (where • is the operator on whichL acts). Letting |ψ i and E i denote the eigenstates and eigenvalues of the HamiltonianĤ, it is straightforward to find the dynamical eigenmodes ofL. Each dynamical eigenmode α corresponds to a pair of energy eigenstates i, j: These modes have a simple interpretation. For a closed system, populations in the energy eigenstate basis are time-independent, corresponding to the zero-eigenvalue modes λ (0) i,i . Further, the coherences in the energy eigenstate basis have a simple undamped oscillatory behaviour, corresponding to the i = j modes. We stress that in the purely closed system case, the dynamics no longer picks out a unique steady state, as any incoherent mixture of energy eigenstates is stationary. The only constraint from the dynamics to gaurantee stationarity is thatρ s be diagonal in the energy eigenstate basis: ρ s = k p k |ψ k ψ k |. As usual, one must then assume the probabilities p k when computing average values and correlation functions. Formally, this assumption corresponds to the usual limit where the dissipation is nonzero but infinitesimally weak; in this limit, dissipation can determine the steady state, but does not impact dynamics.
Using the above eigenmodes in Eq. (9) and defining The first line matches what one would obtain from a direct calculation using Â (t)B(0) = tr e iĤtÂ e −iĤtBρ s . In the second line, we have used the diagonal form ofρ s . For a system in thermal equilbrium, the p k are simple Boltzmann weights; we then recover the usual textbook thermal equilibrium formula (see, e.g., [21,29]).

C. Single particle Green's functions
In the rest of the paper we will focus on the retarded single-particle Green's function of a bosonic system. Lettingâ denote the canonical bosonic annihilation operator, the single-particle Green's function is defined as This correlation function plays an important role in many different contexts. For example, via the Kubo formula, it describes the linear response of the expectation â(t) to a weak, classical field h(t ), which couples linearly toâ † . In the case whereâ describes a photonic cavity mode, G R (t) can be directly measured by weakly coupling the cavity to an input-output waveguide and measuring the reflection of a weak probe tone (see e.g. [34,35]).

Closed stationary system
For a closed system in a time-independent steady state, the Fourier transform of the Lehmann representation of the retarded Green's function is [21] where η is a positive infinitesimal. Consistent with causality, this function is analytic in the upper half plane. It has simple poles with infinitesimal negative imaginary part, and with purely real weights. Of particular interest is the imaginary part of G R (ω), which defines the singleparticle spectral function or density of states A(ω): For a closed system, the spectral function follows directly from Eq. (14): The first equality allows us to give a simple physical interpretation of A(ω) in terms of Golden rule transition rates. The first term is naturally associated with adding a particle to the steady state and creating an excitation with energy ω, whereas the second term is associated with removing a particle and creating an excitation with energy −ω.
The second equality in Eq. (16) also leads to an important result. If we assume that p j ≤ p i whenever E j ≥ E i , then we immediately can conclude: i.e. the sign of the spectral function A(ω) matches the sign of ω. A violation of this condition indicates the existence of population inversion in the steady state: a higher-energy energy eigenstate has a larger population in the steady state than a lower-energy eigenstate. While this is impossible in thermal equilibrium, it is indeed possible in a generic driven-dissipative non-thermal steady state. We will discuss the emergence of population inversion in the spectral function in great detail in the next section, in the context of a specific system.

Open system case
We can use the generalized Lehmann representation, Eq. (9), to derive a corresponding result for the singleparticle retarded Green's function of a Linbdlad open system. We obtain There is clearly some similarity to the closed-system expression Eq. (14). Like the closed-system case, the Green's function is decomposed into a sum of simple poles. However, whereas for the closed system poles occurred at energy differences that were infinitesimally shifted from the real axis, now the poles occur at eigenvalues of the Liouvillian, and will be shifted a finite distance below the real axis.
More intriguingly, the residues w α associated with the poles of G R (ω) are no longer necessarily real (as it must be for a closed system). This has a direct consequence on the spectral function (c.f. Eq. (15)), which now takes the form It follows that the spectral function is no longer simply a sum of Lorentzians. An immediate corollary is that unliked the closed-system case (c.f. Eq. (17)), the sign of the spectral function is not controlled in a simple way by the structure of the steady state distribution. In other words, a driven-dissipative systems can have spectral functions which violate the sign property Eq. (17), without this necessarily coming from an inverted population of the stationary state. We will see this explicitly in Sec. IV in the context of a driven quantum Kerr cavity. We remark that the spectral function in Eq. (19) satisfies the sum rules originating from the commutation relations of operators at equal time: as one can directly verify from Eq. (18). As a result, interactions, driving and dissipation can reshape A(ω), but they cannot change its area. Our spectral decomposition makes it clear that one can extract information on the eigenvalues of the Liouvillian from the frequency dependence of G R (ω). This could be particularly useful in extended systems, e.g., to determine a dissipative phase transition in which the eigenvalue λ α with smallest non-zero real part becomes purely imaginary [36][37][38]. We note that in general it is difficult to compute the spectrum of a Liouvillian, whereas the Green's function may be found via many body techniques (see, e.g., [39]).

D. Effective temperature
As we have seen, Green's functions (via the spectral function A[ω]) can provide information on the effective single-particle density of states of our system. They can also provide information on how these states are occupied, i.e. the effective distribution function or the effective temperature of our system. To obtain this information, one must consider the Keldysh single-particle Green's function: which heuristically describes the fluctuations of the observableâ. If the system was in true thermal equilibrium, the quantum fluctuation-dissipation theorem (FDT) would have required [30,40] where T is the system temperature.
In a non-equilibrium system, there is no well-defined temperature and the FDT does not hold in general. Nonetheless, it is useful to use the LHS of the FDT relation in Eq. (23) to define at each frequency an effective temperature T eff (ω), i.e.: As discussed extensively in Ref. [41], this T eff [ω] has a direct operational meaning and is a useful quantity in many different physical contexts (e.g. the theory of optomechanical cavity cooling using driven resonators [42]). In general, if a second narrow-bandwidth auxiliary system interacts weakly with our main system via exchanging photons, it will equilibrate to a temperature T eff (ω aux ), where ω aux is the frequency of the auxiliary system. As example, the auxiliary system could be a qubit with splitting frequency ω aux , which interacts with the main system via H int ∝ (σ +â + h.c.) [35,41]. Note that for a general non-equilibrium system, there is no requirement that the effective temperature T eff [ω] be positive. The fact that iG K (ω) > 0 implies that the sign of A(ω) dictates that of T eff (ω). In particular, if the sign of A(ω) obeys the equilibrium property Eq. (17), then T eff (ω) > 0 ∀ω. If this is not true, then there will be frequency regions in which T eff (ω) is negative. We thus see that the anomalous sign of the spectral function discussed earlier is directly connected to the existence of negative effective temperatures. We stress that this negative temperature has physical consequences. Again, consider weakly coupling an auxiliary qubit to our system. If the qubit splitting frequency ω aux is such that T eff (ω aux ) < 0, the qubit would thermalize at negative temperature, implying a population inversion (i.e. higher probability for the qubit to be in the excited state rather than its ground state).
Finally, we stress that in general, T eff [ω] for a nonequilibrium system will both be frequency dependent and operator dependent. That is, if one defined T eff [ω] using the FDT relation for Green's functions corresponding to operators other thanâ,â † , one would in general obtain a different function T eff [ω] [15,41].

III. APPLICATION TO THE DRIVEN-DISSIPATIVE VAN DER POL OSCILLATOR
We now use the general results of the previous section to study a specific, non-trivial driven dissipative system. We consider the nonlinear version of the well-known quantum van-der Pol oscillator [22,25]: a bosonic mode with a Kerr nonlinearity subject to incoherent singleparticle driving, and two particle losses. It is described by the master equation ( = 1) Here ω 0 is the cavity frequency and U/2 the strength of the Kerr (or Hubbard) interaction. γ is the two-photon loss rate, while γr is the single photon pumping rate. We will set ω 0 = 0 in the following, as it can be eliminated by moving to a rotating frame. Note first that the unique steady state density matrix of this model has been previously found analytically in Ref. [28,43]. The steady state is an incoherent mixture of photon number Fock states; further, it is completely independent of the interaction strength U , and is only determined by the dimensionless parameter r (ratio of the driving to the nonlinear loss). The photon-number probabilities in the steady state are in fact determined by a classical master equation (i.e. coherences play no role). In Fig. 2 we plot the photon number probabilities p n in the steady state for two different values of r. For the smaller value of r, the probabilities decay monotonically with n, whereas for large values, one obtains a peaked, non-monotonic distribution. As the Hamiltonian H dictates that energy increases with increasing photon number, this latter situation corresponds formally to a population inversion.
A natural question to now ask is whether this inversion effect (which is essentially classical) manifests itself in the cavity's spectral properties. Such a question was first addressed in a series of seminal works [22,28], where the spectral properties of a related quantum van der Pol oscil-lator were discussed (see instead Ref. 44 for the undriven model). More recently, the power spectrum of a coherently driven quantum van der Pol oscillator was computed to investigate signatures of synchronization [24].

A. Liouvillian eigenmodes and symmetry considerations
To understand the Green's functions of our model, it will be useful to first discuss its symmetry properties. Due to driving and dissipation, the system does not conserve photon number. Nonetheless, the Liouvillian is invariant under the U (1) symmetryâ →âe iθ . This implies that the LiouvillianL commutes with the superoperator K = [â †â , •] that generates the symmetry operation. As a result, the eigenvalues k ofK are quantum numbers which label the eigenstates ofL. We can use this to write the Liouvillian in the block-diagonal formL = ⊗ kLk , wherê L k acts only within the eigensubspace ofK corresponding to the (integer) eigenvalue k. We denote the right eigenstates of a particular blockL k bŷ In Fock space, we see that this is a matrix that only has non-zero elements along the kth off-diagonal. The presence of this symmetry greatly reduces the numerical complexity of the problem, as we can diagonalize the different blocks separately. It also gives a simple physical way to label the different eigenmodes ofL. Eigenmodes corresponding to k = 0 describe how diagonal elements of the density matrix (in the Fock basis) decay. Such decay modes conventionally referred to as T 1 relaxation processes. In contrast, eigenmodes corresponding to k = 0 describe how Fock-state coherences decay. These are generically referred to as T 2 relaxation processes.
Note crucially that different correlation functions will only be sensitive to a particular (small) subset of Liouvillian eigenmodes. For example, for the single particle Green's function defined in Eq. (13), it is only the eigenmodes corresponding to k = 1 that contribute. This follows immediately from Eq. (18) and the fact that for k = 1: Analogously, correlation functions like G (2) = â † (t)â † (t)â(0)â(0) would probe Liouvillian eigenmodes with k = 2, i.e. T 2 processes involving coherences between states whose photon number differs by 2. Similarly, a correlation function destroying k bosons at t = 0 and creating k bosons at time t would probe the decay of coherences between states whose photon numbers differ by k. It also follows that if one wishes to probe T 1 processes (i.e. k = 0), one needs to look at density-density correlation functions.

IV. SPECTRAL PROPERTIES OF DRIVEN-DISSIPATIVE VDP OSCILLATOR
We now turn to the spectral properties of the nonlinear driven-dissipative cavity model introduced in Eq. (25). We use the Lehmann representation given in Eq. (18) to compute the spectral functions numerically, by truncating the bosonic Hilbert space to a maximum number of states N max = 15 and diagonalizingL. The cutoff N max = 15 is enough to obtain accurate results, as it is shown by the agreement of the steady state numerical solution with the analytical prediction in Fig. 2. We further checked that the results for the Green's functions are stable by increasing N max and that they satisfy sum properties like Eq. (21).

A. Spectral function and the role of interactions
In Fig. 3 we plot the spectral function A(ω) (c.f. Eq. (15)) of our system for several values of the dimensionless interaction strength U/γ and for two values of the drive/loss ratio r. An immediate result, visible in all four panels, is that the spectral functions strongly depend on the interaction strength. This dependence is remarkably different from the steady state density matrix, which (as discussed in Sec. III) is completely insensitive to U . Heuristically, while the steady state density matrix is completely independent of the system's coherent Hamiltonian dynamics, the system's response to perturbations retains a strong dependence onĤ.
For a more detailed analysis, consider first the regime of relatively weak driving where r = 0.5 (top row of Fig. 3). To understand lineshapes, recall from Sec. III A that the spectral function is probing T 2 decay modes which describe the decay of coherences between Fock states |n and |n + 1 . The oscillation frequency of these coherences is largely determined by the coherent Hamil-tonianĤ. For U = 0, there is no Hamiltonian, and coherences do not oscillate; we thus obtain a single peak in the spectral function. As U/γ is increased, distinct peaks become visible in A(ω) (each approximately Lorentzian), corresponding to different coherences and different decay modes; the peaks become more and more resolved with increasing U/γ. Note that in this weak driving regime, there is no obvious signature of non-equilibrium in the spectral function.
For larger values of the driving parameter r (bottom row of Fig. 3, r = 6), the situation is markedly different. For large driving and large enough interaction U , we find that the spectral function hits zero at a positive finite frequency, and for larger frequencies, becomes negative. We term this negativity of A(ω) at ω > 0 a "negative density of states" (NDoS). This is a clear indicator of non-equilibrium: as discussed in Sec. (II C 1), this cannot happen in a system that is in thermal equilibrium. We also stress that (as discussed in Sec. II D and in [35]), this NDoS corresponds to a negative effective temperature; this is shown in Fig. 4. As also discussed, this negative temperature effect could be directly probed by coupling the cavity weakly to an auxiliary probe qubit.
One might first think that the NDoS effect here is simply a reflection of the population inversion in the steady state photon number distribution, which occurs when r is sufficiently large. This is not the case: while the population inversion in the steady state is independent of U/γ, A(ω) only becomes negative at ω > 0 for sufficiently large U/γ. This is shown explicitly in Fig. 3. The relation between the NDoS effect in the spectral function and population inversion in the steady state is thus not entirely trivial; we will explore this in more detail in the next sections. Note that similar spectral function negativity in presence of a population inversion has previously been identified in a related model of a quantum van der Pol oscillator in presence of negative damping and monochromatic drive [22,28], as well as in parametrically driven bosonic systems [35].

B. Dissipation-Induced Lifetime
The results of the previous section show that the spectral properties of the nonlinear quantum VdP oscillator are remarkably rich. In this section, we investigate the extent to which these can be understood using a perturbative approach where the only dynamical effect of dissipation and driving taken into account is to give a finite lifetime to the Fock-state eigenstates of the system Hamiltonian H.
Our starting point is the open-system Lehmann representation of Eq. (18). We will approximate the eigenstates of the Liouvillian to be the same as those of the closed system, e.g. simple outer products of Fock states (c.f. Eq. 10), ignoring their perturbative corrections. We will however take into account the modification of the Liouvillian's eigenvalues to leading order in the dissipation (i.e. in γ). Formally, this procedure can be implemented using the Lindblad perturbation theory approach introduced in Ref. [45]. We write our full Liouvillian asL =L (0) +D, with the unperturbed LiouvillianL (0) [ρ] = −i[Ĥ,ρ] and the perturba- . WhileL (0) is highly degenerate, one can still employ the simple non-degenerate perturbation theory of Ref. [45] due to the symmetry ofL discussed in Sec. III A. This symmetry prevents mixing between different degenerate sectors.
We expand eigenvalues and eigenvectors of the Liouvillian in powers of γ: with the unperturbed quantities λ α already defined in Eqs. (11), (10). We retain the perturbative corrections to the eigenvalues, while ignoring for the time being any corrections to the eigenstates. The validity of such an approximation and the role of these corrections will be discussed later in the manuscript. Perturbation theory tells us that the leading order correction to the Liouvillian eigenvalues λ α are given by λ α ) . Accordingly, Eq. (18) yields the following approximate form for the spectral function: where Note that the first order correction to the λ α is purely real, implying there is no shift in the position of the spectral function resonances. The approximate spectral funciton in Eq. (31) is exactly the same as the equilibrium expression in Eq. (14), except that the populations p n are non-thermal, and each resonance has a finite width Γ n+1,n . While we have shown how this width can be calculated using formal perturbation theory, it also has a simple physical origin: it is the sum of the Fermi's Golden rule decay rates for the states |n and |n + 1 , i.e.
The first line is the decay rate due to the incoherent driving, the second due to the two-photon loss.
In Fig. 5 we compare the perturbative result with the full calculation obtained with the Lehmann representation for r = 6 and two values of the Kerr interaction. We see that at large U/γ 100 the perturbative approach captures rather well the main features of the spectrum, in particular the location of the peaks, their width and weight. However, upon decreasing the interaction, the agreement deteriorates, as we show for U/γ = 3. This behaviour is of course not surprising, as the perturbative approach is only valid in the small dissipation limit 1 U/γ,r U/γ. As a rough rule of thumb, when resonances in Eq. (31) begin to overlap, perturbation theory starts getting bad, as the spacing between adjacent resonances is of order U and the width of the resonances of order γ.
Taking into account the dissipation-induced lifetime in Eq. (31) allows to uncover a mechanism by which dissipation can mask the effect of a population-inverted density matrix on the spectral function. Indeed, a population inversion in the stationary state, if there were no lifetime broadening, would certainly result in a violation of the Green's functions sign property in Eq. (17), as one can see straight from Eq. (14). On the other hand, the lifetime broadens the resonances, making them overlap and possibly resulting in those with smaller weights to be completely masked by bigger ones. As a result, the spectral function in Eq. (31) does not obey anymore a precise sign rule which is dictated by the behaviour of populations of the density matrix. As a corollary, the presence of population inversion in the stationary state may not be revealed by a change of sign of the spectral function.
In Fig. 6, we summarize the above analysis by presenting a phase diagram, in the (r, U/γ) plane, the regions of parameter space exhibiting the NDoS effect (i.e. spectral function A(ω) is negative at positive frequencies). The region r > 1 (shaded grey) indicates where the steady-state exhibits a population inversion; this boundary can be determined analytically from the exact stationary state solution [43] and we remark that it is independent of U . In contrast, the spectral function is sensitive to both interaction and non equilibrium effects, resulting in a non-trivial value U c (r) above which the negative density of state emerges. We plot this threshold interaction strength both for the numerically exact calculation of the spectral function (red-solid line), and for the approximate perturbative (lifetime broadening) calculation (red-dashed line). In general, the perturbative approach underestimates the NDoS effect; further, it fails to yield any population inversion in the region r < 1. In contrast, the numerically exact calculation reveals that NDoS can occur even for r < 1, i.e. in regions where the steady state photon number exhibits no population inversion. This is a remarkable result, which points toward yet another origin of NDoS, as we are going to further discuss below.

C. Dissipative effects beyond lifetime broadening
As demonstrated above, the simple (perturbative) lifetime broadening of eigenstates introduced in Eq. (31) was able to capture many aspects of the spectral function of our model. It however failed to describe the most interesting aspect of Fig. 6: there are parameter regions where the spectral function exhibits NDoS, even though the steady state density matrix does not exhibit population inversion. As we now show, this effect can also be captured in perturbation theory if we go beyond simply calculating a correction to the Liouvillian eigenvalues due to dissipation, but also calculate the change to the eigenmodes themselves. The leading eigenmode correction can cause the weight factors w α in Eqs. (18) and to acquire an imaginary part, implying that the spectral function is no longer a simple sum of Lorentzians. This provides a new route for NDoS.
Using the same perturbation theory used in Sec. IV B, we can analytically compute the leading-order-in-γ correction to the Liouvillian eigenmodes. Following [45] and using the existence ofL −1 , the first order corrections to the right and left eigenstates are given by: As expected, dissipation mixes the various eigenmodes together with a strength that is inversely proportional to the difference in eigenvalues. Here, the denominator is purely imaginary (as all unperturbed eigenvalues are imaginary). As discussed, for the spectral function, the unperturbed modes of interest correspond to coherences between the |n and |n + 1 Fock states: With dissipation, these modes acquire a real part to their eigenvalues, corresponding to dephasing. The first order correction to the mode wavefunctions take the form: At a physical level, these corrections tell us that dephasing eigemodes of the Liouvillian no longer correspond to a single Fock state coherence; rather, each mode involves three distinct coherences.
In Fig. 7 we show the effect of including these eigenmode corrections in the evaluation of the spectral function. We see that this modified approach is able to capture non-Lorentzian contributions to the spectral function, and to improve qualitatively and quantitatively the agreement with the exact numerical result. In particular, a region of negative density of states now appears at small frequency, an effect which is completely missed by the lifetime broadening approximation.
These eigenmode corrections can also be given a physical interpretation in terms of interference of different dephasing modes. Consider the contributions to the timedomain correlation function in Eq. (9) associated with a particular initial photon number m: n e λn+1,nt l l|âr n+1,n |l m|l † n+1,nâ † |m ρ m,m (40) Recall the interpretation: starting with m photons, we add a photon to the cavity, exciting a dephasing eigen-mode α = (n + 1, n) of the Liouvillian. To 0-th order in dissipation, the time-independent weight factors are necessarily real. This follows from the fact that i) r (0) n+1,n =l (0) n+1,n , and ii) the only non-zero contribution is when l = n = m, i.e. adding a photon to |m excites a single, unique dephasing eigenmode.
Including dissipation to first order, both conditions (i) and (ii) no longer hold. In particular, as the dephasing eigenmodes no longer correspond to a single Fock coherence (c.f. Eq. (38)-(39)), adding a photon to |m can simultaneously excite several distinct dephasing eigenmodes. It is the interference between these processes that give rise to complex weights and hence non-Lorentzian contributions to the spectral functions. The spectral function is thus sensitive to an interference in the dynamics, even though there is no coherence in the steady state density matrix.
Stepping back, we thus see that even for weak dissipation, the spectral function is sensitive to more than just the lifetime-broadening effect of dissipation: the fact that dissipation can also create more complicated dephasing processes also directly impacts the form of A(ω). This gives rise to anti-Lorentzian contributions, and (in our model) negative density of states in regimes where the steady state exhibits no population inversion.

V. CONCLUSIONS
In this work we have studied the spectral properties of driven-dissipative quantum systems, taking the simple case of a quantum van der Pol oscillator as a working example. We have first derived some general results concerning the single particle Green's function of systems described by a Lindblad Master Equation. Using a decomposition in terms of exact eigenstates of the Liouvillian we have derived a Lehmann like representation for the Green's function and compared it to the well known result for closed systems in thermal equilibrium. Such a result, in addition of being of practical relevance for numerical computations whenever the system is sufficiently small to be diagonalized exactly, has also a conceptual value. From one side it connects properties of the Liouville eigenvalues and eigenstates, which are of theoretical interest but often hard to access, to the behavior of the spectral functions, which are of direct experimental relevance. In addition it allows for a more transparent interpretation of spectral features in regimes far from equilibrium, for which a simple intuition is often lacking or misleading. As an example we have shown that the well known sign property of equilibrium Green's functions, changing sign at zero frequency as a result of thermal occupation, can be violated in driven-dissipative systems and it is in general not directly constrained by the structure of the stationary state density matrix.
We have then applied our approach to the case of a Kerr nonlinear oscillator with incoherent driving and two-particle losses. Such a model turns out to be a per-fect case study, since the properties of its stationary density matrix are well known, while its spectral features reveal a number of surprises. In particular the resonator density of state shows a strong dependence from the strength of the Kerr nonlinearity, a feature completely absent in the steady state populations only set by pump/loss ratio. Even more interestingly, in the regime of large interaction and large non-equilibrium imbalance a NDoS emerges, an effect which would not be possible in thermal equilibrium.
We have summarized the behavior of the spectral function of this model in the phase diagram of Figure (6) which shows that NDoS is not necessarily related to an inverted population in the steady state density matrix. In order to build physical intuition and to better understand the origin of this result we have developed a semianalytical approach that starts from the spectral function of the isolated problem and adds a lifetime due to dissipation in the spirit of a Fermi Golden Rule. This method, which turns out to be equivalent to a perturbation theory in the dissipation where only the eigenvalues of the Liouvillian are corrected, was able to partially capture the NDoS effect, at least for sufficiently large interaction and whenever the stationary density matrix shows population inversion. Finally we have shown that including the perturbative correction to the eigenstates of the Liouvillian results into a new mechanism for NDoS, due to the emergence of complex weights in the spectral function. This turns to be crucial to capture NDoS in the regime where the populations of the steady state are not yet inverted.
To conclude we mention that the approach outlined here is rather general and can be used to shed light on the spectral properties of other small driven-dissipative quantum models. Interesting future directions include for example the study of resonance fluorescence lineshapes beyond the two-level system limit [46,47], the spectral features of a coherently driven cavity across a zero-dimensional dissipative phase transition [48][49][50] or applications related to quantum synchronization [51,52].