Third quantization of open quantum systems: new dissipative symmetries and connections to phase-space and Keldysh field theory formulations

The connections between standard theoretical tools used to study open quantum systems can sometimes seem opaque. Whether it is a Lindblad master equation, the equation of motion for the Wigner function or a dissipative Keldysh action, features evident in one formalism are often masked in another. Here, we reformulate the technique of third quantization in a way that explicitly connects all three methods. We first show that our formulation reveals a fundamental dissipative symmetry present in all quadratic bosonic or fermionic Lindbladians. This symmetry can then be used to easily diagonalize these models, and provides a intuitive way to demonstrate the separation of dissipation and fluctations in linear systems. For bosons, we then show that the Wigner function and the characteristic function can be thought of as ''wavefunctions'' of the density matrix in the eigenbasis of the third-quantized superoperators we introduce. The field-theory representation of the time-evolution operator in this basis is then the Keldysh path integral. To highlight the utility of our approach, we apply our version of third quantization to a dissipative non-linear oscillator, and use it to obtain new exact results.


I. INTRODUCTION
Open quantum systems are generically much more difficult to characterize than their closed-system counterparts.Quadratic systems of bosons or fermions coupled to Markovian baths are arguably an exception: computing average values and correlation functions is straightforward using either the master equation, Keldysh field theory or the Heisenberg-Langevin equations.However, it is often useful to understand the full structure of the eigenvalues and eigenvectors of the Lindbladian.This information can help answer questions that are difficult to address using Keldysh or Langevin equations, e.g. the full time evolution of an arbitrarily complicated non-Gaussian initial state.In a set of seminal works, Prosen [1] and Prosen and Seligman [2] introduced the technique of third quantization, which allows one to obtain the spectral decomposition of quadratic multi-mode Lindbladians in a manner that is analogous to diagonalizing a second-quantized quadratic Hamiltonian.For all its virtues and widespread use in the literature [3][4][5][6][7][8][9][10][11], the method as presented does not provide an intuitive picture of open quantum system dynamics.Further, it is not a priori clear how it is related to more standard methods.A more physically-transparent formulation of third quantization which addresses these issues could thus help make it an even more powerful tool.
In this work we present an alternate reformulation of third quantization which satisfies both criteria simultaneously.This is achieved by introducing a new set of canonical superoperators, defined in Eq. (50) for bosons and Eq.(74) for fermions.In this basis, one directly identifies the two pieces of data that determine any quadratic Lindbladian: a non-Hermitian Hamiltonian which acts as an effective dynamical matrix, and a second matrix which encodes the noise.This separation of dissipation and fluctuations further reveals a fundamental symmetry of all such models.We show that this symmetry can be used to effectively gauge away the fluctuation terms in the master equation using a novel similarity transformation implemented at the superoperator level.This provides a simple quantum-oriented way to demonstrate that the eigenvalues of the Lindbladian are independent of noise in linear systems.
The superoperators we introduce also make the connection between the Lindbladian and the Keldysh formalism evident.We make this notion exact by formulating a path integral representation of third quantization over a finite time contour.For bosons, we go further by demonstrating that third quantization is naturally linked to phase-space representations of the density matrix.Namely, we show that the Wigner function and its characteristic function can be thought of as wavefunctions in the basis of what is essentially the thirdquantized equivalent of position and momentum eigenkets.We use this insight to study a dissipative interacting model, a thermally-damped Kerr cavity.Our approach lets us analytically describe the time-evolution of cavity's Wigner function, starting from an arbitrary initial state (see Fig. 1 for an example).To the best of our knowledge, this result has not been previously presented, and highlights the power of our approach.
The paper is organized as follows.In Sec.II we work out in detail the diagonalization procedure of the Lindbladian describing a thermally-damped harmonic oscillator using our formulation of third quantization.We then discuss how our approach is naturally connected to both phase-space and Keldysh formalisms of open quantum systems in Secs.III and IV respectively.The insight gained from this example allows for a straightforward generalization to arbitrary quadratic bosonic Lindbladians as shown in Sec.V. Despite the change in statistics,  103).With no dissipation κ = nth = 0, the resulting state is a superposition of coherent states [16].The formation of such a state is however very sensitive to any amount of dissipation.As explained in Sec.VII, third quantization allows us to obtain the exact evolution of the Wigner function for an arbitrary initial state.a similar approach is presented for the fermionic version of this problem as detailed in Sec.VI.We finish by applying our phase-space formulation of third quantization to the dissipative Kerr oscillator problem in Sec.VII.
Before we proceed, let us highlight the differences between this work and the work of Prosen and Prosen and Seligman.The main dissimilarity is our starting point: we use a distinct set of superoperators to express the third-quantized Lindbladian, see Eq. (50) and Eq.(74).Although a seemingly-minor difference, this alternate choice of basis makes the symmetry we use to gauge away fluctuations apparent.Our superoperators also provide a natural connection between thirdquantization, Keldysh field theory and phase-space representations of open quantum systems, which has not been considered in previous works [1,2,[12][13][14][15].

II. DIAGONALIZING A LIOUVILLIAN USING THIRD QUANTIZATION
A. Setup While the diagonalization method we discuss can be applied to arbitrary multi-mode bosonic quadratic Lindbladians, for clarity we discuss these ideas in the simplest possible setting: a harmonic oscillator coupled to a ther-mal Markovian bath.Our starting point is the equation of motion for the density matrix (with = 1 throughout) where â is the bosonic annihilation operator of the oscillator, ω 0 is its frequency, κ its decay rate and nth is the thermal occupation of the bath.Here, we have defined the usual dissipator as D[ X]ρ = X ρ X † − { X † X, ρ}/2 and L as the Liouvillian superoperator.Our goal is to diagonalize L. It is convenient to think of operators as being elements of a Hilbert space ρ → |ρ with the usual Hilbert-Schmidt inner product [17] X| Ŷ ≡ Tr X † Ŷ (2) for arbitrary operators X and Ŷ [18].Using this notation, superoperators become operators acting on this new space and we thus write them with hats.To differentiate them from operators acting on the Hilbert space of wavefunctions, we write all superoperators using a bold typeface L → L. The inner product Eq.( 2) then allows us to define the adjoint of any superoperator.For instance the adjoint Liouvillian, which controls the time evolution of observables, is in standard second quantized form given by As the Linbladian is non-Hermitian L = L † , to diagonalize it we must find both its right and left eigenvectors.That is, we seek operators and complex numbers which satisfy where, with some foresight, we have labelled the eigenvectors by two non-negative integers µ, ν ≥ 0.
Although the physical interpretation of the right and left eigenvectors may not be straightforward, the eigenvalues have a simple intuitive meaning.The real part of E µ,ν plays the role of energy differences between eigenstates, whereas the imaginary component dictates the rate of decay.For the linear system of interest however, just as in the classical case, the characteristic oscillation and decay rates are independent of temperature.This simply reflects the fact that a linear system's response to noise is always independent of the noise itself, something that is clear e.g. by using Langevin equations to describe our system.We thus conclude that the eigenvalues of L must not depend on nth .A crucial question though still remains: how do we rigorously identify these fluctuations and reach this conclusion directly from the master equation?

B. Classical and Quantum Superoperators
We argue that this can be straightforwardly done by introducing appropriate third-quantized superoperators defined by with which we refer to as classical and quantum superoperators.The naming convention is in analogy with Keldysh field theory, which is justified by writing the Lindbladian in this basis of superoperators This form directly mimics the structure of the Keldysh action for this same model [19,20]; we make this connection more explicit in Sec IV.Just as in the field theory description, the coefficients of â † q âcl and âq â † cl serve as an effective non-Hermitian Hamiltonian, whereas the term â † q âq should be thought of as the fluctuations that must accompany the dissipation [21].Although â † cl and â † q are the conjugate of âcl and âq , they are not the creation operators of true bosonic modes since [â cl , â † cl ] = [â q , â † q ] = 0. Thus, we cannot interpret â † cl âcl and â † q âq as number operators.Rather, despite being non-Hermitian, we should think of â † q âcl and −â q â † cl as operators whose eigenvalues correspond to a number of quanta in a mode.This follows since the only non-vanishing commutation relation are with 1 the identity superoperator.Given their status as number operators, to determine the eigenvectors of â † q âcl and −â q â † cl we must first find the "vacuums" of âcl , â † cl and âq , â † q .These are simply the parity and identity operator which follows since e iπâ † â and 1 respectively anticommute and commute with both â and â † .Using Schur's lemma, these vacuums are unique.In complete analogy with a simple harmonic oscillator, we can obtain the simultaneous eigenvectors of these number operators by repeatedly applying the appropriate creation operators on the vacuums.Thus, we find that are the right and left eigenvectors of â † q âcl and âq â † cl with eigenvalue µ and −ν respectively.
The above eigenvectors are unfortunately not eigenvectors of the Liouvillian: while the first two terms in Eq. ( 6) are proportional to our generalized number operators, the last term is not.We next show that there is a simply way to deal with this.

C. Gauging away the noise
Recall that we have already reasoned that fluctuations in this linear system should never affect the eigenvalues of L, only its eigenvectors.We should thus be able to remove the temperature-dependent noise term −iκ(2n th + 1)â † q âq entirely using a similarity transformation.This can indeed be achieved by defining and using the commutation relations Eq. ( 7) along with the Baker-Campbell-Hausdorff identity to verify that Since this superoperator and L are isospectral, we confirm that as promised, the temperature does not effect the eigenvalues of the Lindbladian.Further, the transformed Lindbladian is now just a sum of our generalized number operators, and hence is diagonalized by the eigenvectors introduced in Eqs.(9a)-(9b).Note that the transformed Lindbladian is not in Lindblad form due to the presence of a negative rate which is however irrelevant for diagonalization purposes [22].
We have thus made a non-trivial mapping at the level of master equations from a noisy damped harmonic oscillator to its fluctuation-free equivalent.Our ability to remove the noise term exactly is a consequence of the usual separation of dynamics and fluctuations in linear systems, yet here there are quantum consequences, seeing as how the eigenmodes of the Liouvillian have a Fock state structure.
We pause to note that one could imagine adding a quadratic term to our Liouvillian that would break our symmetry, i.e. a term proportional to â † cl âcl .The similarity transformation e −(2n th +1)â † q âq would no longer bring L to a diagonal form, since the former commutes with the quantum superoperators but not the classical ones.Such a term is however forbidden by the trace-preserving nature of the dynamics Tr(Lρ) = 0q | L|ρ = 0, which a classical-classical term would violate.At a quantum level, being able to eliminate the noise can be understood as a special kind of symmetry; the combination of linearity and conservation of probability.This structure and an analogous transformation to Eq. (10) will also be present in the multi-mode bosonic case considered in Sec.V. A comparable symmetry, transformation and conclusion will also be presented for fermionic master equations in Sec.VI.In non-linear quantum systems, which are of course also probabilityconserving, fluctuations can impact the dynamics as we show explicitly in Sec.VII D. Eigenvalues, eigenvectors and non-Hermitian quasiparticles With Eq. ( 11) already in diagonal form, the similarity transformation allows us to simply relate the eigenvectors rµ,ν , lµ,ν of L to r µ,ν , l µ,ν , those of â † q âcl and âq â † cl .They are given by with corresponding eigenvalues Here we have identified ρss as the unique zero-eigenvalue eigenvector of L, i.e. the steady state of dissipative dynamics.We have also defined where note that since V is not unitary, â † cl (n th ) is not the conjugate of âcl (n th ) despite the notation.This is however of no concern, seeing as the similarity transformations preserve all commutation relations, e.g.
We can thus identify âcl (n th ), âq and â † cl (n th ), âq as the correct independent "non-Hermitian quasiparticle" superoperators which create and destroy excitations on top of the their vacuum with a definite energy and decay rate.Although the vacuum of the quantum superoperators is the identity, the right vacuum of these new classical superoperators is the steady state which is of course explicitly temperature-dependent.

E. Remaining unsatisfactory issues
This completes the diagonalization of a single thermally-damped oscillator (something that can be easily generalized to the multi-mode case, as we show below).There are however two unsatisfactory features of the procedure presented above.First, the similarity transformation V at this stage does not seem to possess a simple physical interpretation.For instance, it is well-known that the steady state is a Gaussian state, with an average particle number of nth .In contrast, our procedure simply asserts that it is equivalent to e −(2n th +1)â † q âq | 0cl .Is there a direct way to see the connection here, that is It is at this stage not clear how we should think about this result or even explicitly demonstrate this equality.Second, the right eigenvectors are given by nested commutators of â † , â and the steady state; their explicit secondquantized form is tedious to obtain and largely uninformative.Take for example the right eigenvector describing population decay with a characteristic decay rate of 2κ, r2,2 = 1 8n 2 th (n th + 1) 2 n2 − (4n th + 1)n + 2n 2 th ρss (20) with n ≡ â † â.At first glance, there is no obvious intuitive way to think about such an eigenvector, let alone a simple connection to a 2-excitation Fock state.
Table I.The action of third-quantized creation and annihilation superoperators on the density matrix |ρ , Wigner function W ρ(α) or characteristic function Λ ρ(η).These rules can be used to obtain an equivalent set of equations of motion, see Eq. ( 31) for the single linear oscillator case.

III. EQUIVALENCE BETWEEN THIRD-QUANTIZATION AND PHASE-SPACE FORMALISM
While seemingly devoid of physical content, we shall now demonstrate that both the similarity transformation V in Eq. ( 10) and eigenvectors rµ,ν , lµ,ν of Eq. ( 13) have transparent phase-space representations.The continuous nature of phase space arises by asking a natural question: what are the eigenvectors and eigenvalues of the classical and quantum superoperators?As we now show, the answer to this question leads to the main results of this section: the similarity transformation V is equivalent to a phase-space convolution with a Gaussian of width 2n th + 1 and the eigenvectors of L have a functional form equivalent to that of a Fock state.

A. Classical and quantum eigenvectors
Thankfully, the eigenvectors of âcl , â † cl and âq , â † q are already familiar from quantum optics.For any complex number α and η let us define corresponding displaced parity and displacement operators [23] with D(α) ≡ e αâ † −α * â.Using the defining property of the displacement operator â D(α) = D(α)(â + α) and parity operator âe iπâ † â = −e iπâ † ââ, we arrive at The spectrum of each of our third-quantization superoperators is thus the whole complex plane.The corresponding eigenvectors are reminiscent of standard single-mode coherent states, an analogy we can make even stronger by writing in the same way that coherent states are displaced vacuum.Unlike coherent states however, which are not eigenvectors of â † , âcl and â † cl share the same eigenvectors, as do âq and â † q .This a consequence of [â cl , â † cl ] = [â q , â † q ] = 0.For the same reason, αcl and ηq are orthogonal to other eigenstates of the same kind: ηq where the trace was computed using coherent states.In this sense, these operators are more reminiscent to position and momentum eigenkets.Just like their firstquantized counterparts, they can then be used to form a resolution of the identity where d 2 β ≡ dReβ dImβ for any complex β.It follows that we can write the density matrix as with Remarkably, these "wavefunctions" of the density matrix in these bases are well-known objects: W ρ(α) and Λ ρ(η) are precisely the Wigner function of ρ and its characteristic function respectively [24,25].
To recover the result that W ρ(α) and Λ ρ(η) are Fourier transforms of each other, we compute the overlap use the definitions Eqs.(27a)-(27b) and the resolution of the identify Eq. ( 25) to obtain Much like position and momentum eigenkets corresponds to a perfectly spatially-localized or delocalized particle in standard quantum mechanics, the eigenkets |α cl and |η q can be thought of as states in phase space which are perfectly localized at α or delocalized with wavevector η.That α and η are Fourier transform pairs is also consistent with âcl , â † q and âq , â † cl being conjugate to one another.This can be further demonstrated by us-ing Eqs.(23a)-(23b) to write Given any third-quantized Lindbladian, we can then use Eqs.(27a)-(27b) and Eqs.(30a)-(30b) to readily find the equation of motion for either the Wigner or the characteristic function.These rules are summarized in Table I.For example, we can recover the standard equivalent equations of motion [25] i∂ From now on, we will write the classical and quantum superoperators and their action on these eigenvectors as equivalent, e.g.âcl = α = ∂ η * .

B. Interpretation of V and phase-space eigenvectors of L
It is worth stressing that by using the resolution of the identity Eq. ( 25), one can obtain a phase-space representation of an arbitrary operator ) X .This is nothing but the Wigner-Weyl transform [26], yet here it appears much more naturally using the language of standard quantum mechanics.This allows us to easily describe how the similarity transformation V which eliminated the noise acts on an arbitrary density matrix, something which would be extremely challenging using standard formulations and tools of quantum mechanics in phase space, e.g.phase-point operators [27,28].Being diagonal in the quantum basis, V = e −(2n th +1)â † q âq simply acts by multiplication on the characteristic function Λ We can then obtain the Wigner function by Fourier transforming using Eq.(29a) As promised, V is phase-space convolution with a Gaussian of width (2n th + 1), spreading out and delocalizing any distribution.This is also consistent with the similarity transformation e −(2n th +1)â † q âq removing the noise, or equivalently the diffusion term in the Fokker-Planck equation Eq. (31b), which has an identical delocalization effect.Note that if nth = 0, then this is precisely the Glauber-Sudarshan P -function of the oscillator, which is insensitive to vacuum noise [25].
Turning our attention to the right eigenvectors of L, we find that they have a compact form in phase space Focusing on the µ = ν case for simplicity and comparing with the Wigner function of Fock state projectors |µ µ| we obtain where L µ is the µ-th Laguerre polynomial.Unlike the explicit second-quantized form, which from Eq. (13a) requires computing nested commuators of â, â † and ρss , in this representation rµ,ν bear a striking similarity to Fock states.Although they can be made almost identical by re-scaling the phase-space variable W rµ,µ (α) → W rµ,µ ( √ 2n th + 1α), which would make the Gaussian factors match, the argument of the Laguerre polynomial would still differ by a factor of 2. This prop-erty also holds when comparing arbitrary right eigenmodes rµ,ν and outer products of Fock states |µ ν| which we show in Appendix A. There we also give the phasespace representations of the left eigenvectors.
Using the equivalent representations of the master equation Eq. ( 31)), we conclude that the phase-space representations of |r µ,ν and lµ,ν | are also right and left eigenvectors of the classical Fokker-Planck differential operator.The existence of quantized Laguerre-polynomial eigenfunctions is a well-known property of this classical damped oscillator problem [29], but the corresponding quantum analog arguably is not.The quantized nature of these modes in a classical or quantum setting is easily understood for this linear problem: any observable should only oscillate and decay at integer multiples of the fundamental frequency and decay rate respectively.For instance, (â † ) ν âµ and its classical counterpart includes the harmonics ω 0 (µ − ν) and decay κ 2 (µ + ν).Finally, note that these eigenvectors are still correct even when there is no dissipation nth → 0, κ → 0. This might seem surprising at first, as in this limit W rµ,µ (α) = W |µ µ| (α).This discrepancy is easily understood.Without decay, each eigenvalue There is thus no unique way to diagonalize L; in this limit our approach picks out one linearlyindependent set of eigenvectors out of the infinitely-many available ones.

IV. CORRESPONDENCE BETWEEN THIRD QUANTIZATION AND KELDYSH FIELD THEORY
Despite the Lindbladian Eq. ( 6) being reminiscent of the dissipative Keldysh action of the same model, the connection between these approaches we have been suggesting are at this stage not rigorous.In this section, we make the correspondence concrete by showing that Keldysh field theory emerges directly from the structure of the third-quantized Lindbladian expressed in the basis of classical and quantum superoperators.Our approach here complements standard derivations of the Keldysh action for open systems [19,20].Unlike these treatments, our focus is not the derivation of a partition function for the calculation of correlation functions, but rather on obtaining the time evolution of an arbitrary initial state.In this sense, this section introduces an alternate formulation of the Feynman-Vernon influence functional approach [30] and subsequent closely-related descriptions [31,32].
A. Path integral for a thermally-damped harmonic oscillator The Keldysh path integral emerges here by asking how the Wigner or characteristic function evolves in time.Indeed, since e −i Lt is the time evolution superoperator, we can use the resolution of the identity Eq. ( 25) formed using the classical and quantum operators to write where The matrix element K(η, α; t), also referred to as the kernel, thus serves as the mixed phase-space propagator, and in principle contains the same information as the spectral decomposition of L. Using a resolution of the identity which follows directly from Eq. ( 25) and Eq. ( 28), along with the defining property of the classical and quantum eigenvectors, we can obtain a path-integral description of Eq. (36b) using the familiar prescription [33].The central object will be a functional very similar to the standard Keldysh action of a thermally-damped harmonic oscillator.There are however important differences that emerge.We emphasize these here, leaving the technical and mostly-familiar details in Appendix B. Once these details have been taken into account, we are left with a path integral representation of the kernel where D[a cl , a q ] is the usual functional integration measure, the action takes the form and the inverse of the matrix Green's function reads Throughout this work, we will denote vectors and matrices with a blackboard bold font to distinguish them from operators and superoperators.Here we have defined the row vectors a † (t ) ≡ a * cl (t ), a * q (t ) and the source term Note that in most path integrals, from Eq. (36b), one would normally impose boundary conditions on the fields such as a cl (0) = α, a q (t) = η; this is incorrect here as it implies that G −1 has a zero-eigenvalue eigenvector.We stress that one must instead keep track of the boundary term through the source J and work directly the with the discrete representation of the Green's function as we do in Appendix B.
There we show that, although G −1 is symbolically the same regardless of the length of the time contour, the matrix Green's function is sensitive to such changes.While the retarded G R (t , t ) and advanced G A (t , t ) Green's function remains the same regardless of the contour, the Keldysh Green's with Θ(t) the Heaviside step function.Crucially, the discrete-time path integral makes it clear that the correct choice of the Heaviside step function at zero is Θ(0) = 1.
The path integral can be computed in the usual manner when one has a quadratic action and a source term.Making the displacement and, using D[a cl , a q ]e iS0[a cl ,aq] = det(iG) = 1 as in the usual Keldysh formalism, we obtain As expected due to the quadratic nature of the Lindbladian, the propagator is a Gaussian function of α and η.Equation ( 45) can independently be verified to be correct by noting that it must satisfy Eq. (31c) along with the initial condition K(η, α; 0) = e η * α−c.c. .

B. Correlation functions and timer evolution of the Wigner function
For the sake of completeness, in this subsection we demonstrate how, with knowledge of K(η, α; t) and the initial Wigner function W ρ(0) (α), one can straightforwardly obtain both arbitrary correlation functions and the time-evolved Wigner function W ρ(t) (α).Since the characteristic function Λ ρ(t) (η) is the generator of symmetrically-ordered correlation functions, we can simply use Eq.(36a) to obtain which is valid for any Lindbladian, and not simply the damped oscillator considered here.
Using the resolution of the identity Eq. ( 25) we can obtain the classical-classical phase-space propagator which we can then use to write Once again this is valid for an arbitrary Lindbladian, but can be done here exactly due to the Gaussian nature of the problem.

C. Path integral for arbitrary Lindbladians
Although in this example L was quadratic in creation and annihilation superoperators, it is evident that a path integral representation of the propagator can be obtained for arbitrary Lindbladians.Replacing all âq , â † q and âcl , â † cl by their respective fields in the action however is only valid once all quantum superoperators have been placed to the left of the classical ones.This is due to the form of the resolution of the identity Eq. ( 25).This convention is different than the one used in the standard coherent-state path integral, where one must ensure that the Lindlabidan is normal-ordered [19,20].No differences emerge in this linear setting, but the two different prescriptions will produce two different Keldysh actions for an interacting problem.We expand on this point in Sec.VII when we study a dissipative Kerr oscillator.

A. Setup
We now generalize the previous results to an arbitrary quadratic multi-mode bosonic Lindbladian.The model now under consideration is a Lindblad master equation of the form where n and m label independent bosonic modes [â m , â † n ] = δ n,m .Further H nm and K nm are, respectively, arbitrary Hermitian and complex symmetric matrices which describes the coherent particle-conserving and non-conserving process of the isolated system.The coefficients l bm , p * bm characterize the coupling to a set of dissipative baths indexed by b and, consequently, the positive semi-definite Hermitan matrices L nm = b l * bn l bm and P nm = b p * bn p bm capture the phase-insensitive loss and pumping respectively.The complex matrix C nm = b l * bn p * bm describes coherences between these processes.Its phase-sensitivity indicates that not all quadratures are equivalent.

B. Multi-mode classical and quantum eigenvectors
Mimicking the procedure of the single-oscillator case, for each mode m we define a set of classical and quantum creation and annihilation superoperators whose only non-vanishing commutators are A tedious but straightforward calculation gives the thirdquantization representation of the Lindbladian Here we have defined the effective Hamiltonian and twophoton pumping matrices in addition to the phase-insensitive and phase-sensitive noise matrices The form of L directly parallels the single-mode example.The quantum-classical terms encodes the dynamics, the quantum-quantum terms the noise, and there is no classical-classical term due to conservation of probability.
The presence of the matrices K and C breaks a global U (1) symmetry âm → e iθ âm , â † n → e −iθ â † n and, just like in the closed-system case, implies that we have to work with a Nambu structure.For clarity, we first set K = C = 0 and avoid the additional technical details that comes with doubling the size of the matrix one has to diagonalize.At the end of this section we briefly discuss how the presence of these terms changes the diagonalziation procedure.

C. Correlation matrix
To justify referring to H eff and N as the dynamical and noise matrices respecticely, let us compute the equations of motion for the symmetrically-ordered correlation matrix which is given by As promised, H eff serves as the dynamical matrix which determines the fluctuation-free evolution of the system.It is also the dynamical matrix that would appear in the equation of motion for the average values of âm (t).
The noise matrix N only affects the correlations and occupation of the modes.Using the Heisenberg-Langevin equations, it can be shown that N corresponds to the correlation matrix of a set of fluctuating noisy forces acting on the oscillators.Further, note that the dissipation cannot be larger than the noise i(H eff − H † eff ) ≤ N, which ensures that the Heisenberg uncertainty principle is satisfied [21].Finally, we remind the reader that the identification of noise is always contingent on the ordering prescription used to define the corelation matrix.For example, the equations of motion for the normal-ordered covariance matrix takes the same form as that of S but with modified noise terms: N → P.
Let us assume that all eigenvalues of H eff have negative imaginary part, ensuring the existence of a unique steady state.The formal solution to Eq. ( 56) then reads where the steady state correlation matrix can be written as Here, are the frequency-space retarded and advanced Green's function respectively whereas σ and τ label the biorthonormal eigenvectors and eigenvalues of H eff which satisfy It is worth stressing that while Eq.(58c) assumes H eff can be diagonalized, Eqs.(58a)-(58b) are always valid.Thus, even when the dynamical matrix H eff is defective or equivalently at an exceptional point [34,35], the steady state is well defined.

D. Gauging away noise, eigenvectors and eigenvalues
Returning to the question of how to diagonalize L, recall that in the single-mode case linearity and the lack of a classical-classical term is ultimately what let us easily diagonalize the Lindbladian.As we now show, a similar approach works in the multi-mode case.We can again find a multi-mode similarity transformation that effectively gauges away the fluctuations, leaving a transformed Liouvillian that can be written in terms of generalized number operators.The generalization of e −(2n th +1)â † q âq to several oscillators is simply The Liouvillian in the new frame defined by V is, using the Baker-Cambell-Hausdorff identity and Eq.(51), where we have used H eff S ss − S ss H † eff + iN = 0 which follows from Eq. (56).In this new basis, we are then left with two commuting third-quantized non-Hermitian Hamiltonians.Using the same eigenvectors and eigenvalues we used to diagonalize H eff , these two terms can easily be brought to diagonal form (i.e.written as a sum of commuting generalized number operators).Moving back to the original frame, we finally arrive at where the relevant quasiparticle superoperators are with Ψ R m (σ) ≡ m|Ψ R σ and (Ψ L m (σ)) * ≡ Ψ L σ |m .Although â † cl (σ) is not the conjugate of âcl (σ), this is irrelevant since, using the bi-orthonormality of the eigenvectors of H eff , the only non-vanishing commutators amongst the set Eqs. (64a)-(64d) are We can therefore think of â † q (σ)â cl (σ) and −â q (σ)â † cl (σ) as non-Hermitian number operators.The right and left vacuums of both number operators is simply the steady state and identity operator respectively with 0q ≡ 1.The steady state |ρ ss = V| 0cl is Gaussian with two-point correlation matrix determined by S ss .Here 0cl = 2 M e iπ N is the total parity operator with N = m â † m âm , which serves as the joint vacuum of the classical superoperators, and M is the total number of modes.
In complete analogy to a second-quantized quadratic Hamiltonian, we can now write the eigenvectors and eigenvalues of L as with corresponding eigenvalues Here, µ, ν are lists of non-negative integers µ σ , ν σ which characterize the number of "particles" µ σ or "holes" ν σ in mode σ.The right and left eigenvectors are normalized such that the bi-orthonormality condition l µ , ν |r µ, ν = δ µ, µ δ ν, ν is satisfied.

E. Multi-mode Wigner function and path integral
The equivalence between third quantization, phasespace representations of the density matrix and Keldysh field theory presented in Secs.II-IV can be extended to this multi-mode setting in an obvious manner.For any set of complex numbers α m and η m the multi-mode displaced parity and displacement operator are the eigenvectors of the classical and quantum superoperators respectively.The multi-mode Wigner function and characteristic function are given by from which we can obtain a simple physical interpretation of the similarity transformation used to eliminate the noise.Using the basis which diagonalizes the Hermitian covariance matrix S ss , that is the set of orthogonal modes âλ satisfying {â λ , â † λ } ss = (2n λ +1)δ λλ , we have V = exp − λ (2n λ + 1)â † q,λ âq,λ .Given the result of Sec.III we conclude that V is a convolution with a multimode Gaussian in phase spaces along a set of orthogonal axes determined by the modes âλ whose width depends on the average number of quanta in that mode in steady state.
The equations of motion for W ρ( α) and Λ ρ( η) can be easily read off from the third-quantized form of L in Eq. ( 52) and the multi-mode generalization of the rules in Table I e.g.âcl,m = α m = ∂ η * m .We do not reproduce them here for the sake of compactness.We can also find the propagator K( η, α; t) ≡ ηq |e −i Lt |α cl for this differential equation using the Keldysh path integral.The multi-mode matrix Green's function for the finite-time contour of length t is nothing but the generalization of Eq. ( 42) and consequently the propagator is nothing but the multi-mode generalization of Eq. ( 45) where here we have defined the integrated source term

F. Symmetry-breaking terms
In the presence of U (1) symmetry-breaking terms in the Lindbladian Eq. ( 52) due to non-zero K and C, the strategy would be analogous to the one presented above.
One would now have non-vanishing anomalous steadystate covariance matrix elements {â m , ân } ss , which would now enter the similarity transformation used to remove the noise terms.What remains is a 4M × 4M non-Hermitian Hamiltonian.Just as in the closed-system case, one would have to perform an appropriate Bogoliubov transformation when defining third-quantized quasi-particle superoperators.The propagator would also take essentially the same form, the only difference being that the matrix Green's function would be of size 4M × 4M .

A. Setup
In this section we will diagonalize an arbitrary fermionic Lindbladian of the form where m and n label independent fermionic degrees of freedom {ĉ m , ĉ † n } = δ nm .The matrices H, L, P and C have the same form and interpretation as in the bosonic case studied in Sec.V.The only difference is that the complex matrix K now describes the pairing of fermions, and is thus necessarily anti-symmetric.

B. Creation and annihilation superoperators
In analogy with the bosonic case Eq. ( 50), we would like to define a set of fermionic superoperators which satisfy canonical anti-commutation relations.The issue however is that superoperators acting on the left and the right of the density matrix commute with one another; this would seem to preclude introducing superoperators with canonical fermionic anti-commutation relations.The solution is to introduce the parity operator e iπ N in the definition of these modes, where N = m ĉ † m ĉm is the total number operator.Defining one verifies that ĉ † 1,m and ĉ † 2,m are the Hermitian conjugates of ĉ1,m and ĉ2,m as per the definition of the inner product Eq.( 2).The only non-vanishing set of anticommutators are Thus, ĉ1,m and ĉ2,n can be interpreted as genuine annihilation operators of fermionic modes.
From the definitions Eqs (74a)-(74b), the total parity and identity operator serve as the unique joint vacuum for one type of superoperator, while simultaneously serving as the completely-filled state of the other, i.e.where M is the total number of modes.
We stress that these superoperators and states differ from Prosen [1], who instead introduced Majorana superoperators instead of fermionic creation and annihilation superoperators.Our choice allows us to make the connection to the Keldysh formalism, as is apparent when writing the Lindbladian in this basis where for fermions the effective non-Hermitian Hamiltonians and pairing matrices are defined as with the noise matrices taking the form Setting K = Q = 0 for simplicity, we confirm that the nomenclature is justified by defining the antisymmetrized covariance matrix 80) and computing its equation of motion Note that here the dissipation must be larger than the fluctuations i(H eff − H † eff ) ≥ N, whereas the opposite is true for bosons.This ensures that the Pauli exclusion principle is satisfied.The solution and physical content of Eq. ( 81) was already discussed in Sec.V. Our only goal here was to stress that each term in the Lindbladian has a physical interpretation; it describes either dynamics or noise.

C. Gauging away the noise
Just as in the bosonic case, there are terms in the third-quantized representation of L such as ĉ † 2,n ĉ1,m which are not allowed due to conservation of probability 1| L = 01 , 12 | L = 0. We can use this to move to a frame where there are no fluctuations, interpreting the separation of dissipation and noise in this linear system as the underlying symmetry which allows us find such a basis.Defining the anti-symmetrized steady-state covariance matrix (A ss ) mn ≡ [ĉ m , ĉ † n ] ss which satisfies we can use the Baker-Campbell-Hausdorff formlua and the anti-commutations relations Eq. ( 75) to obtain We can now use the eigenvectors and eigenvalues of H eff to bring this superoperator to diagonal form (i.e.written as a sum of commuting generalized number operators).
Moving back to the original frame, we have where we have defined Although ĉ † 1 (σ) and ĉ † 2 (σ) are not the conjugate of ĉ1 (σ) and ĉ2 (σ), the canonical anti-commutation relations are still preserved, e.g., However, as a result of the non-unitary nature of V, the left and right vacuum and completely filled state of these new superoperators are not the same where we can write |ρ ss = 2 −M V| 01 , 12 .In direct analogy with a standard Hermitian Hamiltonian, the eigenvectors of L can then be written as with corresponding eigenvalues Here, µ, ν are lists with entries µ σ , ν σ = {0, 1} and the eigenvectors satisfy l µ , ν |r µ, ν = δ µ, µ δ ν, ν .The inclusion of U (1) symmetry-breaking terms adds no formal degree of complexity; one simply has to diagonalize a larger matrix.The underlying physics, that is the separation of dynamics and fluctuations, still allows one to eliminate the noise using an appropriate similarity transformation analogous to Eq. (82).

D. Third quantization and Keldysh field theory for fermions
Just as we did in the bosonic case, we now develop a path integral representation of the dynamics starting from the third-quantized representation of the Lindbladian in Eq. (77).For simplicity, let us assume that we have a single fermionic mode, such that the thirdquantized Lindbladian takes the form where γ is the decay rate and n = ĉ † ĉ ss the average number of particles in the steady state.Although here the Hamiltonian loss and pumping matrices have been simply replaced by numbers, H → 0 L → γ(1 − n), P → γ n, the generalization to several modes is straightforward.The first step in formulating any fermionic path integral is to introduce the fermionic equivalent of coherent states, for which one must work with the usual Grassmann numbers [33].We thus introduce four Grassmann numbers ψ 1 , ψ 2 , ψ1 , ψ2 which anticommute with each other and all superoperators.Recalling that we have identified the identity operator as the vacuum of type-1 quasiparticlces and the occupied state of type-2 quasiparticles, see Eq. (76), we then define which satisfy and One can then verify that these coherent states can be used as a resolution of the identity from which one can construst a Grassmann-valued "wavefunction" which, in analogy with the Wigner function for bosons, we denote W ρ(ψ) ≡ ψ|ρ .It follows that the matrix element contains all information about the dynamics.A continuous-time version of this propagator can be build using essentially the same procedure presented in Sec.IV.Leaving all details to Appendix C, we arrive at where D[ψ 1 , ψ 2 ] is the standard functional measure, and the action is where the matrix Green's function now reads The matrix Green's function is then where Mimicking the bosonic case, we can make the appropriate translation to the Grassmann variables to eliminate the source term and integrate over the quadratic action which gives unity (see Appendix C for details).We then obtain From Eq. ( 92), by taking derivatives of the Grassmann variables we can use this propagator to obtain expectation values of observables, e.g.
Of course, we could have obtained this answer easily in this singlemode example; the application of this path integral is more useful when there are several fermionic modes.

VII. APPLICATION: DISSIPATIVE NON-LINEAR CAVITY
We now provide an example of how our formalism can be extremely useful even in situations where there are true interactions, and the Lindbladian is not simply quadratic.We consider a paradigmatic problem: a single bosonic mode with a Kerr (or Hubbard) type nonlinearity, and subject to thermal single-photon dissipation.The master equation for this system is with U the strength of the Kerr non-linearity.Since L is now evidently quartic in creation and annihilation superoperators, it would seem as though describing this problem exactly is not feasible.
Surprisingly, well-known analytic expressions for this system exist, derived using a variety of different techniques (see e.g.[36][37][38][39][40][41][42]).While elegant, these solutions typically provide little direct intuition or insight into the underlying physics.Further, using them for practical calculation is often complicated (e.g.Ref. [43] for instance chose to study Eq. ( 103) numerically as opposed to exploiting existing analytic expressions).A more recent work Ref. [44] demonstrated that the solvability of this interacting model is related to a special symmetry.However, the results of that work still do not provide a simple prescription for calculating observable quantities, especially if one is interested in a phase-space picture of the dynamics.
In this section, we will use the equivalence between third-quantization and the Keldysh path integral to show that there exists simple physical picture of the interplay between the non-linearity and damping in this system.We call this "self-dephasing": the fluctuating photon number in the oscillator causes it to dephase itself, i.e. degrade Fock space coherences.Further, using a novel gauge transformation within the path integral, this intuition will lead to simple expressions for physical observables of interest.Namely, we will demonstrate that we can obtain the time evolution of the Wigner function for an arbitrary initial state.We note that we will also implicitely make use of the results in Secs.III-IV extensively.

A. Propagator via Keldysh path integral
In our formulation of third quantization, the nonlinear term in Eq. ( 103) takes the form which, using Table I, in phase space is equivalent to where φ = arg α.The first term in the bracket tells us that the frequency of the mode is amplitude-dependent, whereas the second stems from the quantum noise associated with inherent uncertainty associated with the amplitude |α| 2 .The last term is due the ordering convention: phase-space averages weighted by the Wigner function are equivalent to symmetrically-ordered operator averages [25].Our Hamiltonian on the other-hand is normal-ordered, hence this extra factor.
A path integral representation of the mixed phasespace kernel K(η, α; t) can immediately be read off from Eq. ( 6) and Eq.(105), where the average is with respect to the Keldysh action .
Here we have at times dropped the temporal argument for notational compactness.We have also defined n cl (t ) ≡ a * cl (t )a cl (t ) + a * q (t )a q (t ).The nomenclature is justified by noting that this term comes from the contribution â † cl âcl + â † q âq = {n sym , •} in third quantization where nsym ≡ 1 2 {â, â † } is the symmetrized number operator.Despite the name, our formulation still takes into account all quantum fluctuations.Only by dropping the quantum noise term a * cl (t )a cl (t ) + a * q (t )a q (t ) → a * cl (t )a cl (t ) would the dynamics be rendered classical, since the phase-space equation of motion would be equivalent to that of a classical Fokker-Planck equation, see Eq. (105).
It is worth emphasizing that Eq. ( 107) is the usual Keldysh action for this model written in a non-standard way in that we have singled-out n cl (t ).This way of writing the action is useful however since it makes it evident that the Kerr non-linearity serves as an amplitudedependent frequency.By making the following gauge transformation within the path integral we effectively move to a frame where the numberdependent phase accumulated by the fields has been taken into account.This transformation has two effects: it both removes the non-linearity in the action S → S 0 and introduces a multiplicative phase factor attached to a cl (t) Although a seemingly a more unwieldy expression than our starting one Eq.(106), we now show how this can not can not only be straightforwardly evaluated but has a simple physical interpretation.

B. Self-dephasing through amplitude fluctuations
To remove the phase factor e −i U In going from Eq. (110a) to Eq.(110b), we have lumped the phase proportional to the integral of n cl (t ) in a new quadratic action defined by As shown in Appendix VII, the resulting Gaussian integral for each l can then be performed and we arrive at where we have defined φ = arg α, ϕ = arg η, J l as the l-th Bessel function of the first kind and with Γ l ≡ κ 2 − U 2 l 2 + 2iκU l(2n th + 1).It is worth emphasizing that for l = 0 Eq. ( 111) is not a valid Keldysh action, as it does not vanish when the quantum fields are set to zero [19,20].We should not expect this to be the case however, since this action has been modified to include the effects of the fluctuating particle-dependent phase and is only physically meaningful when computing specific averages.For example, if we recall that the wavefunction in the quantum basis Λ ρ(t) (η) is the generator of symmetrically-ordered correlation functions, we have We thus have a compact expression for the exact evolution of the average photon amplitude, starting with an arbitrary initial state.It can also be used to compute unequal-time correlation functions such as e.g.linear response coefficients [44,45].This result can be directly extended to more general averages.The general prescription is that the relevant quadratic action that one needs to use has a form that explicitly depends on the chosen correlation function.For example, for a correlation function like {â µ (â † ) ν } sym the corresponding action would be S U l 2 with l = µ − ν.We see that our analytically exact results here directly validate the interpretation that the Kerr nonlineary generates a fluctuating frequency which causes dephasing.In Fig. 2 we plot a scaled version of â(t) for various parameters with ρ(0) a coherent state of varying magnitude.We see that increasing either nth or the magnitude of the initial state tends to kill the revival of the average.This is in agreement with our intuition: increasing the noise or |α 0 | causes an increase in particle-number fluctuations and consequently an increase in dephasing.Unlike the linear problem presented in Sec.(II), the temperature changes both average values and in addition has an impact on the oscillation and decay rates.

C. Time-dependent Wigner function
Our third-quantization fueled approach for deriving exact expressions for this nonlinear dissipative cavity problem is useful even beyond calculating specific timedependent averages or correlators.As we now demonstrate, it allows us to analytically describe the time evo- Increasing the thermal occupation nth or the magnitude |α0| of the initial coherent state decreases the amplitude of the revivals at times tU = 2πn for n ∈ N.This is due to an increase in particle number fluctuations and thus an accompanying increase in dephasing.lution of the full Wigner function describing the state.From Eqs. ( 47)-(48) we have access to this information by Fourier-transforming in the η variable to obtain the propagator relating the Wigner function at time t = 0 to that at later times.We give this expression in Appendix D; it has the same functional form as Eq.(113).As an example, let us assume that the initial state is a coherent state ρ(0 The time-evolved Wigner function can be computed using only Gaussian integrals.Leaving these details to Appendix D, we have where Φ = arg α 0 and we have defined To the best of our knowledge, this relatively compact expression was not previously known.The clos-est analogue were results derived in Refs.[40,41] that were more complicated, involving two infinite summations.Other equivalent exact solutions instead work with equally-intricate expressions involving the evolution of Fock states [36]. In Fig. 1 we plot this function for various choice of parameters.We stress that our choice of initial state was purely for illustrative purposes, and one can obtain the Wigner function of ρ(t) for an arbitrary initial state with the propagator given in Eq. (D16)

VIII. CONCLUSION
By starting with a different set of superoperators for both bosons Eq. ( 50) and fermions Eq. ( 74), we have reinterpreted the standard approach to third quantization introduced by Prosen [1] and Prosen and Seligman [2].These new superoperators allows us to naturally relate third quantization, Keldysh field theory and the phasespace formulation of quantum mechanics.Further, our approach let us straightforwardly identify a symmetry of all quadratic Lindbladians, which can then be used to effectively gauge away fluctuations in these models.This leads to a simple and intuitive diagonalization procedure, and provides a quantum-inspired way of demonstrate that dynamics are not affected by noise in linear systems.Finally, we have shown that our formalism provides a simple picture of the dynamics in a paradigmatic dissipative non-linear cavity model, which can be used to provide straightforward exact expressions.
In future work, it would be worth exploring if third quantization can be extended to spins in a useful manner.Further, it is interesting to ask whether the set of non-Hermitian quasiparticles which were used to diagonalize the Lindbladian can be used to build a mean-field theory of interacting open quantum systems.
With the Baker-Campbell-Hausdorff identity, we can move e η * âcl +ηâ † cl past the factor of e (2n th +1)â † q âq and use the defining property of the classical and quantum eigenvectors to obtain In this Appendix, we provide the explicit computations which led us to the continuous-time representation of the propagator K(η, α; t) = ηq |e −i Lt |α cl in Eq. (38).We start with the usual approach [33] and divide the propagator in N − 1 parts e −i Lt = e −i L∆t • • • e −i L∆t with ∆t = t/(N − 1).We then insert the resolution of the identity Eq. ( 37) between each element e −i L∆t and ηq | and |α cl .Expanding to linear order in ∆t and reexponentiating gives âq,j |e where L(a q,j , a cl,j−1 ) ≡ âq,j | L|â cl,j−1 .Using the defining property of the classical and quantum eigenvectors Eqs.(22a)-(22b), L(a q,j , a cl,j−1 ) is obtained by simply replacing the classical and quantum superoperators that appear in L by a cl,j−1 and a q,j respectively.We thus have where the tilde indicates that we are working with discrete-time objects.We have defined the two 2N column vectors where [ G−1 ] R and [ G−1 ] A are Hermitian conjugate N ×N lower and upper-triangular matrices respectively, whose only non-vanishing elements are on the main off-diagonal h j (h * ) j × h − min(j,j ) (h * ) − min(j,j ) − 1 . (B12) The strange delta function terms (1 − δ j1 δ j 1 ) simply ensures that the Keldysh component vanishes when either argument is 1, i.e. at the lower boundary.This is automatically satisfied in the continuum representation where G K (0, t ) = G K (t , 0) = 0.In Eq. (B2), we can then make the displacement ã → ã + i GJ (B13) and perform the integral over the various a cl,j and a q,j using the well-known result for Gaussian inte- In this appendix, we demonstrate how to achieve the continuous-time representation of the path integral of Eq. (97).The procedure is nearly identical to the bosonic case presented in Appendix B and we will thus be rather brief in our derivation.
After writing the propagator e −i Lt = e −i∆t L • • • e −i∆t L as a product of N − 1 terms ∆t = t/(N − 1), using the resolution of the identity formed by the coherent states Eq. (95), expanding to linear order in ∆t, using the defining property of coherent states Eqs.(93a)-(93b) and exponentiating we arrive at a Gaussian integral whose discrete action is nearly identical to the bosonic version Eq. (B2) with Grasmmann variables taking the place of complex variables.If we define the 4N Grassmann variables to be integrated over with two indices ψ1/2,j , ψ1/2,j , then the main difference between the bosonic expression for the propagator is that the source term is e ψ1 ψ 1,N − ψ2,N ψ 2 + ψ where once again the retarded, advanced and Keldysh Green's function are the same as the bosonic version in Eqs.(B10)-(B12) with the aforementioned substitutions to the decay frate and occupation factor.The continuous-time version of the path integral is then given by Eq. (97) in the main text.
In the discrete representation, after eliminating the source term, one can integrate over the remaining Gaussian Grassmann integral using a well-known identity [33], after which one concludes that this term is the identity.The final result Eq. (102) then follows.where Y > and Y < are two matrices that satisfy We thus have four degrees of freedom remaining, as expected from a solution to a first order differential equation of a 2 × 2 matrix, which must be fixed by the appropriate boundary conditions.Instead of going back to the discrete representation of the action to fix these boundary conditions, we can instead make use of Dyson's equation Computing the matrix exponential only requires diagonalizing a 2 × 2 matrix and thus can be done exactly.With G U l 2 (t , t ) in hand, we can then eliminate the source term in Eq. (D1) as usual by making the appropriate displacement.
After performing the displacement, we must still compute the source-free functional integral, which is equal to the inverse of the functional determinant of the matrix Green's function [19]  where recall n cl (t ) ≡ a * cl (t )a cl (t ) + a * q (t )a q (t ).Taking the derivative with respect to U on both sides and using standard Gaussian integral identities we find

Fig. 1 .
Fig. 1.Exact Wigner function of a dissipative non-linear cavity mode at time t = πU (with U the Kerr constant), where in each case the initial state is a coherent state ρ(0) = |α α|, α = √ 2(1 + i).Evolution is generated by the master equation Eq. (103).With no dissipation κ = nth = 0, the resulting state is a superposition of coherent states[16].The formation of such a state is however very sensitive to any amount of dissipation.As explained in Sec.VII, third quantization allows us to obtain the exact evolution of the Wigner function for an arbitrary initial state.

2 e
iã † G−1 ã = det i G[33].Since the lower-right block of G vanishes and because the retarded and advanced Green's functions are lower and upper triangular respectively, we have det i G = det GR det GA = 1.We are then left withK(η, α; t) = exp −i J † GJ + O(∆t) = exp −i |η| 2 ( GK ) N N − η * α( GR ) N 1 − ηα * ( GA ) 1N + O(∆t) (B14)Using Eqs.(B10)-(B12) and taking the N → ∞ limit while keeping the length of the contour fixed ∆t(N − 1) = t, we recover Eq. (45).As we have emphasized in the main text, Eqs.(B10)-(B11) indicate that the correct normalization of the Heaviside step function at time t = 0 as they appear in the retarded and advanced Green's function is Θ(0) = 1.This stems from the ordering of the classical and quantum variables within the path integral; in the discrete representation one always takes the classical field one time step before the quantum fields.Appendix C: Continuous-time representation of the fermionic path integral -Details

2 (D1) where S U l 2 is 2 G U l 2 ( 2 ω 0 + i κ 2 ω 0 − i κ 2 U l 2 − 2 (2 Y > e it D U l 2 X+ 2 X
Appendix D: Path integral approach to the non-linear Kerr oscillatorIn this appendix, we provide all the technical steps to obtain the results in Sec.(VII).Let us first show that Eq. (110b) and Eq.(112) are equivalent by computing exp η * e −iθ a cl (t ) + αa * q (0) − c.c. S U l the non-standard quadratic action defined in Eq. (111).It follows that to compute the functional integral, we must first find the matrix Green's function G U l 2 (t , t ) corresponding to the action S U l 2 which by definition satisfiesi∂ t X − D U l t , t ) = δ(t − t ) (D2)where X is the first Pauli matrix iκ(2n th + 1) t , t ) = − iΘ(t − t )e −it XD U l iΘ(t − t )e −it XD U l 2 Y < e it D U l (D5)

d dU e −i U l 2 t0 2 t0 2 t 0 dt Tr G U l 2 (
dt n cl (t ) S0 e −i U l dt n cl (t ) S0 = l t , t ) (D14) To be explicit, | 11 , 02 defines a state where all type-1 quasiparticle states are filled, and all type-2 quasiparticle states are empty; | 01 , 12 is interpreted in a similar fashion.Note that these effective states are not normalized to unity but rather 11 , 02 | 11 , 02 = 01 , 12 | 01 , 12 = 2 M