# Field Theory of Many-Body Lindbladian Dynamics

Foster Thompson1  Alex Kamenev1,2
1 School of Physics and Astronomy, University of Minnesota, Minneapolis,
55414, MN, USA
2 William I. Fine Theoretical Physics Institute, University of Minnesota,
Minneapolis, 55414, MN, USA

###### Abstract

We review and further develop the Keldysh functional integral technique for
the study of Lindbladian evolution of many-body driven-dissipative quantum
systems. A systematic and pedagogical account of the dynamics of generic
bosonic and fermionic Lindbladians is presented. Our particular emphasis is on
unique properties of the stationary distribution function, determined by the
Lyapunov equation. This framework is applied to study examples of Lindbladian
dynamics in the context of band theory, disorder, collisionless collective
modes, and mean-field theory.

###### Contents

  1. 1 Introduction
  2. 2 Formalism
    1. 2.1 Lindblad and Keldysh
    2. 2.2 Bosons
    3. 2.3 Lyapunov Equation
    4. 2.4 Observables and Response
    5. 2.5 Exceptional Points
    6. 2.6 Fermions
  3. 3 Examples
    1. 3.1 Parametrically Driven Oscillator
    2. 3.2 Non-Hermitian Band Theory
    3. 3.3 Disordered Fermions
    4. 3.4 Lindbladian Gas
    5. 3.5 Mean Field Theory
  4. 4 Conclusion
  5. A Keldysh-Nambu Diagonalization
  6. B Superoperator Quantization
  7. C Disconnected Diagrams

##  1 Introduction

Non equilibrium dynamics of open quantum systems has attracted a lot of
attention in recent years [1, 2, 3, 4, 5]. The interest is mostly stimulated
by a rapid progress in design and manufacturing of prototypical quantum
computers with dozens and hundreds of qubits [6, 7, 8, 9, 10, 11, 12, 13, 14,
15, 16]. The very essence of qubits as controllable quantum systems dictates
both their non-equilibrium nature as well as their coupling to extensive
number of external degrees of freedom.

An economic and, in many cases, justified way of treating such driven
dissipative quantum systems is to employ the Markovian (i.e. time-local)
approximation. Under such assumption dynamics of a reduced density matrix is
given by a Lindblad equation [17, 18]. Historically, the study of Lindbladian
dynamics was primarily restricted to systems with a few degrees of freedom,
with most of the focus coming from the quantum optics literature [19, 20, 21,
22]. For example, Lindbladian evolution of a single two-level system is fully
equivalent to the set Bloch equations. Other examples include the dynamics of
parametrically driven oscillators [23, 24, 25], cavity QED and coupled cold
atom-cavity systems [26, 27], etc. Their considerations lead to a number of
insightful physical results and powerful theoretical approaches.

The modern quantum computation platforms, such as Josephson or ionic traps,
fall squarely into the realm of many-body systems. In addition to these, we
also mention driven-dissipative quantum fluids and Bose condensates [28, 29,
30, 31, 32, 33, 34, 35, 36], dynamics of large networks of coupled parametric
oscillators and optical cavities [37, 38, 39, 40, 41, 42, 43], and monitored
dynamics of spatially extended systems [44, 45, 46, 47], all of which
implicitly involve extensively large numbers of degrees of freedom. Indeed,
already $N=50$ connected qubit devices gives rise to the Hilbert space
dimension ${\cal N}=2^{50}$, which is well beyond traditional single-particle
matrix manipulation techniques. This calls for the developing of many-body
field theoretical techniques geared towards description of the Lindbladian (as
opposed to von Neumann) evolution.

An important step in this direction is consideration of many-body bosonic or
fermionic systems, traditionally described in the occupation number basis via
the algebra of creation/annihilation operators. In this approach both many-
body effective Hamiltonian and a set of many-body quantum jump operators are
all expressed as polynomials of such creation/annihilation operators. It is
important to remember, though, that despite of a deceptively simple appearance
all these operators act in the exponentially large (in the number of the
degrees of freedom) Hilbert space. Therefore a brute force numerical solution
of the corresponding Lindblad equation requires diagonalization of ${\cal
N}^{2}\times{\cal N}^{2}$ matrix (for, e.g., fermionic case). Clearly this is
not a productive direction.

Various techniques have been developed for studying dynamics of quadratic
Lindbladians, i.e. those with the Hamiltonian given by a quadratic form, while
all quantum jump operators by linear forms of the creation/annihilation
operators. One such approach is provided by the so-called “third quantization”
technique [48, 49, 50, 51, 52], based on the use of algebras of bosonic or
fermionic superoperators. This approach shows that ${\cal N}^{2}$ eigenvalues
of a quadratic many-body Lindbladian may be constructed from $N$ complex
eigenvalues of a certain $N\times N$ non-Hermitian matrix, using conventional
bosonic or fermionic occupation numbers.

Let us also mention the notable topological classification of Lindbladian
fermions [53, 54, 55] and various results pertaining to both bosonic and
fermionic Gaussian states [56, 57, 58, 59, 60, 61, 62, 63]. These techniques
have been variously applied to the study numerous problems, including the
Bogoliubov spectrum of driven-dissipative condensates [29], exact solutions of
nonlinear integrable systems [64], and the study topological properties of
various low-dimensional dissipative systems [65, 66, 67, 68, 69].

The purpose of this manuscript is to review and further develop an alternative
apparatus, based on coherent state functional integral field-theoretical
treatment of the Lindbladian dynamics, pionered by Sieberer, Buchhold, and
Diehl [3]. This technique originates from the Keldysh theory [70] of the
underlying Von Neumann dynamics of the interacting system-bath pair. Upon
integrating out the bath degrees of freedom and adopting Markovian
approximation for the bath-induced self-energy, one ends up with a time-local
effective action. The latter is fully equivalent to the many-body operator
Lindblad equation [3]. It constitutes, however, a more convenient starting
point for calculation of various observables, correlation functions, linear
response characteristics, collective modes, etc. It is also indispensable for
generalizations beyond the quadratic theory.

In the quadratic approximation this approach naturally reproduces the results
of the third quantization for the spectra of the Lindbladian superoperators.
Its advantage is in making unmistakably clear that this information is only
part of the whole picture. While in equilibrium systems, both statistical
weights and the dynamics are determined by the same set of energies, this is
not the case for driven-dissipative Lindbladian dynamics. In this case the
complex spectrum of the dynamical relaxation is not directly related to (real)
statistical weights of the stationary (but non-equilibrium) density matrix.
The latter is determined by the stationary distribution function
$\check{F}_{\mathrm{st}}$, which naturally emerges as one of the main building
blocks of the Keldysh treatment.

One of the main goals of this text is to draw distinctions between the
transient relaxation spectra (derived from the eigenvalues of a certain non-
Hermitian $N\times N$ matrix $\check{H}$) and a Hermitian $N\times N$
stationary distribution $\check{F}_{\mathrm{st}}$. Already for quadratic
Lindbladians finding $\check{F}_{\mathrm{st}}$ requires solving a linear
kinetic equation. The latter takes the form of the so-called continuous-time
Lyapunov matrix equation, well-known in the dynamical systems literature in
the context of stability and control of linear systems [71]. We show that, on
the one hand, $\check{F}_{\mathrm{st}}$ determines a host of observables and
correlation function of physics interest. On the other hand, its properties
are often qualitatively different from those of $\check{H}$. While the latter
are frequently emphasized in the literature, the former are undeservedly
overlooked. For example, certain non-analyticities (associated with the
exceptional points) in the $\check{H}$ spectra, do not show up in the
$\check{F}_{\mathrm{st}}$ spectra. Similarly, while the band structure of
$\check{H}$ often exhibits interesting topological characteristics, the band
structure of the corresponding $\check{F}_{\mathrm{st}}$ may be topologically
trivial.

The structure of the manuscript is rather straightforward. In section 2 we
develop major aspects of the field-theoretical treatment of the Lindbladian
dynamics for bosonic and fermionic many-body systems. Here we emphasize the
role of the stationary distribution and explain the origin of the Lyapunov
equation. We also derive generic expressions for observables and linear
response and consider exceptional points. Section 3 is devoted to a number of
pedagogic examples illustrating various aspects of many-body Lindbladian
dynamics. Some technicalities are delegated to appendices.

##  2 Formalism

This section discusses the general formalism of quadratic theories of both
bosons and fermions. It begins with an introduction to the Keldysh formalism
as it relates to Lindbladians. With this established, one can obtain the many-
body Lindbladian spectrum and stationary density matrix using the Keldysh
Green’s function formalism.

###  2.1 Lindblad and Keldysh

In the operator formalism approach to open quantum systems, one studies the
reduced density matrix, $\rho$, of the system of interest, resulting from
tracing out the environment Hilbert space. The tracing out of environmental
degrees of freedom coupled to the system generates additional terms in the
evolution equation for $\rho$ alongside the standard von-Neumann part,
resulting in an effective non-equilibrium dynamics. In situations where memory
effects may be neglected, time evolution of $\rho$ is described by the
Lindblad master equation [19],

$\partial_{t}\rho=\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\mathcal{L}}}}\rho,$
(1a)
$\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\mathcal{L}}}}=-i[\hat{\mathcal{H}},\cdot]+\sum_{v}\Big{(}\hat{\mathcal{L}}_{v}\cdot\hat{\mathcal{L}}_{v}^{\dagger}-\frac{1}{2}\\{\hat{\mathcal{L}}_{v}^{\dagger}\hat{\mathcal{L}}_{v},\cdot\\}\Big{)},$
(1b)

where the latter equation defines the Lindbladian superoperator. The Hermitian
operator $\hat{\mathcal{H}}$ is an effective (possibly renormalized by the
environment) Hamiltonian of the system. The jump operators
$\hat{\mathcal{L}}_{v}$ (in general non-Hermitian) specify channels through
which the system is coupled to its environment.

The Lindbladian plays a role analogous to the Hamiltonian in closed quantum
systems in that determining its eigenvectors and eigenvalues provides complete
knowledge of the system’s dynamics. One thus seeks to solve the superoperator
eigenvalue problem
$\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\mathcal{L}}}}\rho_{\Lambda}=\Lambda\rho_{\Lambda}$.
A given density matrix $\rho$ will be a superposition of the Lindbladian
eigenvectors $\rho_{\Lambda}$. The corresponding Lindbladian eigenvalues will
generically be complex, coming in complex-conjugate pairs, with the real and
imaginary parts corresponding to the rates of decay and coherent (phase)
rotation of the $\rho_{\Lambda}$ component of $\rho$.

Dynamical stability at long times requires all $\Lambda$ to have non-positive
real parts. This requirement is comparable to a Hamiltonian spectrum being
bounded from below in dynamics of a closed quantum system. The
$\rho_{\Lambda}$ with purely imaginary $\Lambda$ play a special role: they do
not dissipate. There is always at least one zero eigenvalue, corresponding to
a stationary state $\rho_{\mathrm{st}}$. In general there may be multiple
stationary states spanning a multidimensional operator subspace. The structure
and dimension of the space of stationary states is determined by the
symmetries of the system [72, 73, 74, 75, 76, 77]. For a generic Lindbladian
without additional symmetry however, the stationary state is unique.

The focus of this manuscript is on many-body Lindbladians, in which the
Hilbert space of states is either a bosonic or fermionic Fock space. It is
thus convenient to use the creation/annihilation operator basis $\hat{a}_{j}$
and $\hat{a}_{j}^{\dagger}$, where $j\leq N$ is a generic internal index
encompassing e.g. spin, flavor, orbital number, space or momentum, etc. In
this basis, $\hat{\mathcal{H}}=\mathcal{H}(\hat{a}_{j}^{\dagger},\hat{a}_{j})$
and $\hat{\mathcal{L}}_{v}=\mathcal{L}_{v}(\hat{a}_{j}^{\dagger},\hat{a}_{j})$
where $\mathcal{H}$ and $\mathcal{L}_{v}$ without hats are polynomial
functions. The focus below is on $\mathcal{H}$ quadratic and $\mathcal{L}_{v}$
linear. In such theories, the Lindbladian dynamics is akin to non-interacting
Hamiltonian dynamics: it is possible to solve exactly the full many-body
problem from studying single-particle quantities. In particular, a generic
many-body eigenvalue takes the form:

$\Lambda_{n_{1}...n_{2N}}=-i\sum_{s}n_{s}\epsilon_{s},$ (2)

where the $2N$ quantum numbers $n_{s}$ are integer-valued occupation numbers
and $\epsilon_{s}$ are the solutions to a single-particle non-Hermitian
eigenvalue problem. The stationary state density matrix can be obtained from
stationary solution to the single-particle quantum kinetic equation.

The formal machinery used to extract this information is the Keldysh path
integral and the corresponding theory of non-equilibrium Green’s functions
[70]. Using the formalism of [3], the dynamics encoded in the Lindblad
equation can be mapped to a coherent state path integral. In broad strokes,
this is achieved by expressing the density matrix at time $t$ in terms its
value at an initial time $t_{0}$ via a time evolution superoperator
$\rho(t)=\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\mathcal{U}}}}_{t}\rho(t_{0})$,
defined by
$\exp(t\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\mathcal{L}}}})$.
One may then introduce the Keldysh partition function as the superoperator
analogue of the propagator,

$Z=\operatorname{tr}\Big{(}\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\mathcal{U}}}}_{t}\rho(t_{0})\Big{)}.$
(3)

which is always identically equal to 1 due to the density matrix
normalization. $Z$ can be brought into the form a path integral by cutting the
time interval into infinitesimal slices via the Trotter formula, so that one
may write
$\exp(\delta_{t}\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\mathcal{L}}}})\simeq
1+\delta_{t}\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\mathcal{L}}}}$
on each time slice where $\delta_{t}$ is the duration of a slice. By inserting
factors of a coherent state resolution of identity in-between each slice on
both sides of the density matrix, operators are converted into fields. For a
single bosonic mode, one uses the bosonic coherent states defined by
$\hat{a}\ket{\phi}=\phi\ket{\phi}$ where $\phi$ is a complex number
(generalization to $N$ bosons is automatic, for details on fermionic Keldysh
integrals, see Section 2.6). Each coherent state component
$\ket{\phi^{+}_{1}}\bra{\phi^{-}_{1}}$ of the density matrix at time $t$,
being acted upon by the superoperator
$\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\mathcal{L}}}}$, leads to
the following matrix elements:

$\bra{\phi^{+}_{2}}\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\mathcal{L}}}}(\ket{\phi^{+}_{1}}\bra{\phi^{-}_{1}})\ket{\phi^{-}_{2}}=-ie^{\bar{\phi}^{+}_{2}\phi^{+}_{1}+\bar{\phi}^{-}_{1}\phi^{-}_{2}}\mathcal{K}(\bar{\phi}^{+}_{2},\phi^{+}_{1},\bar{\phi}^{-}_{1},\phi^{-}_{2}),$
(4)

where the Keldysh “Hamiltonian” $\mathcal{K}$ has the same form as the
Lindbladian upon replacing creation operators $\hat{a}$, that multiply the
density matrix on the left/right, with $\phi^{\pm}$ respectively. Provided the
functional form of the Hamiltonian $\mathcal{H}$ and jump operators
$\mathcal{L}_{v}$ are normal ordered, one may write:

$\mathcal{K}(\bar{\phi}^{+},\phi^{+},\bar{\phi}^{-},\phi^{-})=\mathcal{H}^{+}-\mathcal{H}^{-}+i\mathcal{D}(\bar{\phi}^{+},\phi^{+},\bar{\phi}^{-},\phi^{-}),$
(5a)
$\mathcal{D}(\bar{\phi}^{+},\phi^{+},\bar{\phi}^{-},\phi^{-})=\sum_{v}\Big{(}\bar{\mathcal{L}}_{v}^{-}\mathcal{L}_{v}^{+}-\frac{1}{2}\bar{\mathcal{L}}_{v}^{+}\mathcal{L}_{v}^{+}-\frac{1}{2}\bar{\mathcal{L}}_{v}^{-}\mathcal{L}_{v}^{-}\Big{)},$
(5b)

where $\mathcal{H}^{\pm}=\mathcal{H}(\bar{\phi}^{\pm},\phi^{\pm})$ and
$\mathcal{L}_{v}^{\pm}=\mathcal{L}_{v}(\bar{\phi}^{\pm},\phi^{\pm})$. Note
that the second equality holds for quadratic theories considered here because
the normal ordering of the product
$\hat{\mathcal{L}}_{v}^{\dagger}\hat{\mathcal{L}}_{v}$ is equivalent to the
product of the normal ordering up to a trivial constant.

Upon re-exponentiation in the limit $\delta_{t}\to 0$, this retrieves a
functional integral in terms of two sets of fields
$\phi^{\pm}(t),\bar{\phi}^{\pm}(t)$ corresponding to multiplication of the
density matrix by the creation/annihilation operators
$\hat{a},\hat{a}^{\dagger}$ on the left or right side in the Lindbladian. It
is conventional to present this functional integral using the Keldysh rotated
basis of “classical” and “quantum” fields,
$\phi^{\mathrm{c},\mathrm{q}}=(\phi^{+}\pm\phi^{-})/\sqrt{2}$. All together,
one has for the partition function [3],

$Z=\int\text{D}\bar{\phi}^{\alpha}\text{D}\phi^{\alpha}e^{iS[\bar{\phi}^{\alpha},\phi^{\alpha}]},$
(6)

with $\alpha=\mathrm{c},\mathrm{q}$. The Keldysh action $S$ is the given by
the time integral of a Lagrangian defined by the Keldysh Hamiltonian
$\mathcal{K}$ defined as a function of $\phi^{\pm}$ by eq. (5a),

$S=\int\text{d}t\Big{(}\bar{\phi}^{\mathrm{q}}i\partial_{t}\phi^{\mathrm{c}}+\bar{\phi}^{\mathrm{c}}i\partial_{t}\phi^{\mathrm{q}}-\mathcal{K}(\bar{\phi}^{\alpha},\phi^{\alpha})\Big{)}.$
(7)

The value of density matrix $\rho(t_{0})$ at the initial time is contained in
the boundary conditions of the path integral.

With the Keldysh path integral in hand, $n$-point correlation functions can be
calculated via operator insertion at different times on either side of the
density matrix. By extension, the expectation of a quadratic observable
$\hat{\mathcal{O}}=\hat{\mathcal{A}}^{\dagger}\check{O}\hat{\mathcal{A}}$,
where $\hat{\mathcal{A}}=[\hat{a}\ \hat{a}^{\dagger}]$ is the Nambu space
spinor and $\check{O}$ is a matrix, at a finite time can be computed as the
expectation inside the function integral,

$\braket{\hat{\mathcal{O}}(t)}=\int\text{D}\bar{\phi}^{\alpha}\text{D}\phi^{\alpha}e^{iS}\mathcal{O}\big{(}\bar{\phi}^{-}(t),\phi^{+}(t)\big{)},$
(8)

where $\mathcal{O}(\bar{\phi},\phi)=\bar{\Phi}\check{O}\Phi$ is the classical
function of operators replaced with fields, with $\Phi=[\phi\ \bar{\phi}]$ the
Nambu space vector. The integral can performed by expressing $\mathcal{O}$ as
a function of the Keldysh fields. As an example which will be relevant below,
one may consider the single-particle bosonic covariance matrix
$\langle\\{\hat{\mathcal{A}},\hat{\mathcal{A}}^{\dagger}\\}\rangle$. This
Nambu space matrix-valued expectation reduces to the equal time two-point
expectation of classical fields
$\langle\\{\hat{\mathcal{A}},\hat{\mathcal{A}}^{\dagger}\\}\rangle(t)=\langle\Phi^{\mathrm{c}}(t)\bar{\Phi}^{\mathrm{c}}(t)\rangle$.

One of the main advantages of the Keldysh formalism is that this procedure is
straight-forwardly generalized to expectations of fields with different time
arguments. For a quadratic theory, the two-point functions are the most
important as they can be used to compute all higher-order $n$-point functions
via Wick’s theorem. Combining all four fields together into a single Keldysh-
Nambu vector $\Phi=[\phi^{\mathrm{c}}\ \bar{\phi}^{\mathrm{c}}\
\phi^{\mathrm{q}}\ \bar{\phi}^{\mathrm{q}}]$, the two-point functions define
the spectral and Keldysh Green’s functions,

$\begin{bmatrix}i\check{G}^{\mathrm{K}}(t,t^{\prime})&i\check{G}^{\mathrm{R}}(t,t^{\prime})\\\
i\check{G}^{\mathrm{A}}(t,t^{\prime})&0\end{bmatrix}=\braket{\Phi(t)\bar{\Phi}(t^{\prime})}.$
(9)

Note that Green’s functions defined using this convention are matrices on the
Nambu space so as to account for the possibility of non-zero anomalous
expectations e.g. $\langle\phi^{\alpha}(t)\phi^{\beta}(t^{\prime})\rangle$.
The Keldysh Green’s function $\check{G}^{\mathrm{K}}$ contains information
about the distribution function of the system; at equal times one can see that
it is just the covariance matrix from the example above,
$\check{G}^{\mathrm{K}}(t,t)=\langle\\{\hat{\mathcal{A}},\hat{\mathcal{A}}^{\dagger}\\}(t)\rangle$.
It is conventional to represent the Green’s functions diagrammatically, as
shown in fig. 1. This relation generalizes to $N$ particles and, as discussed
below, can be used to compute the stationary state density matrix. The
spectral Green’s functions $\check{G}^{\mathrm{R},\mathrm{A}}$ contain purely
dynamical information and are independent of the state of the system. As
discussed in the following section, they are related to the single-particle
eigenvalue spectrum.

Figure 1: Diagrammatic representations of the three Green’s functions from the
Keldysh theory. The solid lines denotes the classical field
$\phi^{\mathrm{c}}$ and the dashed line, the quantum field
$\phi^{\mathrm{q}}$.

###  2.2 Bosons

A generic quadratic Lindbladian for a system of $N$ Bosons is defined by a
quadratic Hamiltonian and a set of linear Jump operators,

$\hat{\mathcal{H}}=\sum_{i,j}^{N}\big{(}\Delta_{ij}\hat{a}^{\dagger}_{i}\hat{a}_{j}+\lambda_{ij}\hat{a}_{i}\hat{a}_{j}+\lambda^{*}_{ij}\hat{a}^{\dagger}_{i}\hat{a}^{\dagger}_{j}\big{)},$
(10a)
$\hat{\mathcal{L}}_{v}=\sum_{j}^{N}(\mu_{vj}\hat{a}_{j}+\nu_{vj}\hat{a}^{\dagger}_{j}),$
(10b)

where the hat $\hat{a}_{j}$’s are bosonic creation operators,
$[\hat{a}_{j},\hat{a}^{\dagger}_{i}]=\delta_{ij}$. One can in principle allow
the parameter matrices to vary as functions time, but for simplicity here only
time-independent models are considered. The corresponding Keldysh action of
eq. (7) can be arranged into a quadratic form:

$S=\frac{1}{2}\int\text{d}t\
\bar{\Phi}\begin{bmatrix}0&\check{\tau}^{3}i\partial_{t}-\check{H}_{0}-i\check{Q}\\\
\check{\tau}^{3}i\partial_{t}-\check{H}_{0}+i\check{Q}&i\check{D}\end{bmatrix}\Phi,$
(11)

where $\check{\tau}^{i}$ are the Pauli matrices acting in Nambu space and the
fields $\phi^{\alpha}$ are understood as vectors with $N$ entries
$\phi^{\alpha}_{j}$ so that the Keldysh-Nambu spinor $\Phi$ has a total of
$4N$ entries.

The operators $\check{H}_{0}$, $\check{Q}$, and $\check{D}$ are Hermitian
$2N\times 2N$ Nambu space matrices:

$\check{H}_{0}=\begin{bmatrix}\Delta&2\lambda^{\dagger}\\\
2\lambda&\Delta^{\mathrm{T}}\end{bmatrix},$ (12a)
$\check{Q}=\frac{1}{2}\begin{bmatrix}\gamma-\tilde{\gamma}&\eta_{\mathrm{a}}^{\dagger}\\\
\eta_{\mathrm{a}}&\tilde{\gamma}^{\mathrm{T}}-\gamma^{\mathrm{T}}\end{bmatrix},$
(12b)
$\check{D}=\begin{bmatrix}\gamma+\tilde{\gamma}&\eta_{\mathrm{s}}^{\dagger}\\\
\eta_{\mathrm{s}}&\gamma^{\mathrm{T}}+\tilde{\gamma}^{\mathrm{T}}\end{bmatrix},$
(12c)

where

$\gamma_{ij}=\sum_{v}\mu^{*}_{vi}\mu_{vj},\quad\tilde{\gamma}_{ij}=\sum_{v}\nu_{vi}\nu^{*}_{vj},\quad\eta_{ij}=\sum_{v}\nu^{*}_{vi}\mu_{vj},$
(13)

and $\eta_{\mathrm{s},\mathrm{a}}=\eta\pm\eta^{\mathrm{T}}$. The $N\times N$
parameter matrices have index symmetries:

$\Delta=\Delta^{\dagger},\quad\lambda=\lambda^{\mathrm{T}},\quad\gamma=\gamma^{\dagger},\quad\tilde{\gamma}=\tilde{\gamma}^{\dagger},\quad\eta_{\mathrm{s},\mathrm{a}}=\pm\eta_{\mathrm{s},\mathrm{a}}^{\mathrm{T}}.$
(14)

Note that the matrix $\check{H}_{0}$ is nothing more than the single-particle
Hamiltonian,
$\hat{\mathcal{H}}=\frac{1}{2}\hat{\mathcal{A}}^{\dagger}\check{H}_{0}\hat{\mathcal{A}}$.
The matrices $\check{Q}$ and $\check{D}$ are determined entirely by the
coupling to the environment through the jump operators.

The dynamic matrix $\check{H}=\check{\tau}^{3}(\check{H}_{0}-i\check{Q})$ is a
non-Hermitian matrix that replaces the single-particle Hamiltonian in coherent
many-body systems. Its $2N$ eigenvalues are the obtained by solving the non-
Hermitian eigenvalue problem,

$\check{H}\ket{s}=\epsilon_{s}\ket{s},$ (15)

where $\epsilon_{s}$ are the eigenvalues of $\check{H}$ obtained through
diagonalization by a generically non-Unitary matrix $\check{U}$,

$\big{(}\check{U}\check{H}\check{U}^{-1}\big{)}_{ss^{\prime}}=\epsilon_{s}\delta_{ss^{\prime}}.$
(16)

Note that the complex eigenvalues of $\check{H}$ come in complex-conjugate
pairs on due to the Nambu space particle-hole symmetry. The $\epsilon_{s}$ act
as eigenvalues of the single-particle sector of the Lindbladian. In the
standard quantum theory, the absence of interactions implies the energies of
the full many-body system to be determined by filling the single particle
states with various numbers of particles. This is essentially true in the
Lindbladian theory as well, except that the single-particle states $\ket{s}$
have a finite lifetime on account of the complex single-particle “energies”
$\epsilon_{s}$ having migrated into the bottom half of the complex plane, see
Fig. (2). The expression in eq. (2) for a many-body eigenvalue is just the
assigning of bosonic occupation numbers to this single-particle sector.
Demonstrating this fact explicitly can be achieved either by semiclassical
quantization of the Keldysh action (see Appendix A) or through the third
quantization formalism [48] (see Appendix B for connections between the
Keldysh and third quantization formalisms).

To gain further intuition about the nature of the dynamics, consider the the
classical mechanics of the Keldysh action eq. (11). The equation of motion of
the classical field $\Phi^{\mathrm{c}}=[\phi^{\mathrm{c}}\
\bar{\phi}^{\mathrm{c}}]$ with the quantum field set to zero
$\Phi^{\mathrm{q}}=[\phi^{\mathrm{q}}\ \bar{\phi}^{\mathrm{q}}]=0$ is

$i\partial_{t}\Phi^{\mathrm{c}}=\check{H}\Phi^{\mathrm{c}}.$ (17)

This is a non-Hermitian Schrödinger equation where the classical field acts as
the single-particle wave function. This can equivalently be conceptualized in
first-quantized language as an equation of motion for the coordinate on the
$N$-particle phase space. Due to the non-Hermiticity of the dynamic matrix,
the classical mechanics this equation encodes is dissipative: the phase
portrait will consist of spiraling paths centered at the origin. The dynamics
are only stable when all of the eigenvalues of $\check{H}$ have non-positive
imaginary part, so that all phase space trajectories fall into the origin
rather than running to infinity. This behaviour is not ensured for generic
choices of the parameter matrices and can fail if the magnitudes of $\lambda$
or $\tilde{\gamma}$ are large compared to other parameters. These situations
are unphysical, being associated with either an unstable Hamiltonian in which
the potential of one or more coordinates in the phase space is inverted or
with a situation where the rate of particle gain is greater than loss,
resulting in an uncontrolled pumping of quanta into the system. At the
threshold of such an instability the eigenvalues of $\check{H}$ can be purely
imaginary, resulting in a coherent orbiting around the origin. This
corresponds to a closing of the dissipative gap in the Lindbladian spectrum
and stable long-time dynamics beyond a single stationary state.

Figure 2: Example spectra of $\check{H}$ plotted in the complex plane of the
eigenvalues $\epsilon_{s}$. The eigenvalues labeled in red correspond to
eigenvectors that evolve in a purely dissipative way. The purple eigenvalues
correspond to modes with simultaneous dissipation and coherent rotation. The
black eigenvalue at the origin denotes the stationary state(s), which can be
compared to energy eigenstates in a closed quantum system. The blue
eigenvalues aligned along the real axis correspond to limiting cycles which do
not decay at long times; these are analogous to coherent superpositions of
different energy eigenstates.

The spectral Green’s functions
$\check{G}^{\mathrm{R},\mathrm{A}}(t,t^{\prime})$ contain the same dynamical
information in their pole structure. They can be read off as the off-diagonal
blocks of the inverse of the quadratic form in eq. (11). This is equivalent to
inverting the differential operator in eq. (17). The spectral Green’s
functions are independent of the distribution of the system and are thus
always functions of the difference of their time arguments,

$\check{G}^{\mathrm{R}}(t,t^{\prime})=-i\theta(t-t^{\prime})e^{-i(t-t^{\prime})\check{H}}\check{\tau}^{3},\quad\quad\check{G}^{\mathrm{A}}(t,t^{\prime})=i\theta(t^{\prime}-t)\check{\tau}^{3}e^{-i(t-t^{\prime})\check{H}^{\dagger}}.$
(18)

Upon Fourier transform with respect to the difference of the time arguments
$t-t^{\prime}$, they are adopt a simple form of the resolvent of the dynamic
matrix,

$\check{G}^{\mathrm{R}}(\epsilon)=\frac{1}{\epsilon-\check{H}}\check{\tau}^{3},\quad\quad\check{G}^{\mathrm{A}}(\epsilon)=\check{\tau}^{3}\frac{1}{\epsilon-\check{H}^{\dagger}}.$
(19)

The poles of the Green’s functions are located at the eigenvalues
$\epsilon_{s}$. This can be compared to the association between the energies
of single-particle states and poles in standard quantum theory.

###  2.3 Lyapunov Equation

In this section, the stationary state of the Lindbladian is discussed.
Contrasting to the spectral Green’s functions discussed above, the Keldysh
Green’s function $\check{G}^{\mathrm{K}}$ depends on the distribution of the
system. Conventionally, one parameterizes the Keldysh Green’s function in
terms of the spectral Green’s functions and a Hermitian matrix
$\check{F}(t,t^{\prime})$,

$\check{G}^{\mathrm{K}}=\check{G}^{\mathrm{R}}\circ\check{\tau}^{3}\check{F}-\check{F}\check{\tau}^{3}\circ\check{G}^{\mathrm{A}},$
(20)

where the composition $\circ$ denotes matrix composition both in the time
argument and the Nambu space. Note the additional factor of $\check{\tau}^{3}$
here compared to the standard convention in [70] is a consequence of the
symplectic structure of the bosonic Nambu space. The matrix $\check{F}$ acts
as a single-particle distribution matrix. Acting on this equation on the left
by $\check{\tau}^{3}\check{G}^{\mathrm{R}-1}$ and on the right by
$\check{G}^{\mathrm{A}-1}\check{\tau}^{3}$ retrieves the quantum kinetic
equation for $\check{F}$,

$[\partial_{t}\overset{\circ}{,}\check{F}]=-i\big{(}\check{H}\check{F}-\check{F}\check{H}^{\dagger}\big{)}+\check{\tau}^{3}\check{D}\check{\tau}^{3},$
(21)

where the $\partial_{t}$ and $\check{D}$ are understood to be diagonal
functions of their two time arguments, i.e. coming with factors of
$\delta(t-t^{\prime})$. Note the use of the relation
$i\check{D}=-\check{G}^{\mathrm{R}-1}\circ\check{G}^{\mathrm{K}}\circ\check{G}^{\mathrm{A}-1}$
obtained by inverting the quadratic form in the action eq. (11) to derive this
equation. This is valid because the path integral in eq. (6) is Gaussian on
the bulk of the time contour, except possibly at the initial time $t_{0}$ in
the case of a non-Gaussian initial density. The initial density lives on the
boundary of the time contour and determines the boundary condition of the
quantum kinetic equation.

In a stationary state, $\check{G}^{\mathrm{K}}$ and by extension $\check{F}$
are independent of the the central time $t+t^{\prime}$. This nullifies the
left-hand side of eq. (21), meaning the right-hand side must be independently
nullified by a stationary solution. This is achieved by a time-diagonal
ansatz, $\check{F}(t-t^{\prime})=\check{F}_{\mathrm{st}}\delta(t-t^{\prime})$,
where $\check{F}_{\mathrm{st}}$ is a time-independent matrix in the orbital
and Nambu spaces obeying the relation:

$0=-i\big{(}\check{H}\check{F}_{\mathrm{st}}-\check{F}_{\mathrm{st}}\check{H}^{\dagger}\big{)}+\check{\tau}^{3}\check{D}\check{\tau}^{3}.$
(22)

This is a complex Lyapunov equation. There is a unique solution provided that
all of the eigenvalues of $\check{H}$ have finite imaginary parts [71]. This
implies that the Lindblad equation possesses a unique stationary state that
arbitrary initial conditions converge towards at long times. The existence of
additional stationary states for bosonic quadratic Lindbladians thus occurs
only at the brink of dynamical instability, when a parameter is tuned so that
the dissipative gap closes.

Provided the stationary state is unique, one can solve eq. (22) in the
eigenbasis of $\check{H}$,

$\big{(}\check{U}\check{F}_{\mathrm{st}}\check{U}^{\dagger}\big{)}_{ss^{\prime}}=\frac{-i}{\epsilon_{s}-\epsilon^{*}_{s^{\prime}}}\big{(}\check{U}\check{\tau}^{3}\check{D}\check{\tau}^{3}\check{U}^{\dagger}\big{)}_{ss^{\prime}}.$
(23)

Note that $\check{D}$ is generically not diagonal in this basis, meaning that
off-diagonal elements of $\check{F}_{\mathrm{st}}$ are finite. This can be
compared to the equilibrium theory, in which the preferred stationary
$\check{F}$ is the thermal distribution, which is diagonal in Nambu space in
the eigenbasis of the single-particle Hamiltonian. The relation in eq. (20) is
equivalent to the Fluctuation-Dissipation theorem for each particle species.
In the Lindbladian setting, $\check{F}_{\mathrm{st}}$ is
$\epsilon$-independent and generically develops off-diagonal elements in the
eigenbasis of $\check{H}$ and so is more naturally thought of as a matrix. An
alternative but equivalent interpretation of $\check{F}$ is given by
integrating eq. (20), giving the relation
$i\check{G}^{\mathrm{K}}(t,t)=\check{F}_{\mathrm{st}}$. That is,
$\check{F}_{\mathrm{st}}$ is equivalent to the covariance matrix
$\braket{\\{\hat{\mathcal{A}},\hat{\mathcal{A}}^{\dagger}\\}}$ discussed in
section 2.1.

As an alternative to eq. (23), one can instead express
$\check{F}_{\mathrm{st}}$ in its eigenbasis. Letting $\check{U}_{F}$ be a
diagonalizing transformation of $\check{\tau}^{3}\check{F}_{\mathrm{st}}$, one
has [78]:

$\check{U}_{F}\check{\tau}^{3}\check{F}_{\mathrm{st}}\check{U}_{F}^{\dagger}=\mathrm{diag}\big{(}\coth(\beta_{1}/2),...,-\coth(\beta_{1}/2),...\big{)},$
(24)

where the $N$ numbers $\beta_{j}$ parametrize the eigenvalues of
$\check{F}_{\mathrm{st}}$ and act as effective inverse temperatures for the
$j$th eigenvector. For a dynamically stable theory, $0<\beta_{j}\leq\infty$.
As mentioned above, this is in general not the same basis in which the dynamic
matrix $\check{H}$ is diagonal, $\check{U}\neq\check{U}_{F}$. This is in stark
contrast to the equilibrium theory of quadratic Hamiltonians, in which the
thermal state is the Gaussian state with
$\hat{\mathcal{H}}_{\mathrm{st}}=\hat{\mathcal{H}}$. In equilibrium, the bases
in which the dynamics and the distribution are diagonal are the same. Out of
equilibrium, as is the case for the Lindbladian theory, this is generically
untrue.

The form of the stationary density matrix $\rho_{\mathrm{st}}$ can be obtained
using the identity of $\check{F}_{\mathrm{st}}$ as the covariance matrix. The
covariance matrix is a central object in the theory of Gaussian states and is
known to be equivalent to full knowledge of such a state [56, 57, 58]. A
Gaussian state is a state with a density matrix given by the exponentiation of
some quadratic operator of the form of eq. (10a). For a quadratic Lindbladian
with a unique stationary state, the state will be Gaussian. As such, one can
write the stationary density matrix $\rho_{\mathrm{st}}$ in terms of an
effective Hermitian Hamiltonian,

$\rho_{\mathrm{st}}\propto\exp(-\hat{\mathcal{H}}_{\mathrm{st}}),$ (25a)
$\hat{\mathcal{H}}_{\mathrm{st}}=\frac{1}{2}\hat{\mathcal{A}}^{\dagger}\check{H}_{\mathrm{st}}\hat{\mathcal{A}},$
(25b)

where the proportionality is determined by the normalization
$\operatorname{tr}(\rho)=1$ and $\check{H}_{\mathrm{st}}$ is generically not
the same as $\check{H}_{0}$. The effective Hamiltonian can be found from the
stationary distribution matrix through the relation [59]:

$\check{F}_{\mathrm{st}}\check{\tau}^{3}=\coth(\check{\tau}^{3}\check{H}_{\mathrm{st}}/2).$
(26)

In the eigenbasis of $\check{F}_{\mathrm{st}}$, it adopts a particularly
simple form,

$\hat{\mathcal{H}}_{\mathrm{st}}=\sum_{j}^{N}\beta_{j}\hat{b}^{\dagger}_{j}\hat{b}_{j},$
(27)

where the diagonal basis bosons are defined by $[\hat{b}\
\hat{b}^{\dagger}]=\check{U}_{F}[\hat{a}\ \hat{a}^{\dagger}]$. The $\beta_{j}$
determine the average populations of the $\hat{b}$ bosons in the stationary
state. Note that in situations where there is not a unique stationary state,
some stationary states may be non-Gaussian and the value of the Keldysh
Green’s function at long times depends on the initial conditions.

###  2.4 Observables and Response

With the stationary distribution in hand, one can compute the stationary
expectation of observables. As an example, consider a quadratic observable
$\hat{\mathcal{O}}=\hat{\mathcal{A}}^{\dagger}\check{O}\hat{\mathcal{A}}$. The
expectation at arbitrary finite time after reaching the stationary state is:

$\braket{\hat{\mathcal{O}}}=\frac{1}{2}\operatorname{tr}\Big{(}\check{O}\big{(}\check{F}_{\mathrm{st}}-\check{\tau}^{3}\big{)}\Big{)}.$
(28)

The quantity $\frac{1}{2}(\check{F}_{\mathrm{st}}-\check{\tau}^{3})$ thus acts
like am effective single-particle density matrix. Alternatively, naming the
classical and quantum parts of the observable
$\mathcal{O}^{\mathrm{c},\mathrm{q}}=\frac{1}{2}(\bar{\Phi}^{+}\check{O}\Phi^{+}\pm\bar{\Phi}^{-}\check{O}\Phi^{-})$,
one has

$\braket{\mathcal{O}^{\mathrm{c}}}=\frac{1}{2}\operatorname{tr}\big{(}\check{O}\check{F}_{\mathrm{st}}\big{)}$
(29)

Thus, single-particle traces with the distribution matrix by itself generates
moments of the classical (Weyl-ordered) parts of observables. Correlations of
different observables at different times and observables containing products
of more than two field operators can be obtained using Wick’s theorem.

The response of the system in its stationary state can be studied by
introducing perturbations. It is assumed that the system has been prepared
then allowed to relax to its stationary state. Formally, this amounts to
pushing the initial time into the infinite past $t_{0}\to\infty$ so that the
system retains no memory of its initial condition. Then, at a later finite
time $t_{i}$, some potentially time-dependent perturbations are switched on,
$\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\mathcal{L}}}}\to\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\mathcal{L}}}}+\delta\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\mathcal{L}}}}(t)\theta(t-t_{i})$.
One may of course consider perturbing the system directly by modifying the
Hamiltonian
$\hat{\mathcal{H}}\to\hat{\mathcal{H}}+\delta\hat{\mathcal{H}}(t)$. For a
quadratic perturbation, this is equivalent to changing the single-particle
Hamiltonian by the inclusion of a matrix-valued classical source
$\check{H}_{0}\to\check{H}_{0}+\delta\check{H}_{0}(t)$. In a Lindbladian
problem, one may additionally consider variations to the dissipative part of
the evolution, either by introducing a new jump operator or by varying an
existing jump operator. In both cases, this leads to a modification of the
other two parameter matrices $\check{Q}\to\check{Q}+\delta\check{Q}(t)$ and
$\check{D}\to\check{D}+\delta\check{D}(t)$. In the prior case,
$\delta\check{Q}$ and $\delta\check{D}$ are of the same form as in eq. (12).
In the latter, one may take perturbations to the jump operators by modifying
$\mu\to\mu+\delta\mu(t)$ and $\nu\to\nu+\delta\nu(t)$ in eq. (13) and keeping
only whatever order in $\delta\mu$ and $\delta\nu$ is required.

In the path integral formalism, this is equivalent to introducing a
perturbation to the action $S\to S+\delta S$. This translates to a
perturbation of the Keldysh Hamiltonian in eq. (5a),
$\mathcal{K}\to\mathcal{K}+\delta\mathcal{K}(t)$. The perturbations from
varying $\check{H}_{0}$, $\check{Q}$, and $\check{D}$ respectively are given
by:

$\delta_{H_{0}}\mathcal{K}(t)=\frac{1}{2}\check{\sigma}^{1}_{\alpha\beta}\bar{\Phi}^{\alpha}\delta\check{H}_{0}(t)\Phi^{\beta}$
(30a)
$\delta_{Q}\mathcal{K}(t)=-\frac{1}{2}\check{\sigma}^{2}_{\alpha\beta}\bar{\Phi}^{\alpha}\delta\check{Q}(t)\Phi^{\beta}$
(30b)
$\delta_{D}\mathcal{K}(t)=-\frac{i}{2}\bar{\Phi}^{\mathrm{q}}\delta\check{D}(t)\Phi^{\mathrm{q}}$
(30c)

where $\check{\sigma}^{i}$ denotes the Pauli matrices in the space of Keldysh
indices. For weak perturbations, it suffices to keep only the first order
correction to the measure. This gives the linear response, for which one
obtains for the expectation of an observable $\hat{\mathcal{O}}$ at any finite
time,

$\braket{\hat{\mathcal{O}}}(t)\simeq\braket{\hat{\mathcal{O}}}_{\mathrm{st}}-i\int_{t_{i}}^{t}\text{d}t^{\prime}\braket{\mathcal{O}^{\mathrm{c}}(t)\delta\mathcal{K}(t^{\prime})}_{\mathrm{st}}.$
(31)

This is nothing but the Kubo formula generalized to the Lindbladian context.
The first term in this expression is given by eq. (28). The latter term
includes only the classical part $\mathcal{O}^{\mathrm{c}}$ because only
$\check{F}$ in eq. (28) receives perturbative corrections. It can be computed
using Wick’s theorem, which leads to bubble diagram contributions depicted in
fig. 3. Note that one must be careful to only keep terms corresponding to
fully connected diagrams, see Appendix C.

Figure 3: Bubble diagrams contributing to the linear response of a Lindbladian
perturbation. A Hamiltonian perturbation as in eq. (32) corresponds to the
difference of the retarded and advanced polarization bubbles, as depicted by
the leftmost two diagrams. Dissipative perturbations modifying $\check{Q}$ as
in eq. (33) are comparably given by the of these two bubbles. The rightmost
diagram appears in perturbations modifying $\check{D}$ as in eq. (34). This
diagram in contrast plays no role in the coherent linear response theory.

For a purely Hamiltonian perturbation, one finds a correction in the standard
form of the retarded response function,

$\displaystyle\begin{split}\braket{\mathcal{O}^{\mathrm{c}}(t)\delta_{H_{0}}\mathcal{K}(t^{\prime})}_{\mathrm{st}}&=\braket{\mathcal{O}^{\mathrm{c}}(t)\delta\mathcal{H}^{\mathrm{q}}(t^{\prime})}_{\mathrm{st}}\\\
&=-\frac{1}{2}\operatorname{tr}\Big{(}\check{G}^{\mathrm{R}}(t,t^{\prime})\big{(}\check{\tau}^{3}\check{F}_{\mathrm{st}}\delta\check{H}_{0}(t^{\prime})-\delta\check{H}_{0}(t^{\prime})\check{F}_{\mathrm{st}}\check{\tau}^{3}\big{)}\check{G}^{\mathrm{A}}(t^{\prime},t)\check{O}\Big{)}.\end{split}$
(32)

Perturbations to the dissipative couplings cannot be expressed as expectations
of products of the quantum and classical parts of observables. They do however
admit simple expressions in terms of the distribution matrix,

$\braket{\mathcal{O}^{\mathrm{c}}(t)\delta\mathcal{K}_{Q}(t^{\prime})}_{\mathrm{st}}=-\frac{i}{2}\operatorname{tr}\Big{(}\check{G}^{\mathrm{R}}(t,t^{\prime})\big{(}\check{\tau}^{3}\check{F}_{\mathrm{st}}\delta\check{Q}(t^{\prime})+\delta\check{Q}(t^{\prime})\check{F}_{\mathrm{st}}\check{\tau}^{3}\big{)}\check{G}^{\mathrm{A}}(t^{\prime},t)\check{O}\Big{)},$
(33)
$\braket{\mathcal{O}^{\mathrm{c}}(t)\delta\mathcal{K}_{D}(t^{\prime})}_{\mathrm{st}}=\frac{i}{2}\operatorname{tr}\Big{(}\check{G}^{\mathrm{R}}(t,t^{\prime})\delta\check{D}(t^{\prime})\check{G}^{\mathrm{A}}(t^{\prime},t)\check{O}\Big{)}.$
(34)

Note that both expressions appropriately have a retarded causality, despite
not being expectations of classical and quantum observables like eq. (32).
Also note that the linear response theory for Lindbladians was studied in [79]
using the superoperator formalism. The above formulas are comparable to the
results presented there specialized to many-body bosons.

To go beyond the linear response, rather than keeping higher orders in the
perturbation theory one can instead solve the quantum kinetic equation eq.
(21) to determine the full the non-stationary $\check{F}$. The assumption that
the system reached its stationary state before the perturbation is encoded in
the boundary condition $\check{F}(t,t^{\prime})=\check{F}_{\mathrm{st}}$ for
all $t,t^{\prime}<t_{i}$. Because the perturbation to the Lindbladian is
assumed to be local in time, one can always seek a solution that is time-
diagonal, with $\check{F}(t,t^{\prime})=\delta(t-t^{\prime})\check{F}(t)$.
With this ansatz, the kinetic equation adopts the local form,

$\partial_{t}\check{F}(t)=-i\big{(}\check{H}(t)\check{F}(t)-\check{F}(t)\check{H}^{\dagger}(t)\big{)}+\check{\tau}^{3}\check{D}(t)\check{\tau}^{3}.$
(35)

Assuming one can solve this equation, the expectations of quadratic
observables at finite times can be computed using the appropriate
generalizations of eq.s (28) and (28),

$\braket{\hat{\mathcal{O}}}(t)=\frac{1}{2}\operatorname{tr}\Big{(}\check{O}\big{(}\check{F}(t)-\check{\tau}^{3}\big{)}\Big{)},$
(36a)
$\braket{\mathcal{O}^{\mathrm{c}}(t)}=\frac{1}{2}\operatorname{tr}\big{(}\check{O}\check{F}(t)\big{)}.$
(36b)

The classical parts of observables containing products more than two field
operators can again be obtained using Wick’s theorem. This is valid even with
the non-stationary distribution because the path integral is still a Gaussian
functional integral with the time dependent perturbations; the density matrix
remains a Gaussian state as it evolves in time.

As a check, one can compare the two approaches by combining and rearranging
eq.s (32), (33), and (34). The correction to $\braket{\hat{\mathcal{O}}}$ in
eq. (31) is given by the trace of $\check{O}$ multiplied by the object:

$\int_{t_{i}}^{t}\text{d}t^{\prime}\check{G}^{\mathrm{R}}(t,t^{\prime})\check{\tau}^{3}\Big{(}-i\big{(}\delta\check{H}(t^{\prime})\check{F}_{\mathrm{st}}-\check{F}_{\mathrm{st}}\delta\check{H}^{\dagger}(t^{\prime})\big{)}+\check{\tau}^{3}\delta\check{D}(t^{\prime})\check{\tau}^{3}\Big{)}\check{\tau}^{3}\check{G}^{\mathrm{A}}(t^{\prime},t),$
(37)

which is just the leading-order perturbative solution to eq. (35).

###  2.5 Exceptional Points

This section addresses subtleties that emerge due non-Hermiticity that have
thus far been ignored. The dynamic matrix $\check{H}$, and by extension the
Lindbladian itself, may be non-diagonalizable. This occurs at so-called
exceptional points of the parameter space, at which two or more eigenvalues
merge. This occurs in a fundamentally different way than in standard Hermitian
quantum mechanics, in which the crossing of energy levels is generally avoided
and degeneracies are traditionally understood to be a consequence of some
underlying symmetry. The coalescing of eigenvalues at an exceptional point
should instead marks a bifurcation in the dynamics and is unrelated to
dynamical symmetry.

The prototypical example for how this occurs is the collision of two
eigenvalues on the real axis, see Fig. (4). The eigenvalues of the matrix
$-i\check{H}$, and by extension eigenvalues of the Lindbladian, come in
complex-conjugate pairs. An exceptional point on the real axis thus
corresponds to the spontaneous breaking of this ‘particle-hole symmetry.’ This
signals an under-damped to over-damped bifurcation, in which the corresponding
eigenmodes undergoing damped coherent rotation before the collision experience
pure dissipation after. There is a resonant damping at the exceptional point,
resulting in the transient algebraic gain of one eigenmode. This transient
gain is the generic signature of exceptional points.

Figure 4: Example of tuning a parameter through an exceptional point
separating an under-damped to over-damped dynamical bifurcation. As a
parameter is tuned, a pair of eigenvalues initially located as a conjugate
pair in the left half of the complex plane begin moving toward the real axis
toward one another. After colliding, they will remain stuck on the real axis
but split and move in opposite directions.

To see how this works, suppose the dynamic matrix $\check{H}$ has an
exceptional point where $M$ eigenvalues have collided at the value
$\epsilon_{s}$. The dynamic matrix is non-diagonalizable but it can be brought
to Jordan canonical form by a non-unitary similarity transformation
$\check{U}$,

$\check{U}\check{H}\check{U}^{-1}=\begin{bmatrix}\ddots&\phantom{0}&\phantom{0}&\phantom{0}&\phantom{0}&\phantom{0}\\\
\phantom{0}&\epsilon_{s}&1&\phantom{0}&\phantom{0}&\phantom{0}\\\
\phantom{0}&\phantom{0}&\epsilon_{s}&\ddots&\phantom{0}&\phantom{0}\\\
\phantom{0}&\phantom{0}&\phantom{0}&\ddots&1&\phantom{0}\\\
\phantom{0}&\phantom{0}&\phantom{0}&\phantom{0}&\epsilon_{s}&\phantom{0}\\\
\phantom{0}&\phantom{0}&\phantom{0}&\phantom{0}&\phantom{0}&\ddots\end{bmatrix}.$
(38)

In this basis, $\check{H}$ is almost diagonal except on the $M\times M$ block
for the eigenvalue $\epsilon_{s}$, for which there are factors of $1$ above
the upper diagonal. As a consequence, there are fewer than $2N$ total
eigenvectors of $\check{H}$. As a technical replacement for the missing
eigenvectors, it is convenient to introduce additional basis vectors of the
Jordan block. Letting $\ket{s,1}$ denote an eigenvector of $\check{H}$ for the
eigenvalue $\epsilon_{s}$, one may introduce $\ket{s,n}$ with $n\leq M$
defined through the relation,

$\check{H}\ket{s,n}=\epsilon_{s}\ket{s,n}+\ket{s,n-1}.$ (39)

The vectors $\ket{s,n}$ comprise a complete basis spanning the single-particle
Hilbert space.

In the presence of such a Jordan block the appearance the spectral Green’s
functions develop higher-order poles. In particular, for $s\leq
s_{1}<s_{2}\leq s+M$ in the same Jordan block, there appears the factor:

$\Big{(}\check{U}\frac{1}{\epsilon-\check{H}}\check{U}^{-1}\Big{)}_{s_{1}s_{2}}=\Big{(}\frac{1}{\epsilon-\epsilon_{s}}\Big{)}^{1+s_{2}-s_{1}}.$
(40)

Written as functions of the time $t$, these off-diagonal components possess
polynomial coefficients in front of the exponent
$|t|^{s_{2}-s_{1}}\exp(-i\epsilon_{s}t)$, resulting in transient algebraic
gain of certain initial correlations. This behavior does not survive away from
the exceptional point: generic perturbations restore the diagonalizability of
$\check{H}$. There is an extreme sensitivity to perturbations at an
exceptional point. This manifests in the analytic structure of the
eigenvalues, which develop fractional power law non-analyticities [80]. This
is anomalous compared to conventional, fully analytic Hermitian perturbation
theory. Consequences of these non-analyticities are explored in some of the
examples in section 3.

It is natural to wonder if there is some analytic signature of an exceptional
point present in the stationary density. Examining the Lyapunov equation eq.
(22), one can see the answer to be negative. Both the matrices $\check{H}$ and
$\check{D}$ are analytic functions in neighborhoods of exceptional points in
parameter space [80]. By expanding all terms in series and matching powers,
one can see that only integer powers are permitted for
$\check{F}_{\mathrm{st}}$. As a consequence, $\check{F}_{\mathrm{st}}$ and by
extension $\check{H}_{\mathrm{st}}$ are analytic functions on the parameter
space even at exceptional points. Thus, there will generically be no residual
signature of the anomalous nature of the dynamics left over at long times.

###  2.6 Fermions

In this section, the above formalism is adapted to study fermionic
Lindbladians. To begin, one needs a fermionic version of the Keldysh path
integral. This is obtained in essentially the same way as its bosonic
counterpart, though some additional care must be taken with respect to the
ordering of the anti-commuting fields. The fundamental building block for the
path integral is the fermionic coherent state defined by
$\hat{c}\ket{\psi}=\psi\ket{\psi}$, where $\psi$ is a complex Grassmann
number. The relevant overlap formula for the Lindbladian action on two sets of
Grassmann coherent states is the mirror of eq. (4),

$\bra{\psi^{+}_{2}}\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\mathcal{L}}}}(\ket{\psi^{+}_{1}}\bra{\psi^{-}_{1}})\ket{\psi^{-}_{2}}=-ie^{\bar{\psi}^{+}_{2}\psi^{+}_{1}+\bar{\psi}^{-}_{1}\psi^{-}_{2}}\mathcal{K}(\bar{\psi}^{+}_{2},\psi^{+}_{1},\bar{\psi}^{-}_{1},\psi^{-}_{2}),$
(41)

where the Keldysh Hamiltonian is given by eq. (5) with Grassmann fields
$\psi^{\pm}$ in place of the bosonic fields. Note that in eq. (5b) the
ordering of fields in the first term is non-trivial, chosen so that backwards
fields $(\bar{\psi}^{-},\psi^{-})$ always appear before the forwards fields
$(\bar{\psi}^{+},\psi^{+})$ in the dissipative term $\mathcal{D}$. With this,
one can massage the partition function eq. (3) into the form of a fermionic
functional integral. It is standard to use the Larkin-Ovchinnikov convention,
in which the Keldysh-rotated fields defined
$\psi^{1,2}=(\psi^{+}\pm\psi^{-})/\sqrt{2}$ and
$\bar{\psi}^{1,2}=(\bar{\psi}^{+}\mp\bar{\psi}^{-})/\sqrt{2}$. The resulting
partition function is:

$Z=\int\text{D}\bar{\psi}^{a}\text{D}\psi^{a}e^{iS[\bar{\psi}^{a},\psi^{a}]},$
(42a)
$S=\int\text{d}t\Big{(}\bar{\psi}^{1}i\partial_{t}\psi^{1}+\bar{\psi}^{2}i\partial_{t}\psi^{2}-\mathcal{K}(\bar{\psi}^{a},\psi^{a})\Big{)}.$
(42b)

Grouping all four fields together into the Keldysh-Nambu vector
$\Psi=[\Psi^{1}\ \Psi^{2}]$, where $\Psi^{1}=[\psi^{1}\ \bar{\psi}^{2}]$ and
$\Psi^{2}=[\psi^{2}\ \bar{\psi}^{1}]$, the matrix of two-point functions
defines the fermion Green’s functions:

$\begin{bmatrix}i\check{G}^{\mathrm{R}}(t,t^{\prime})&i\check{G}^{\mathrm{K}}(t,t^{\prime})\\\
0&i\check{G}^{\mathrm{A}}(t,t^{\prime})\end{bmatrix}=\braket{\Psi(t)\bar{\Psi}(t^{\prime})}.$
(43)

These play the same role as in the bosonic theory.

A quadratic Lindbladian system of $N$ fermions is defined by the Hamiltonian
and jump operators:

$\hat{\mathcal{H}}=\sum_{i,j}^{N}\big{(}\varepsilon_{ij}\hat{c}^{\dagger}_{i}\hat{c}_{j}+\Delta_{ij}\hat{c}_{i}\hat{c}_{j}+\Delta_{ij}^{*}\hat{c}^{\dagger}_{j}\hat{c}^{\dagger}_{i}\big{)},$
(44a)
$\hat{\mathcal{L}}_{v}=\sum_{j}^{N}(\mu_{vj}\hat{c}_{j}+\nu_{vj}\hat{c}^{\dagger}_{j}),$
(44b)

where the $\hat{c}_{j}$’s are fermion creation operators,
$\\{\hat{c}_{i},\hat{c}_{j}^{\dagger}\\}=\delta_{ij}$. The corresponding
Keldysh action is:

$S=\frac{1}{2}\int\text{d}t\bar{\Psi}\begin{bmatrix}i\partial_{t}-\check{H}_{0}+i\check{Q}&i\check{D}\\\
0&i\partial_{t}-\check{H}_{0}-i\check{Q}\end{bmatrix}\Psi,$ (45)

The operators $\check{H}_{0}$, $\check{Q}$, and $\check{D}$ are Hermitian
$2N\times 2N$ matrices in Nambu space:

$\check{H}_{0}=\begin{bmatrix}\varepsilon&2\Delta^{\dagger}\\\
2\Delta&-\varepsilon^{\mathrm{T}}\end{bmatrix},$ (46a)
$\check{Q}=\frac{1}{2}\begin{bmatrix}\gamma+\tilde{\gamma}&\eta_{\mathrm{s}}^{\dagger}\\\
\eta_{\mathrm{s}}&\gamma^{\mathrm{T}}+\tilde{\gamma}^{\mathrm{T}}\end{bmatrix},$
(46b)
$\check{D}=\begin{bmatrix}\gamma-\tilde{\gamma}&\eta_{\mathrm{a}}^{\dagger}\\\
\eta_{\mathrm{a}}&\tilde{\gamma}^{\mathrm{T}}-\gamma^{\mathrm{T}}\end{bmatrix},$
(46c)

where $\gamma$, $\tilde{\gamma}$, and $\eta$ are defined the same as in the
bosonic theory as per eq. (13). The $N\times N$ parameter matrices have the
index symmetries:

$\varepsilon=\varepsilon^{\dagger},\quad\Delta=-\Delta^{\mathrm{T}},\quad\gamma=\gamma^{\dagger},\quad\tilde{\gamma}=\tilde{\gamma}^{\dagger},\quad\eta_{\mathrm{s,a}}=\pm\eta_{\mathrm{s,a}}^{\mathrm{T}}.$
(47)

The matrix $\check{H}_{0}$ is the single-particle Hamiltonian,
$\hat{\mathcal{H}}=\frac{1}{2}\hat{\mathcal{C}}^{\dagger}\check{H}_{0}\hat{\mathcal{C}}$,
where $\hat{\mathcal{C}}=[\hat{c}\ \hat{c}^{\dagger}]$ is the Nambu space
vector of fermion creation operators.

Note that the definitions of the matrix $\check{Q}$ and $\check{D}$ are
reversed compared to the bosonic theory. The consequence is that dynamical
stability is guaranteed for all choices of the parameters. To see why, define
the fermion single-particle dynamic matrix
$\check{H}=\check{H}_{0}-i\check{Q}$. Analogous to the bosonic theory,
dynamical stability of the theory requires all eigenvalues of this matrix to
have non-positive imaginary parts. To see why this is guaranteed, first
introduce the Nambu space vectors $\ket{v_{1v}}=[\mu_{v}^{*}\ \nu^{*}_{v}]$
and $\ket{v_{2v}}=[\nu_{v}\ \mu_{v}]$. Then
$\check{Q}=\frac{1}{2}\sum_{a,v}\ket{v_{av}}\bra{v_{av}}$ is just a sum of
one-dimensional projection operators. As a consequence, all of its eigenvalues
are strictly non-negative. Now consider an eigenvector $\ket{\epsilon}$ of
$\check{H}$. Then the imaginary part of the corresponding eigenvector is
$\mathrm{Im}(\epsilon)=-\braket{\epsilon}{\check{Q}}{\epsilon}\leq 0$. The
physical reason behind this is Pauli exclusion: there is a limit to the number
of fermions each state can hold, preventing an uncontrolled number of
particles from entering the system even when the rate of particle pumping is
greater than the rate of loss.

As for bosons, the many-body Lindbladian eigenvalues are given by assigning
occupation numbers to each eigenvalue of $\check{H}$ via eq. (2). For
fermions, this must be done with the understanding that the $n_{s}$ are
fermionic occupation numbers, equal only to either 0 or 1.

The fermionic spectral Green’s functions are given by:

$\check{G}^{\mathrm{R}}(t-t^{\prime})=-i\theta(t-t^{\prime})e^{-i(t-t^{\prime})\check{H}},\quad\quad\check{G}^{\mathrm{A}}(t-t^{\prime})=i\theta(t^{\prime}-t)e^{-i(t-t^{\prime})\check{H}^{\dagger}}.$
(48)

In frequency space, this becomes the resolvent of $\check{H}$:

$\check{G}^{\mathrm{R}}(\epsilon)=\frac{1}{\epsilon-\check{H}},\quad\quad\check{G}^{\mathrm{A}}(\epsilon)=\frac{1}{\epsilon-\check{H}^{\dagger}}.$
(49)

The poles of the Green’s functions are located at the complex eigenvalues of
$\check{H}$.

The Keldysh Green’s function is parametrized through the distribution matrix
$\check{F}(t,t^{\prime})$ through:

$\check{G}^{\mathrm{K}}=\check{G}^{\mathrm{R}}\circ\check{F}-\check{F}\circ\check{G}^{\mathrm{A}}.$
(50)

The distribution matrix obeys the quantum kinetic equation:

$[\partial_{t}\overset{\circ}{,}\check{F}]=-i\big{(}\check{H}\check{F}-\check{F}\check{H}^{\dagger}\big{)}+\check{D},$
(51)

where $\check{D}$ and $\partial_{t}$ come with factors of
$\delta(t-t^{\prime})$. In the stationary limit, one has
$\check{F}(t,t^{\prime})=\check{F}_{\mathrm{st}}\delta(t-t^{\prime})$ where
$\check{F}_{\mathrm{st}}$ obeys the Lyapunov equation,

$0=-i\big{(}\check{H}\check{F}_{\mathrm{st}}-\check{F}_{\mathrm{st}}\check{H}^{\dagger}\big{)}+\check{D}.$
(52)

Paralleling eq. (23), in the eigenbasis of $\check{H}$, one may solve for the
components of $\check{F}_{s}s$ as:

$\big{(}\check{U}\check{F}_{\mathrm{st}}\check{U}^{\dagger}\big{)}_{ss^{\prime}}=\frac{-i}{\epsilon_{s}-\epsilon^{*}_{s^{\prime}}}\big{(}\check{U}\check{D}\check{U}^{\dagger}\big{)}_{ss^{\prime}}.$
(53)

The stationary distribution matrix $\check{F}_{\mathrm{st}}$ is equivalent to
the stationary equal time Keldysh Green’s function, which itself is equal to
the fermion covariance matrix,
$\check{F}_{\mathrm{st}}=i\check{G}^{\mathrm{K}}(t,t)=\braket{[\hat{\mathcal{C}},\hat{\mathcal{C}}^{\dagger}]}$.

When the stationary state is unique, the stationary density matrix is a
fermionic Gaussian state of the form of eq. (25a), with

$\hat{\mathcal{H}}_{\mathrm{st}}=\frac{1}{2}\hat{\mathcal{C}}^{\dagger}\check{H}_{\mathrm{st}}\hat{\mathcal{C}}.$
(54)

The effective Hamiltonian is related to the distribution matrix through [60]:

$\check{F}_{\mathrm{st}}=\tanh(\check{H}_{\mathrm{st}}/2).$ (55)

Diagonalizing $\check{F}_{\mathrm{st}}$, one finds $2N$ eigenvalues of
$\check{H}_{\mathrm{st}}$ that determine the population numbers of
$\rho_{\mathrm{st}}$ in its diagonal basis [62]:

$\check{U}_{F}\check{F}\check{U}^{-1}_{F}=\mathrm{diag}\big{(}\tanh(\beta_{1}/2),...,-\tanh(\beta_{1}/2),...\big{)}.$
(56a)
$\hat{\mathcal{H}}_{\mathrm{st}}=\sum_{j}^{N}\beta_{j}\hat{d}_{j}^{\dagger}\hat{d}_{j},$
(56b)

with $[\hat{d}\ \hat{d}^{\dagger}]=\check{U}[\hat{c}\ \hat{c}^{\dagger}]$. In
contrast to the bosonic theory, there are no bounds on the inverse effective
temperatures $\beta_{j}$; they may be negative.

The expectations of observables, in the presence of sources, can be evaluated
by solving the fermionic analogue eq. (35) and using that of eq. (36),

$\partial_{t}\check{F}(t)=-i\big{(}\check{H}(t)\check{F}(t)-\check{F}(t)\check{H}^{\dagger}(t)\big{)}+\check{D}(t).$
(57a)
$\braket{\hat{\mathcal{O}}}(t)=\frac{1}{2}\operatorname{tr}\Big{(}\check{O}\big{(}1-\check{F}(t)\big{)}\Big{)},$
(57b)
$\braket{\mathcal{O}^{\mathrm{c}}(t)}=-\frac{1}{2}\operatorname{tr}\big{(}\check{O}\check{F}(t)\big{)}.$
(57c)

For weak perturbations from the stationary state, the linear response formulas
which replace eq.s (32), (33), and (34) are:

$\braket{\mathcal{O}^{\mathrm{c}}(t)\delta_{H_{0}}\mathcal{K}(t^{\prime})}_{\mathrm{st}}=\frac{1}{2}\operatorname{tr}\Big{(}\check{G}^{\mathrm{R}}(t,t^{\prime})\big{[}\check{F}_{\mathrm{st}},\delta\check{H}_{0}(t^{\prime})\big{]}\check{G}^{\mathrm{A}}(t^{\prime},t)\check{O}\Big{)},$
(58a)
$\braket{\mathcal{O}^{\mathrm{c}}(t)\delta\mathcal{K}_{Q}(t^{\prime})}_{\mathrm{st}}=\frac{i}{2}\operatorname{tr}\Big{(}\check{G}^{\mathrm{R}}(t,t^{\prime})\big{\\{}\check{F}_{\mathrm{st}},\delta\check{Q}(t^{\prime})\big{\\}}\check{G}^{\mathrm{A}}(t^{\prime},t)\check{O}\Big{)},$
(58b)
$\braket{\mathcal{O}^{\mathrm{c}}(t)\delta\mathcal{K}_{D}(t^{\prime})}_{\mathrm{st}}=-\frac{i}{2}\operatorname{tr}\Big{(}\check{G}^{\mathrm{R}}(t,t^{\prime})\delta\check{D}(t^{\prime})\check{G}^{\mathrm{A}}(t^{\prime},t)\check{O}\Big{)}.$
(58c)

##  3 Examples

Below, the formalism presented above is developed to study examples of
Lindbladian band theory, semiclassical kinetics, and mean-field theory. While
by no means an exhaustive list, these examples demonstrate how both quadratic
and nonlinear Lindbladians may be studied using the Keldysh language and serve
to illustrate important differences compared to the equilibrium theory.

###  3.1 Parametrically Driven Oscillator

As a warm up, consider first a model with one degree of freedom: a single
linear bosonic oscillator in contact with a thermal bath and subjected to a
parametric drive. This simple model is prototypical of much of the phenomena
unique to Lindbladian dynamics described above. The Hamiltonian and jump
operators for the parametrically driven oscillator in the rotating frame of a
drive are given by:

$\hat{\mathcal{H}}=\Delta\hat{a}^{\dagger}\hat{a}+\lambda(\hat{a}^{\dagger
2}+\hat{a}^{2}),$ (59a)
$\hat{\mathcal{L}}_{1}=\hat{a},\qquad\hat{\mathcal{L}}_{2}=\hat{a}^{\dagger}.$
(59b)

The jump operators describe the loss and gain of quanta to and from the
environment at the corresponding rates $\gamma$ and $\tilde{\gamma}$. The
rates can be related to the strength of the coupling between the system and
the bath $\kappa$, the bath temperature $T$, and the natural frequency of the
system $\omega_{0}$ by:

$2\kappa=\gamma-\tilde{\gamma},\qquad\coth(\omega_{0}/2T)=\frac{\gamma+\tilde{\gamma}}{\gamma-\tilde{\gamma}}.$
(60)

Note the appearance of the Bose function
$2n_{\mathrm{B}}+1=\coth(\omega_{0}/2T)$ at the bath temperature and system
frequency.

The dynamic matrix $\check{H}$ is given by:

$\check{H}=\begin{bmatrix}\Delta-i\kappa&2\lambda\\\
-2\lambda&-\Delta-i\kappa\end{bmatrix}.$ (61)

The spectral Green’s functions are given by eq. (19),

$\check{G}^{\mathrm{R},\mathrm{A}}(\epsilon)=\frac{1}{(\epsilon\pm
i\kappa)^{2}-\Omega^{2}}\begin{bmatrix}\epsilon+\Delta\pm i\kappa&2\lambda\\\
2\lambda&\epsilon-\Delta\pm i\kappa\end{bmatrix},$ (62)

where $\Omega^{2}=\Delta^{2}-4\lambda^{2}$ is the Bogoliubov frequency of the
Hamiltonian part of $\check{H}$. The eigenvalues of $\check{H}$ match the
poles of the spectral Green’s functions from the numerator in the above
expression,

$\epsilon_{1,2}=-i\kappa\mp\Omega.$ (63)

From this, one can see that the system is stable so long as the coupling
$\kappa$ is positive, which occurs when the rate of loss of quanta is greater
than the rate of gain, $\gamma>\tilde{\gamma}$. Additionally, the model is
only stable when the drive strength is small enough,
$2\lambda\leq\sqrt{\Delta^{2}+\kappa^{2}}$. When the bound is saturated, one
of the two eigenvalues is tuned to zero and the theory is at the brink of
instability. The matrix $\check{U}$ performing the diagonalization is
equivalent to the coherent Bogoliubov rotation matrix in the absence of the
coupling to the bath:

$\check{U}=\frac{1}{\sqrt{2}}\begin{bmatrix}\sqrt{\frac{\Delta}{\Omega}+1}&\sqrt{\frac{\Delta}{\Omega}-1}\\\
\sqrt{\frac{\Delta}{\Omega}-1}&\sqrt{\frac{\Delta}{\Omega}+1}\end{bmatrix}$
(64)

The matrix $\check{D}$ is proportional to the identity matrix by the factor
$2\kappa\coth(\omega_{0}/2T)$. The Keldysh Green’s function is given by:

$\check{G}^{\mathrm{K}}(\epsilon)=\frac{-2i\kappa\coth(\omega_{0}/2T)}{\big{(}(\epsilon+i\kappa)^{2}-\Omega^{2}\big{)}\big{(}(\epsilon-i\kappa)^{2}-\Omega^{2}\big{)}}\\\
\times\begin{bmatrix}(\epsilon+\Delta)^{2}+\kappa^{2}+4\lambda^{2}&-4\lambda(\Delta+i\kappa)\\\
-4\lambda(\Delta-i\kappa)&(\epsilon-\Delta)^{2}+\kappa^{2}+4\lambda^{2}\end{bmatrix}.$
(65)

Solving eq. (21) gives the stationary distribution matrix:

$\check{F}_{\mathrm{st}}=\frac{\coth(\omega_{0}/2T)}{\kappa^{2}+\Omega^{2}}\begin{bmatrix}\Delta^{2}+\kappa^{2}&-2\lambda(\Delta+i\kappa)\\\
-2\lambda(\Delta-i\kappa)&\Delta^{2}+\kappa^{2}\end{bmatrix}.$ (66)

The eigenvalues of $\check{F}_{\mathrm{st}}$ are related to the diagonal
frequency $\beta$ of the effective Hamiltonian
$\hat{\mathcal{H}}_{\mathrm{st}}$ through eq. (24),

$\coth(\beta/2)=\coth(\omega_{0}/2T)\sqrt{\frac{\kappa^{2}+\Delta^{2}}{\kappa^{2}+\Omega^{2}}}.$
(67)

The diagonalizing transformation $\check{U}_{F}$ is given by:

$\check{U}_{F}=\frac{1}{\sqrt{2}}\begin{bmatrix}e^{-i\theta}\sqrt{\sqrt{\frac{\kappa^{2}+\Delta^{2}}{\kappa^{2}+\Omega^{2}}}+1}&e^{i\theta}\sqrt{\sqrt{\frac{\kappa^{2}+\Delta^{2}}{\kappa^{2}+\Omega^{2}}}-1}\\\
e^{-i\theta}\sqrt{\sqrt{\frac{\kappa^{2}+\Delta^{2}}{\kappa^{2}+\Omega^{2}}}-1}&e^{i\theta}\sqrt{\sqrt{\frac{\kappa^{2}+\Delta^{2}}{\kappa^{2}+\Omega^{2}}}+1}\end{bmatrix},$
(68)

where $\theta=\arg(\Delta+i\kappa)$. As discussed above, the bases in which
the distribution and dynamic matrices are diagonal are not the same,
$\check{U}\neq\check{U}_{F}$.

To gain more intuition for the model, consider the limits of strong versus
weak dissipation. When the dissipation is very weak compared to all other
scales in the problem $\kappa\to 0$, the two Lindbladian eigenvalues
$-i\epsilon_{s}$ are complex conjugates with a small negative real part. The
dynamics in this limit are weakly under-damped coherent rotation. Up to
corrections of order $\kappa$ the stationary state is diagonal in the basis
Bogoliubov quasi-particles of the coherent problem. That is,
$\check{U}_{F}=\check{U}+O(\kappa)$. With this, the stationary state effective
Hamiltonian in its diagonal basis is given by:

$\hat{\mathcal{H}}_{\mathrm{st}}=\beta\hat{b}^{\dagger}\hat{b},$ (69)

with $[\hat{b}\ \hat{b}^{\dagger}]=\check{U}_{F}[\hat{a}\ \hat{a}^{\dagger}]$.
This is equal to the original Hamiltonian $\hat{\mathcal{H}}$ up to a
proportionality constant $\beta=\Omega/T_{\mathrm{eff}}$,

$\coth(\Omega/2T_{\mathrm{eff}})=\frac{\Delta}{\Omega}\coth(\omega_{0}/2T).$
(70)

The stationary state is thus a thermal state of the coherent dynamics, but
with a different effective temperature than the bath temperature. Note that
even in the limit $T\to 0$, the effective temperature $T_{\mathrm{eff}}$ is
finite. This phenomenon is known as quantum heating [23, 24, 25].

Alternatively, one may consider the the limit of strong dissipation,
$\kappa\to\infty$. In this limit, the drive strength can stably be larger than
the detuning, $2\lambda>\Delta$, past the point of coherent stability. In this
situation, the Lindbladian eigenvalues $-i\epsilon_{s}$ are purely imaginary
and the dynamics are that of over-damped pure dissipation without any coherent
rotation. Up to corrections of the order of the Hamiltonian parameters
$\Delta$ and $\lambda$, the dynamic matrix $\check{H}$ is proportional to the
identity with a factor of $\kappa$. Thus, the diagonal basis for the dynamics
is the starting basis of the problem, $\check{U}=\check{1}$. Similarly,
examining eq. (66) in this limit, one sees that $\check{F}_{\mathrm{st}}$ is
also proportional to the identity, so that $\check{U}_{F}=\check{1}$. The
stationary state effective Hamiltonian is given by:

$\hat{\mathcal{H}}_{\mathrm{st}}=\frac{\omega_{0}}{T}\,\hat{a}^{\dagger}\hat{a}.$
(71)

The stationary state is a thermal state of the un-driven Hamiltonian
$\omega_{0}\hat{a}^{\dagger}\hat{a}$ with a temperature equal to the bath
temperature.

Thus, tuning the dissipation strength between the extreme limits the extreme
limits $\kappa\to 0$ and $\kappa\to\infty$ changes the diagonal basis of both
the dynamics and stationary state between the Bogoliubov basis of the
$(\hat{b},\hat{b}^{\dagger})$ bosons and the original basis of the
$(\hat{a},\hat{a}^{\dagger})$ bosons. In both of these limits, these basis are
the same, $\check{U}=\check{U}_{F}$; this will cease the be the case in
between these two extremes. Interstitial between these two regimes, there is
an exceptional point at which the two Lindbladian eigenvalues $-i\epsilon_{s}$
coalesce to the value $-\kappa$ and the dynamics is a resonantly damped
dissipation. This occurs at the threshold of coherent instability
$2\lambda=\Delta$. At this point, the matrix $\check{H}$ is non-diagonalizable
and is brought to Jordan canonical form by the similarity transformation:

$\check{U}=\begin{bmatrix}0&-1\\\ 2\lambda&2\lambda\end{bmatrix}$ (72)

Note that this expression is not the equal to the limiting form of eq. (64),
which itself does not exist. In contrast, the distribution matrix is a smooth
function of the parameters even at this point, as are the diagonal frequencies
of the effective Hamiltonian and the diagonal basis of bosons determined by
$\check{U}_{F}$. These are given by their limiting forms in terms of the above
general expressions. The exceptional point is reflected in the limiting form
of the spectral Green’s functions as second-order pole:

$\check{G}^{\mathrm{R},\mathrm{A}}(\epsilon)=\frac{1}{(\epsilon\pm
i\kappa)^{2}}\begin{bmatrix}\epsilon+2\lambda\pm i\kappa&2\lambda\\\
2\lambda&\epsilon-2\lambda\pm i\kappa\end{bmatrix}.$ (73)

The non-diagonalizability results in the linear gain of certain non-stationary
densities.

To exemplify this, consider increasing $\Delta$ slightly above or below the
critical value of $2\lambda$. This removes the eigenvalue degeneracy, thus
eliminating the polynomial gain of any initial correlations. The resulting
decay at rates are determined by differences of the perturbed eigenvalues
$\epsilon_{s}-\epsilon_{s^{\prime}}^{*}$. Writing $\Delta-2\lambda=\delta$,
expanding eq. (63) around the exceptional point $\delta=0$ gives a series in
powers of $\delta^{1/2}$:

$\epsilon_{1,2}=-i\kappa\mp 2\sqrt{\lambda\delta}+O(\delta).$ (74)

As a consequence, introducing a small detuning causes initial densities with
off-diagonal terms in the perturbed basis of $\check{H}$ eigenvalues to
coherently rotate at a rate $4\sqrt{\lambda\delta}$ that is a non-analytic
function of the deviation $\lambda$. This demonstrates a stronger sensitivity
to perturbations at the exceptional point than the usual $O(\delta)$
corrections in Hermitian systems. Further implications of the square root
singularity in eigenvalues near exceptional points is explored in the
following section.

###  3.2 Non-Hermitian Band Theory

This section examines two simple Lindbladian tight binding models. For
simplicity the focus is restricted to one-dimensional chains, though the
general principles discussed hear can naturally be extended to higher
dimensions. Much like clean coherent models, Lindbladian lattice systems with
translational invariance are described in terms of the band theory. The bands
in a Lindbladian system are the momentum-dependent eigenvalues of the dynamic
matrix $\check{H}$, and as such are generically complex. Beyond having an
imaginary part, there are additional complications that emerge due to the non-
Hermiticity of $\check{H}$, specifically relating to the potential existence
of exceptional points. This subtlety is showcased in two simple models below.
Note that the theory of non-Hermitian bands in connection with Lindbladian
dynamics is still under construction and the following discussion is far from
exhaustive, see [81, 82, 83, 84, 85, 86, 87].

As a first example, consider a chain of identical parametric oscillators from
the preceding section coupled linearly to their nearest neighbors. The
Hamiltonian and jump operators are:

$\hat{\mathcal{H}}=\sum_{j=1}^{N}\Big{(}\Delta_{0}\hat{a}^{\dagger}_{j}\hat{a}_{j}+\tau(\hat{a}^{\dagger}_{j+1}\hat{a}_{j}+\hat{a}^{\dagger}_{j}\hat{a}_{j+1})+\lambda_{0}(\hat{a}_{j}^{2}+\hat{a}_{j}^{\dagger
2})\Big{)},$ (75a)
$\hat{\mathcal{L}}_{1j}=\hat{a}_{j},\qquad\hat{\mathcal{L}}_{2j}=\hat{a}^{\dagger}_{j}.$
(75b)

With periodic boundary conditions, it is convenient to change to the momentum
representation,

$\phi^{\alpha}(k)=\frac{1}{\sqrt{N}}\sum_{j}e^{ijk}\phi^{\alpha}_{j},$ (76)

where $k\in[0,2\pi)$ is the crystal momentum. In momentum space, writing
$\Phi^{\alpha}(k)=[\phi^{\alpha}(k)\ \bar{\phi}^{\alpha}(-k)]$ brings the
Keldysh action to the standard form of eq. (11), where the parameter matrices
are diagonal functions of $k$,

$S=\frac{1}{2}\int\text{d}t\text{d}k\
\bar{\Phi}\begin{bmatrix}0&\check{\tau}^{3}i\partial_{t}-\check{H}_{0}(k)-i\check{Q}(k)\\\
\check{\tau}^{3}i\partial_{t}-\check{H}_{0}(k)+i\check{Q}(k)&i\check{D}(k)\end{bmatrix}\Phi.$
(77)

The parameter matrices are given by a $k$-dependent version of eq. (12) with
$\eta=0$ and the other parameters,

$\Delta(k)=\Delta_{0}+2\tau\cos(k),\qquad\lambda(k)=\lambda_{0},$ (78a)
$\gamma(k)=2\kappa(n_{\mathrm{B}}+1),\qquad\tilde{\gamma}(k)=2\kappa
n_{\mathrm{B}},\qquad\eta(k)=0.$ (78b)

with $n_{\mathrm{B}}$ defined by $2n_{\mathrm{B}}+1=\coth(\omega_{0}/2T)$ and
$\kappa$, $\omega_{0}$, and $T$ the dissipative coupling, system natural
frequency, and bath temperature defined in eq. (60).

With this, one can read off the Green’s functions from the previous section
using eq.s (62) and (65) and replacing the appropriate quantities with their
$k$-dependent analogues. There are two bands of eigenvalues of $\check{H}(k)$,
given by an $k$-dependent version of eq. (63),

$\epsilon_{1,2}=-i\kappa\mp\sqrt{\big{(}\Delta_{0}+2\tau\cos(k)\big{)}^{2}-4\lambda_{0}^{2}}.$
(79)

This defines a pair of complex-valued bands. As $\lambda_{0}$ is tuned from
being small to large, the two bands go from being complex with a constant
imaginary part to purely imaginary. These two regimes correspond to completely
under-damped and over-damped dissipation. The transition between these two
limits occurs via an extended intermediate regime in which the bands touch and
both over- and under-damped dissipation occurs in different momentum ranges.
This is depicted in fig. (5). The points in the Brillouin zone where the bands
touch define the exceptional momenta $k_{\mathrm{EP}}$, here given by:

$k_{\mathrm{EP}}=\arccos\bigg{(}\frac{2\lambda_{0}-\Delta_{0}}{2\tau}\bigg{)},$
(80)

which has two solutions when $\Delta_{0}-2\tau<2\lambda_{0}<\Delta_{0}+2\tau$.
Unlike in a Hermitian band touching, the touching of complex bands occurs at
exceptional points in the parameter space of $\check{H}$. This results in a
resonant damping at the exceptional momenta $k_{\mathrm{EP}}$.

Figure 5: Plots of the real and imaginary parts of the complex bands eq. (79)
with $\kappa=5$, $\Delta=2$, $\tau=1/2$. The three columns show the bands in
the completely under-damped, mixed, and completely over-damped regime,
corresponding to $\lambda=0,1$, and $1.7$, going from left to right.

One can solve the kinetic equation for the stationary state by looking for
translationally invariant solutions. This reduces eq. (22) to a $k$-dependent
Lyapunov equation,

$0=-i\big{(}\check{H}(k)\check{F}_{\mathrm{st}}(k)-\check{F}_{\mathrm{st}}(k)\check{H}^{\dagger}(k)\big{)}+\check{\tau}^{3}\check{D}(k)\check{\tau}^{3}.$
(81)

The solution to this equation can be read off from the solution to the single
parametric oscillator in eq. (66) by appropriately replacing parameters by
their $k$-dependent generalization. Just like the single parametric
oscillator, the stationary distribution has no signature of the exceptional
points: there are no non-analyticities at the exceptional momenta. Response to
translationally-invariant perturbations can be computed using the formalism
from section 2.4 and replacing matrices with their $k$-dependent versions.
Variations that are not spatially homogeneous but which vary smoothly over
large distances can be studied using a semiclassical approach; this is
discussed in section 3.4.

The above example exemplifies that the degeneration of bands in Lindbladian
models occurs at exceptional points and thus differs from Hermitian band
touching. This naively suggests that as long as there are no spectral
degeneracies, the band theory of Lindbladian systems is similar to
conventional Hermitian systems save for the additional imaginary part of each
band. This intuition however is badly wrong. Due to the non-analytic structure
of eigenvalues near exceptional points, non-Hermitian bands can exhibit
unusual structures even if they do not cross. In the simplest situation with a
$2\times 2$ non-Hermitian matrix, an exceptional point corresponds to a
degeneration of two eigenvalues and, as demonstrated in the preceding
sections, leads to a square root singularity. Revolving around the exceptional
point in parameter space induces a monodromy in which the two eigenvalues are
exchanged upon a single winding rather than returning to where they started.
In a one-dimensional tight binding model with two bands, one can imagine
tuning a parameter past the regime of band touching so that the bands remain
intertwined and mutually wind around the exceptional point, exhibiting
monodromy over one period of the Brillouin zone.

This clearly does not occur in the chain of parametric oscillators discussed
above. One can construct a simple model that demonstrates this phenomenon
using a Lindbladian generalization of the SSH model. Similar models were
studied in [82, 83], though not in the context of the Lindbladian dynamics.
The SSH Hamiltonian is defined by assigning two species of fermions,
$\hat{c}_{\mathrm{A}}$ and $\hat{c}_{\mathrm{B}}$, to each site on a one-
dimensional chain. Hopping is allowed between species on-site and between
nearest-neighbors,

$\hat{\mathcal{H}}=\sum_{j=1}^{N}\Big{(}\varepsilon_{0}(\hat{c}_{\mathrm{A}j}^{\dagger}\hat{c}_{\mathrm{A}j}+\hat{c}_{\mathrm{B}j}^{\dagger}\hat{c}_{\mathrm{B}j})+\tau_{1}(\hat{c}_{\mathrm{A}j}^{\dagger}\hat{c}_{\mathrm{B}j}+\hat{c}_{\mathrm{B}j}^{\dagger}\hat{c}_{\mathrm{A}j})+\tau_{2}(\hat{c}^{\dagger}_{\mathrm{A}j+1}\hat{c}_{\mathrm{B}j}+\hat{c}^{\dagger}_{\mathrm{B}j}\hat{c}_{\mathrm{A}j+1})\Big{)}.$
(82)

In addition, consider on-site gain and loss of particles via a superposition
of the A and B fermions,

$\mathcal{L}_{1j}=\sqrt{2\kappa}(\hat{c}_{\mathrm{A}j}+e^{i\theta_{1}}\hat{c}_{\mathrm{B}j}),\qquad\mathcal{L}_{2j}=\sqrt{2\kappa}(\hat{c}_{\mathrm{A}j}^{\dagger}+e^{-i\theta_{2}}\hat{c}_{\mathrm{B}j}^{\dagger}),$
(83)

The restriction to purely loss and gain processes in combination with the lack
of pair-creation terms in the Hamiltonian gives the model an overall $U(1)$
symmetry, with the associated charge being the total number of particles. This
symmetry is weak [76], and so does not imply the conservation of particle
number in time. A consequence of this symmetry is that the anomalous off-
diagonal blocks of the Nambu space Green’s functions vanish. As such, it is
convenient to define by $G$ without a check the upper left block of
$\check{G}$ and similarly for other quantities.

The momentum space Keldysh action is then given by:

$S=\int\text{d}t\text{d}k\
\bar{\psi}\begin{bmatrix}i\partial_{t}-H_{0}(k)+iQ(k)&iD(k)\\\
0&i\partial_{t}-H_{0}(k)-iQ(k)\end{bmatrix}\psi,$ (84)

where $\psi=[\psi^{1}\ \psi^{2}]$, $\psi^{a}=[\psi^{a}_{\mathrm{A}}(k),\
\psi^{a}_{\mathrm{B}}(k)]$, and the un-checked parameter matrices are the
upper-left blocks of their Nambu-space counterparts in the main text eq. (45).
The parameter matrices are:

$H_{0}(k)=\begin{bmatrix}\varepsilon_{0}&\tau_{1}+\tau_{2}e^{ik}\\\
\tau_{1}+\tau_{2}e^{-ik}&\varepsilon_{0}\end{bmatrix},$ (85a)
$Q(k)=\kappa\begin{bmatrix}2&e^{-i\theta_{1}}+e^{-i\theta_{2}}\\\
e^{i\theta_{1}}+e^{i\theta_{2}}&2\end{bmatrix},$ (85b)
$D(k)=2\kappa\begin{bmatrix}0&e^{-i\theta_{1}}-e^{-i\theta_{2}}\\\
e^{i\theta_{1}}-e^{i\theta_{2}}&0\end{bmatrix}$ (85c)

Due to the symmetry in the problem it suffices to study only the two
eigenvalue bands of $H=H_{0}-iQ$; the other two single-particle eigenvalue
bands are given by their complex conjugates.

For a simple demonstration of the concept discussed above, one can choose
different phases for loss and gain, $\theta_{1}=\pi/2$ and $\theta_{2}=0$.
Fixing in addition $\tau_{1}=\kappa$, the eigenvalues of the dynamic matrix
$H(k)$ are given by:

$\epsilon_{1,2}(k)=\varepsilon_{0}-2i\kappa\mp\sqrt{\tau_{2}^{2}-\kappa^{2}+2\kappa\tau_{2}\cos(k)-2i\kappa\Big{(}\kappa+\tau_{2}\big{(}\cos(k)+\sin(k)\big{)}\Big{)}}$
(86)

From this expression one can see that for $\tau_{2}=\kappa$, the eigenvalues
degenerate to an exceptional point at $k_{\mathrm{EP}}=3\pi/2$. For
$\tau_{2}<\kappa$, the eigenvalues are periodic functions of $k$ and the two
bands are disconnected. For $\tau_{2}>\kappa$ however, the sign of the
argument of the square root develops a negative real part for $k$ near
$3\pi/2$. As a consequence, sweeping over the Brillouin zone smoothly
traverses from one branch of the square root Riemann sheet to another,
introducing the monodromy $\epsilon_{1,2}(k\to 2\pi^{-})=\epsilon_{2,1}(k\to
0^{+})$. These three scenarios are depicted in Fig. (6). In other words, the
square root singularity of the exceptional point effects the analytic
structure of the eigenvalue bands even when they do not collide with it. From
encircling of the exceptional point there are no longer two distinct
disconnected bands, but rather a single continuous band that double-covers the
Brillouin zone.

Figure 6: Three regimes of the complex bands of eq. (86) plotted as paths in
the complex plane. The solid points indicate the values of $\epsilon_{j}(0)$
and the arrows show the directed path traced out by varying $k$ from $0$ to
$2\pi$. The three regimes correspond to the two disconnected untwisted loops,
intersecting loops, and a single untwisted loop for $\tau_{2}$ less than,
equal to, and greater than $\kappa$ respectively.

Recent literature has understood this sort of atypical feature of non-
Hermitian band theory in the context of braids and knots [84, 85, 86, 87]. The
eigenvalues of the dynamic matrix $H(k)$ define paths in the complex plain
parametrized by $k$. Due to the periodicity of $k$, the paths must close and
so define loops, which can braid together and link up. This phenomenon does
not have a Hermitian analogue, as the eigenvalues of Hermitian matrices live
on the real line, a space of too low dimension for the knotting of paths. In
the model considered above, the transition through the exceptional point can
be understood as a transition from two disconnected unknots to one. In the
language of braids, this is a transition from two upbraided paths to a pair of
paths with a single twist, separated by a point where the paths collide.

In more general models on one-dimensional lattices, eigenvalues may braid
multiple times with one another in a more complicated fashion. This can result
in the paths of eigenvalues tracing out collections of linked knots. In the
language of braids, an $N$ band model in a given range of parameters will
correspond to an element of the braid group with $N$ generators. Transitions
between different braid group elements/knot configurations occurs via the
joining and separating of paths by moving through exceptional points. Higher-
dimensional generalizations of this phenomenon are harder to visualize, as
they entail the mapping of higher-dimensional Brillouin zones (tori) into the
complex plane. This problem has received some recent attention (see for
example [86]) but in general warrants future study.

As discussed above, the stationary solution to the kinetic theory,
$\check{F}_{\mathrm{st}}$, displays no signature of exceptional points. This
remains true of the model considered here: the stationary distribution, while
non-trivial, varies smoothly across the crossing of the exceptional point and
does not differ qualitatively on either side of the transition. The full
expressions for the bands $\beta_{j}(k)$ of stationary eigenvalues of
$\check{F}_{\mathrm{st}}(k)$ via eq. (55) are somewhat cumbersome and so are
not reproduced here; they are plotted below, at, and above the transition in
Fig. (7). There are no apparent signatures these topological features present
in the different regimes in the stationary distribution.

Figure 7: Plot of the two bands of the bands of stationary eigenvalues
$\beta_{j}(k)$ of the Lindbladian SSH chain, showing the three different
dynamical regimes. The non-analytic structure of the dynamic eigenvalues
$\epsilon_{j}(k)$ is not retained by the stationary distribution; at all
values of $k$, the curves vary smoothly across the transition.

###  3.3 Disordered Fermions

This section examines a simple model of disordered Lindbladian fermions. A
$U(1)$ symmetric model of $N$ flavors of fermions that are otherwise
featureless is studied. In contrast to the preceding sections discussing
‘clean’ systems, the model considered here gives insight into the behaviour of
generic quadratic Lindbladians. The $U(1)$ symmetry provides a meaningful
distinction between loss and gain even in the absence of additional symmetry.
In preparing this manuscript, a paper [88] appeared up on arXiv which
discusses ideas very similar to those presented here. They focus on Majorana
fermions instead of the Dirac fermions considered here and provide a more
detailed discussion of the spectrum and level statistics of both dynamic
matrix and steady state.

The model examined here can be compared to various studies on Lindbladians in
which the Hamiltonian and jump operators are random matrices, which are aimed
at understanding generic Lindbladian dynamics in the absence of any additional
structure [89, 90, 91, 92, 93, 94]. As will be shown below, in certain limits
the single-particle quantities of the random quadratic Lindbladian match the
pure random matrix results, but generically they different. This runs counter
to the usual intuition from coherent systems. For equilibrium disordered
fermions the many-body spectrum is determined by the single-particle
Hamiltonian, which in turn is characterized by Hermitian random matrix theory.
Solving the many-body problem is achieved by solving the pure random matrix
quantum mechanics. In the Lindbladian framework, the single-particle dynamic
matrix and stationary distribution are separate quantities and do not define
any sort of ‘single-particle Lindbladian.’ In this sense, the random quadratic
problem is not equivalent to the random matrix Lindblad problem.

The Hamiltonian and jump operators defining the model are:

$\hat{\mathcal{H}}=\sum_{i,j}^{N}(H_{0})_{ij}\hat{c}^{\dagger}_{i}\hat{c}_{j},$
(87a)
$\hat{\mathcal{L}}_{1v}=\sum_{j}^{N}\mu_{vj}\hat{c}_{j},\qquad\hat{\mathcal{L}}_{2v}=\sum_{j}^{N}\nu_{vj}\hat{c}^{\dagger}_{j},$
(87b)

where there are a total of $M$ jump operators, $v\leq M$. The parameters
$H_{0}$, $\mu$, and $\nu$ are be Gaussian random variables. The disorder
averaging is defined by:

$[...]=\int\text{D}H_{0}\text{D}\mu\text{D}\nu
e^{-\frac{N}{2}\operatorname{tr}(H_{0}^{2})-\frac{M}{\gamma}\sum_{j,v}|\mu_{vj}|^{2}-\frac{M}{\tilde{\gamma}}\sum_{j,v}|\nu_{vj}|^{2}}(...).$
(88)

Of interest here is the large $N$ limit where the ratio $m=M/N$ is held
finite. In this limit, the model has a total of three parameters: $m$ and the
two parameters that control the ratio of loss and gain vs. the strength of the
Hamiltonian, $\gamma$ and $\tilde{\gamma}$.

In the corresponding coherent model, one solves the random matrix theory
problem of the single-particle Hamiltonian. This is exactly solvable in the
large $N$ limit, with the density of states
$\varrho(\epsilon)\equiv[\operatorname{tr}\delta(\epsilon-H_{0})]$ given by
the well-known Wigner semicircle distribution. In the Lindbladian setting, one
can determine the spectrum by solving the non-Hermitian random matrix problem
for the eigenvalue distribution
$\varrho_{H}(\epsilon)=[\operatorname{tr}\delta(\epsilon-H)]$, which is the
probability distribution of the single particle eigenvalues $\epsilon_{j}$ in
the complex $\epsilon$ plain. This problem turns out to involve more subtle
features than its Hermitian counterpart. In the large $N$ limit, the limit
shape of the distribution changes shape as a function of the model parameters,
with different phases defined by the number of connected components. In
addition to the spectrum, one can also analyze the stationary state, in which
the object of interest is the density of states of the stationary effective
Hamiltonian, $\varrho_{\mathrm{st}}(\beta)=[\operatorname{tr}\delta(\epsilon-
H_{\mathrm{st}})]$, which is the probability distribution of the real-valued
$\beta_{j}$. For this one must solve a random Lyapunov equation. Compared to
the spectral random matrix problem, this problem is hard to solve exactly,
even in the large $N$ limit. Simple numerical computations show that the
transitions in the dynamics are mirrored by transitions in the stationary
state.

Consider first the simple case of pure random loss $\gamma\to\infty$, so that
the $H_{0}$ and $\nu$ are dropped. The single-particle dynamic matrix is a
Wishart matrix, given by a sum of one-dimensional projection operators:

$H_{ij}=-i\sum_{v}^{M}\mu^{*}_{vi}\mu_{vj}.$ (89)

In the large $N$ limit, its eigenvalue distribution is given by the well-known
Marchenko–Pastur law [95], with the eigenvalue distribution
$\varrho_{H}(\epsilon)$ given by,

$\varrho_{H}(\epsilon)=-i\Bigg{(}\frac{m}{\gamma\epsilon}\sqrt{\frac{\gamma^{2}}{m}-\frac{1}{4}\bigg{(}\epsilon-\frac{\gamma}{m}(1+m)\bigg{)}^{2}}+\pi(1-m)\theta(1-m)\delta(\epsilon)\Bigg{)},$
(90)

where $\theta$ is the Heaviside function. When $m>1$, all eigenvalues have
finite imaginary part and the stationary state is unique, being given by the
Fock vacuum state with zero particles. When $m<1$, there are always a subset
of modes that do not dissipate, reflected in $\varrho_{H}(\epsilon)$ by the
delta measure at the origin. The stationary state is not unique. The critical
point $m=1$ divides the two regimes. At this point, the dissipative gap closes
and the uniqueness of the stationary state breaks down. Note that for all
values of $m$, the eigenvalues are purely imaginary, so that
$\varrho_{H}(\epsilon)$ is supported only on the imaginary axis. This is
qualitatively different from the results reported for general disordered
Lindbladians with random matrix jump operators and no Hamiltonian, in which
the distribution of Lindbladian eigenvalues occupies a lemon shaped region in
the complex plane with a finite area [90].

Adding in the Hamiltonian part but keeping $\nu=0$, the Marchenko-Pastur
distribution is deformed into the complex plane. The distribution
$\varrho_{H}$ is supported on compact subsets of the complex $\epsilon$ plane
with a finite area. For large enough $\gamma$, the large and small $m$ phases
deform into phases in which the support of $\varrho_{H}(\epsilon)$ has one and
two connected components respectively. Example spectra in the different phases
are shown in Fig. (8). The shapes of the eigenvalue distribution can be
determined for large $N$ using non-Hermitian random matrix theory. Similar
problems have been studied in the context of chaotic scattering [96, 97, 98];
for a more detailed discussion of the shape of the distribution for the
Lindbladian problem one may refer to [88]. With the Hamiltonian part, all
eigenvalues now have a finite imaginary part for all $m$. In contrast to the
above situation without a Hamiltonian discussed above, with a Hamiltonian the
dissipative gap becomes finite even for small $m$. As a consequence, the
stationary state is always unique. It is easy to check that the stationary
state is still the zero particle Fock vacuum for all values of $m$.

Figure 8: Example spectra of $H$ for different values of $m$, with $N=500$ and
$\gamma=1$ and $\tilde{\gamma}=0$. The leftmost plot depicts the small $m$
phase in which there are two disconnected connected components of
$\varrho_{H}(\epsilon)$. Note that in the small $m$ limit, one of the two
regions of eigenvalues is very close to real axis, but still possesses a
finite imaginary part. The second to the left is near the critical $m$ (which
generically differs from 1 with a Hamiltonian term), in which the two regions
are starting to merge. The two on the right show the large $m$ phase, for
which $\varrho_{H}(\epsilon)$ possesses a single connected component. The
rightmost plot shows the large $m$ limit, in which the limiting shape of the
support of $\varrho_{H}(\epsilon)$ can be seen to approach a Ginibre disk
centered around $1/2$ on the imaginary axis. This limit matches the known
behaviour for random matrix Lindbladians with both Hamiltonian and jump
operators [90].

Turning finally to the general model with gain and loss, one finds similar
geometric transitions depending on the value of $m$. There are multiple
different phases with a maximum of three connected components, which can merge
and split in various ways depending on the relative strengths of gain vs. loss
vs. Hamiltonian. The stationary state is generically a non-trivial mixed
state. For large $m$, $\varrho_{\mathrm{st}}(\beta)$ is a Wigner semicircle,
with a width that scales with $1/m$ and is off-centered from zero by an amount
determined by $\tilde{\gamma}/\gamma$. For smaller $m$,
$\varrho_{\mathrm{st}}(\beta)$ has support on multiple disconnected regions of
the real line, each of which deviates from the semi-circle law. A detailed
classification of the various phases is not presented here, but several
different examples of $\varrho_{M}$ and the corresponding
$\varrho_{\mathrm{st}}(\beta)$ are depicted in Fig. (9).

Figure 9: Example spectra of $H$ plotted above the corresponding
$\varrho_{\mathrm{st}}(\beta)$ for finite $\gamma$ and $\tilde{\gamma}$, with
$N=500$. The leftmost plot depicts a phase with in which the support of
$\varrho_{H}(\epsilon)$ has three connected components, which occurs when $m$
is small and gain is meaningfully smaller or larger than loss. In this phase,
$\varrho_{\mathrm{st}}(\beta)$ is supported on two disconnected regions, with
the density mostly to the left of the origin. The middle plot shows a two-
component phase with large balanced loss and gain. In this regime,
$\varrho_{\mathrm{st}}(\beta)$ is supported on three disconnected regions,
with most measure centered at zero. In the rightmost plot, $m$ is large and
$\varrho_{H}(\epsilon)$ is approaching the Ginibre disk limit. The stationary
state distribution is nearly a Wigner semicircle with a positive finite mean,
corresponding to the fact that $\gamma>\tilde{\gamma}$.

###  3.4 Lindbladian Gas

This section discusses a Fermi gas in $d$ dimensions subject to loss and gain
of particles through Markovian exchange with a thermal bath. This provides a
simple example of a theory on a spatial continuum. Like the preceding section,
here a $U(1)$ symmetry is imposed to avoid complexity from the Nambu space.
For simplicity also, no additional matrix structure due to spin, orbitals,
flavor, etc. is considered. Note that even though a fermionic gas is
considered here, because of the symmetry, differences between the bosonic and
fermionic Nambu spaces do not enter and so the details are essentially the
same for the Lindbladian Bose gas.

In terms of the fermionic creation/annihilation operators
$\hat{c}(\mathbf{r})$ and $\hat{c}^{\dagger}(\mathbf{r})$, the many-body
Hamiltonian can be expressed in terms of the single-particle Hamiltonian
$H_{0}(\mathbf{r},\mathbf{r^{\prime}})$ as:

$\hat{\mathcal{H}}=\int\\!\text{d}\mathbf{r}\,\text{d}\mathbf{r^{\prime}}\,\hat{c}^{\dagger}(\mathbf{r})H_{0}(\mathbf{r},\mathbf{r^{\prime}})\hat{c}(\mathbf{r^{\prime}}).$
(91)

There are two families of jump operators:

$\hat{\mathcal{L}}_{1v}(\mathbf{r})=\int\text{d}\mathbf{r^{\prime}}\,\mu_{v}(\mathbf{r},\mathbf{r^{\prime}})\hat{c}(\mathbf{r^{\prime}}),\quad\hat{\mathcal{L}}_{2v}(\mathbf{r})=\int\text{d}\mathbf{r^{\prime}}\,\nu_{v}(\mathbf{r},\mathbf{r^{\prime}})\hat{c}^{\dagger}(\mathbf{r^{\prime}}),$
(92)

in terms of which the single-particle dissipation matrices are given by:

$Q(\mathbf{r},\mathbf{r^{\prime}})=\frac{1}{2}\sum_{v}\int\text{d}\mathbf{r_{1}}\text{d}\mathbf{r_{2}}\Big{(}\mu_{v}^{*}(\mathbf{r},\mathbf{r_{1}})\mu_{v}(\mathbf{r^{\prime}},\mathbf{r_{2}})+\nu_{v}(\mathbf{r},\mathbf{r_{1}})\nu_{v}^{*}(\mathbf{r^{\prime}},\mathbf{r_{2}})\Big{)},$
(93a)
$D(\mathbf{r},\mathbf{r^{\prime}})=\sum_{v}\int\text{d}\mathbf{r_{1}}\text{d}\mathbf{r_{2}}\Big{(}\mu_{v}^{*}(\mathbf{r},\mathbf{r_{1}})\mu_{v}(\mathbf{r^{\prime}},\mathbf{r_{2}})-\nu_{v}(\mathbf{r},\mathbf{r_{1}})\nu_{v}^{*}(\mathbf{r^{\prime}},\mathbf{r_{2}})\Big{)}.$
(93b)

The action for the Lindbladian Bose gas is given by:

$S=\int\text{d}t\,\text{d}\mathbf{r}\,\text{d}\mathbf{r^{\prime}}\begin{bmatrix}\bar{\psi}^{1}&\bar{\psi}^{2}\end{bmatrix}_{\mathbf{r}}\begin{bmatrix}i\partial_{t}-H_{0}+iQ&iD\\\
0&i\partial_{t}-H_{0}-iQ\end{bmatrix}_{\mathbf{r},\mathbf{r^{\prime}}}\begin{bmatrix}\psi^{1}\\\
\psi^{2}\end{bmatrix}_{\mathbf{r^{\prime}}}.$ (94)

In the simplest situation, the single-particle matrices are local differential
operators, so that
$H_{0}(\mathbf{r},\mathbf{r^{\prime}})=\delta(\mathbf{r}-\mathbf{r^{\prime}})H_{0}(\mathbf{r},-i\partial_{\mathbf{r}})$
and similarly for the dissipative matrices. When this is the case, the action
is local in spacetime,

$S=\int\text{d}x\begin{bmatrix}\bar{\psi}^{1}(x)&\bar{\psi}^{2}(x)\end{bmatrix}\begin{bmatrix}i\partial_{t}-H_{0}+iQ&iD\\\
0&i\partial_{t}-H_{0}-iQ\end{bmatrix}\begin{bmatrix}\psi^{1}(x)\\\
\psi^{2}(x)\end{bmatrix}.$ (95)

where $x=(t,\mathbf{r})$ denotes a spacetime coordinate.

One obtains a semi-classical kinetic equation by Wigner transform and
truncation of the gradient expansion of matrix products. Upon Wigner
transform, the parameter matrices become functions on the single-particle
phase space,

$H_{0}(\mathbf{r},\mathbf{k})=\int\text{d}\mathbf{r^{\prime}}e^{-i\mathbf{k}\mathbf{r^{\prime}}}H_{0}\bigg{(}\mathbf{r}+\frac{\mathbf{r^{\prime}}}{2},\mathbf{r}-\frac{\mathbf{r^{\prime}}}{2}\bigg{)},$
(96)

and similarly for $Q(\mathbf{r},\mathbf{k})$ and $D(\mathbf{r},\mathbf{k})$.
In the limit of slow variations, it is often appropriate to truncate the
Wigner expansion to first order in gradients. This amounts to a semiclassical
treatment of the problem and is known as the Wigner approximation. In this
approximation, kinetic equation becomes a Boltzmann equation with a linear
collision integral,

$\partial_{t}F-\partial_{\mathbf{r}}H_{0}\partial_{\mathbf{k}}F+\partial_{\mathbf{k}}H_{0}\partial_{\mathbf{r}}F=-2QF+D.$
(97)

Note that in a bosonic system, this kinetic equation will be of exactly the
same form. Going beyond this approximation by incorporating higher-order
terms, one finds the leading quantum corrections to the collision integral,

$-\frac{1}{2}\partial_{\mathbf{r_{1}}}\partial_{\mathbf{r_{2}}}Q\partial_{\mathbf{k_{1}}}\partial_{\mathbf{k_{2}}}F-\frac{1}{2}\partial_{\mathbf{k_{1}}}\partial_{\mathbf{k_{2}}}Q\partial_{\mathbf{r_{1}}}\partial_{\mathbf{r_{2}}}F+\partial_{\mathbf{r_{1}}}\partial_{\mathbf{k_{2}}}Q\partial_{\mathbf{k_{1}}}\partial_{\mathbf{r_{2}}}F.$
(98)

The inclusion of these terms brings the Boltzmann equation to the form of a
Fokker-Plank equation. This contrasts to the purely coherent situation, in
which there are no leading quantum corrections to the Boltzmann equation in
the absence of multiple bands.

For concreteness, specialize to a local single-particle Hamiltonian that is a
generalized Schrödinger operator,
$H_{0}(\mathbf{r},-i\partial_{\mathbf{r}})=\varepsilon(-i\partial_{\mathbf{r}})+V(\mathbf{r})$.
Then the phase space function
$H_{0}(\mathbf{r},\mathbf{k})=\varepsilon(\mathbf{k})+V(\mathbf{r})$ will be
of the form of a dispersion relation plus a potential. In addition, suppose
that the coupling to bath is translationally invariant, so that the
dissipation matrices are only functions of $\mathbf{k}$ in phase space. Then
the kinetic equation is identifiable as a Boltzmann equation in the relaxation
time approximation,

$\big{(}\partial_{t}+\mathbf{v}_{\mathbf{k}}\partial_{\mathbf{r}}-(\partial_{\mathbf{r}}V)\partial_{\mathbf{k}}\big{)}F=\frac{F_{0}-F}{\tau},$
(99)

where $\mathbf{v}_{\mathbf{k}}=\partial_{\mathbf{k}}\varepsilon$ is the group
velocity and $\tau(\mathbf{k})=1/2Q(\mathbf{k})$ is the relaxation time. The
distribution $F_{0}(\mathbf{k})=D(\mathbf{k})/2Q(\mathbf{k})$ takes the place
of the equilibrium distribution; in the absence of an external potential
$V=0$, one finds $F_{\mathrm{st}}(\mathbf{k})=F_{0}(\mathbf{k})$.

###  3.5 Mean Field Theory

This section illustrates how the above formalism can be extended to treat non-
linear systems. Using a mean-field approach, the semiclassical kinetics of the
previous section can be extended to self-consistently accommodate
interactions. For specificity, bosonic systems will be focused on, though as
with the previous section the details are similar for the fermionic analogue.

Like the previous section, no additional matrix structure is considered beyond
the spatial degrees of freedom. Only terms respecting the $U(1)$ particle
number symmetry are considered, so that no Nambu space structure has to be
dealt with. Additionally, the Hamiltonian and jump operators are assumed to be
invariant under spatial translations and rotations. The non-interacting
Hamiltonian and jump operators are given by the bosonic versions of eq.s (91)
and (92). Because of translational symmetry, in the Wigner representation
$H_{0},\ Q$ and $D$ are only functions of $\mathbf{k}$,

$Q(\mathbf{k})=\frac{1}{2}\sum_{v}\Big{(}\mu_{v}^{*}(\mathbf{k})\mu_{v}(\mathbf{k})-\nu_{v}(\mathbf{k})\nu^{*}_{v}(\mathbf{k})\Big{)},$
(100a)
$D(\mathbf{k})=\sum_{v}\Big{(}\mu_{v}^{*}(\mathbf{k})\mu_{v}(\mathbf{k})+\nu_{v}(\mathbf{k})\nu^{*}_{v}(\mathbf{k})\Big{)}.$
(100b)

Non-linearity can occur either on the level of the Hamiltonian or in the jump
operators. For simplicity, consider a contact interaction,

$\hat{\mathcal{H}}_{\mathrm{int}}=\frac{U}{2}\int\text{d}\mathbf{r}\,\hat{a}^{\dagger}(\mathbf{r})\hat{a}^{\dagger}(\mathbf{r})\hat{a}(\mathbf{r})\hat{a}(\mathbf{r}).$
(101)

For the jump operators, consider a local two-body loss and gain,

$\hat{\mathcal{L}}_{3}(\mathbf{r})=\frac{M}{\sqrt{2}}\,\hat{a}(\mathbf{r})\hat{a}(\mathbf{r}),\qquad\hat{\mathcal{L}}_{4}(\mathbf{r})=\frac{N}{\sqrt{2}}\,\hat{a}^{\dagger}(\mathbf{r})\hat{a}^{\dagger}(\mathbf{r}).$
(102)

Note that the corresponding Lindbladian defined this way is similar to various
models for driven-dissipative condensates [3, 29, 30, 31]. Here the nonlinear
interactions are consider only as a weak perturbation to a stable linear
theory. A condensate occurs when the theory is unstable on the linear level
and is not treated here.

The Keldysh action has the form $S=S_{0}+S_{\mathrm{int}}$, where quadratic
part of the action $S_{0}$ is given by the bosonic version of eq. (95),

$S_{0}=\int\text{d}x\begin{bmatrix}\bar{\phi}^{\mathrm{c}}(x)&\bar{\phi}^{\mathrm{q}}(x)\end{bmatrix}\begin{bmatrix}0&i\partial_{t}-H_{0}-iQ\\\
i\partial_{t}-H_{0}+iQ&iD\end{bmatrix}\begin{bmatrix}\phi^{\mathrm{c}}(x)\\\
\phi^{\mathrm{q}}(x)\end{bmatrix}.$ (103)

The non-linear part of the action is:

$S_{\mathrm{int}}=-\frac{1}{2}\int\text{d}x\Big{(}U(\bar{\phi}^{\mathrm{c}}\phi^{\mathrm{c}}+\bar{\phi}^{\mathrm{q}}\phi^{\mathrm{q}})(\bar{\phi}^{\mathrm{q}}\phi^{\mathrm{c}}+\bar{\phi}^{\mathrm{c}}\phi^{\mathrm{q}})\\\
-iJ(\bar{\phi}^{\mathrm{c}}\phi^{\mathrm{c}}-\bar{\phi}^{\mathrm{q}}\phi^{\mathrm{q}})(\bar{\phi}^{\mathrm{q}}\phi^{\mathrm{c}}-\bar{\phi}^{\mathrm{c}}\phi^{\mathrm{q}})-iR\bar{\phi}^{\mathrm{c}}\phi^{\mathrm{c}}\bar{\phi}^{\mathrm{q}}\phi^{\mathrm{q}}\Big{)},$
(104)

where $J$ and $R$ are defined similarly to $Q$ and $D$,

$J=\frac{1}{2}\big{(}|N|^{2}-|M|^{2}\big{)},\qquad
R=2\big{(}|N|^{2}+|M|^{2}\big{)}$ (105)

Diagrammatically, these interactions are four-point vertices as represented in
Fig. (10).

Figure 10: The left two figures depict the vertices depicting the interactions
in the nonlinear action eq. (104). The leftmost is comparable to one of the
vertices that appears in the equilibrium theory, but come with additional
imaginary factor. The $R$ vertex has no equilibrium analogue and is purely
dissipative. The left two figures are diagrammatic representations of the
mean-field contributions to the action from the interactions. The collective
field $\varphi$ enters as a fully renormalized bubble formed from two
classical legs. The left and right respectively renormalize the single-
particle dispersion and distribution. Note that the latter is uniquely a
feature of the Lindbladian dynamics; in equilibrium theory there is no
modification to the Keldysh component of eq. (104) on the mean-field level.

While in principle the full range of diagrammatic techniques can be applied to
study this model, here a mean-field treatment is discussed as an extension of
the quadratic formalism. To achieve this, one should replace factors of
$\bar{\phi}^{\mathrm{c}}\phi^{\mathrm{c}}$ in the action with their
expectation value. This reduces the nonlinearities to a quadratic coupling to
the collective field
$\varphi(x)=\frac{1}{2}\braket{\phi^{\mathrm{c}}(x)\bar{\phi}^{\mathrm{c}}(x)}$
whose value can be determined self-consistently from this definition. To be
specific, all three of the parameter matrices are modified,

$\delta H_{0}(x)=U\varphi(x),\qquad\delta Q(x)=J\varphi(x),\qquad\delta
D(x)=R\varphi(x).$ (106)

Following section 2.4, one can seek a solution to the now time-dependent
kinetic equation eq. (35). In the Wigner approximation, this will be of the
same form as eq. (99),

$\big{(}\partial_{t}+\mathbf{v}_{\mathbf{k}}\partial_{\mathbf{r}}-(\partial_{\mathbf{r}}V)\partial_{\mathbf{k}}\big{)}F=\frac{F_{0}-F}{\tau},$
(107)

with $\mathbf{v}_{\mathbf{k}}=\partial_{\mathbf{k}}H_{0}(\mathbf{k})$,
$V(x)=U\varphi(x)$, $\tau^{-1}=2Q(\mathbf{k})+2J\varphi(x)$ and
$F_{0}/\tau=D(\mathbf{k})+R\varphi(x)$.

Solutions to the kinetic equation determine $F$ as a function of $\varphi$.
This in turn can be fed into the definition of $\varphi$ to determine its
value self-consistently. To be precise, in the Wigner approximation one can
express the spectral Green’s function as:

$G^{\mathrm{R}}(x,p)\simeq\frac{1}{\epsilon-
H_{0}(\mathbf{k})+iQ(\mathbf{k})-(U-iJ)\varphi(x)},$ (108)

where $p=(\epsilon,\mathbf{k})$. Comparably, the Keldysh Green’s function at
equal spacetime points is given approximately by:

$iG^{\mathrm{K}}(x,x)\simeq
i\int\frac{\text{d}\epsilon}{2\pi}\frac{\text{d}\mathbf{k}}{(2\pi)^{d}}F(x,\mathbf{k})\big{(}G^{\mathrm{R}}(x,p)-G^{\mathrm{A}}(x,p)\big{)}.$
(109)

The integral over $\epsilon$ is just the frequency integral of the spectral
function and is thus equal to $1$. Note that unlike in the equilibrium theory,
there is no assumption that the spectral function is sharply peaked. A quasi-
particle approximation is generally not valid, as for general Lindbladian
systems due to the presence of dissipative terms which are not generically
small. Despite this, the frequency dependence of the distribution function may
anyways be dropped in this mean-field treatment due to the Markovian nature of
the dynamics. This gives the self-consistency condition for the collective
field $\varphi(x)=\frac{i}{2}G^{\mathrm{K}}(x,x)$,

$\varphi(x)=\frac{1}{2}\int\frac{\text{d}\mathbf{k}}{(2\pi)^{d}}F(x,\mathbf{k}).$
(110)

This together with the Boltzmann equation eq. (107) constitute a closed system
of two equations for $F$ and $\varphi$.

In the stationary limit, $F_{\mathrm{st}}$ and $\varphi_{\mathrm{st}}$ are
fully isotropic. The stationary solution can be read off from the Boltzmann
equation as $F_{0}$. Together with the self-consistency condition, this gives
the set of equations:

$F_{\mathrm{st}}(\mathbf{k})=\frac{D(\mathbf{k})+R\varphi_{\mathrm{st}}}{2Q(\mathbf{k})+2J\varphi_{\mathrm{st}}},\qquad\varphi_{\mathrm{st}}=\frac{1}{2}\int\frac{\text{d}\mathbf{k}}{(2\pi)^{d}}F_{\mathrm{st}}(\mathbf{k}).$
(111)

Slow relaxation can be studied by linearizing the kinetic equation around this
stationary solution. Writing $F=F_{\mathrm{st}}+\delta F$ and
$\varphi=\varphi_{\mathrm{st}}+\delta\varphi$, one has the closed system of
equations:

$\big{(}\partial_{t}+\mathbf{v}_{\mathbf{k}}\partial_{\mathbf{r}}+2Q+2J\varphi_{\mathrm{st}}\big{)}\delta
F=\partial_{\mathbf{r}}\delta\varphi\partial_{\mathbf{k}}F_{\mathrm{st}}+(R-2JF_{\mathrm{st}})\delta\varphi,$
(112a)
$\delta\varphi(x)=\frac{1}{2}\int\frac{\text{d}\mathbf{k}}{(2\pi)^{d}}\delta
F(x,\mathbf{k}).$ (112b)

By Fourier transforming in the spacetime coordinate $x=(t,\mathbf{r})$ to
$(\omega,\mathbf{q})$, one may algebraically solve for $\delta F$. In doing
so, the self-consistency condition gives the condition for a non-trivial
solution $\varphi$,

$1+\frac{1}{2}\int\frac{\text{d}\mathbf{k}}{(2\pi)^{d}}\frac{U\mathbf{q}\partial_{\mathbf{k}}F_{\mathrm{st}}-i(R-2JF_{\mathrm{st}})}{\omega-\mathbf{q}\mathbf{v}_{\mathbf{k}}+2i(Q+J\varphi_{\mathrm{st}})}=0.$
(113)

This relation fixes $\omega$ as a function of $\mathbf{q}$, which specifies
the dispersion of collective modes which govern the relaxation the density
field $\varphi$.

In the absence of dissipation $Q=0=D$ and $J=0=R$, this equation determines
the dispersion of a coherent sound mode $\omega(\mathbf{q})\simeq
c|\mathbf{q}|$ at small momenta, where the speed of sound $c$ is determined
from the relation

$1+\frac{U}{2}\int\frac{\text{d}\mathbf{k}}{(2\pi)^{d}}\frac{\mathbf{q}\partial_{\mathbf{k}}F_{\mathrm{st}}}{c|\mathbf{q}|-\mathbf{q}\mathbf{v}_{\mathbf{k}}}=0,$
(114)

where in the equilibrium theory $F_{\mathrm{st}}$ is the equilibrium quasi-
particle distribution function. This is the familiar equilibrium zero-sound
mode. To examine the how this is modified due to the presence of dissipation,
consider first the simpler situation of a weak purely linear dissipation that
is independent of $\mathbf{k}$, so that $R=0=J$ and $Q(\mathbf{k})=1/2\tau$.
Then for small $\mathbf{q}$ one finds,

$\omega(\mathbf{q})\simeq-i/\tau+c|\mathbf{q}|,$ (115)

where $c$ is the speed of sound determined from eq. (114) with
$F_{\mathrm{st}}(\mathbf{k})=\tau D(\mathbf{k})$. Thus, for momenta small
compared to the inverse relaxation time, $|\mathbf{q}|\ll 1/c\tau$, the sound
mode becomes over-damped.

In the more generic setting with non-linear dissipation, it is difficult to
make general statements without a specific form of the single-particle
dispersion. However, one can see that the above behaviour is generic for small
momenta. The zero momentum limit of eq. (113) gives the relation:

$1+\frac{1}{2}\int\frac{\text{d}\mathbf{k}}{(2\pi)^{d}}\frac{2JF_{\mathrm{st}}-R}{2(Q+J\varphi_{\mathrm{st}})-\Gamma}=0,$
(116)

where $\Gamma=i\omega(\mathbf{q}=0)$, demonstrated by this equation to be
purely real. Thus, while the specific form of the collective mode dispersion
$\omega(\mathbf{q})$ for finite $\mathbf{q}$ depends on the microscopic
details of the single-particle dispersion, it is always over-damped for
sufficiently small momenta.

##  4 Conclusion

We have presented a tutorial treatment of a many-body Lindbladian dynamics of
driven-dissipative systems. We have employed the functional formalism, which
naturally follows from the generic closed time contour formalism, under the
assumption of Markovian (i.e. time-local) bath correlators. As demonstrated,
it allows one to evaluate local observables, various correlation functions,
linear response characteristics, and collective modes spectra.

One of the major goals of this review is to emphasize the existence of two
distinct quantities, characterizing dynamics of these non-equilibrium models:
the complex effective Hamiltonian, $\check{H}$, and the stationary
distribution function $\check{F}_{\mathrm{st}}$. The complex effective
Hamiltonian, determines the transient relaxation spectrum as well as the
linear response to certain perturbations. On the other hand,
$\check{F}_{\mathrm{st}}$ dictates steady-state observables, shows up in
spectra of collective modes, and participates in some linear responses. While
in equilibrium the two are rigidly related through the fluctuation-dissipation
theorem, they are essentially independent within the Lindbladian dynamics
framework. Moreover, as is repeatedly demonstrated above, they exhibit
qualitatively different properties. For example, complex spectra of the
effective Hamiltonian generically exhibit exceptional points, where two or
more eigenvalues collide. The relaxation characteristics feature non-analytic
behavior in the vicinity of such exceptional points in the parameter space.
Yet, $\check{F}_{\mathrm{st}}$ and the long time stationary properties are
completely smooth. We provide other examples illustrating qualitative
differences between the two quantities in the non-equilibrium setting.

While spectra (and to some extent eigenfunctions) of the complex Hamiltonian
received some attention, the stationary distribution went largely unexplored.
We have shown here that it is determined by the effective kinetic equation. In
the particular case of linear systems, such kinetic theory acquires the form
of the so-called Lyapunov equation of the matrix algebra. Although there
aren’t many standard analytic tools to deal with it, it may be treated with
stable and efficient numerical algorithms.

The tools, outlined here, allow one to completely solve quadratic many-body
Lindbladian problems by diagonalizing $N\times N$ complex Hamiltonian and
solving $N\times N$ Lyapunov equation for the stationary distribution. Notice
that the dimensionality of the corresponding fermionic Hilbert space is
$2^{N}$. Therefore one achieves the exponential reduction in the problem’s
complexity. An immediate extension of the quadratic theory is the mean-field
approximation, which deals with the linearized treatment near a certain (self-
consistent) state.

Still a lot to be done yet for better understanding truly non-linear many-body
Lindbladian dynamics. In our opinion, the techniques presented here are
indispensable for this goal. One of the most exciting applications of the
functional methods is in the study of non-perturbative (instanton) effects
[99], which provide, eg., an ultimate floor for the qubit decoherence rate.
Other examples of essentially non-linear phenomena include studies of non-
equilibrium phase transitions [100, 101, 102], and various applications of the
functional renormalization group to driven-dissipative problems [3, 31, 103,
104, 105].

##  Appendix A Keldysh-Nambu Diagonalization

This appendix examines the classical mechanics of the Keldysh action. One can
diagonalize the quadratic form by means of a complex coordinate
transformation. This is a generalized form of Bogoliubov rotation, extended to
the Keldysh phase space space. In this new set of coordinates, the dynamics
can be understood through semiclassical quantization.

For bosons, one must perform a complex canonical transformation on the Keldysh
phase space. To this end, one should write the Keldysh action in the form of a
Hamiltonian-Lagrangian. Ignoring the constant term, the action from eq. (11)
is:

$S=\int\text{d}t\Big{(}\Phi^{\mathrm{q}}i\partial_{t}\Phi^{\mathrm{c}}-\frac{1}{2}\Phi^{\alpha}\check{K}_{\alpha\beta}\Phi^{\beta}\Big{)},$
(117)

where here $\Phi^{\mathrm{q}}=[\bar{\phi}^{\mathrm{q}}\ -\phi^{\mathrm{q}}]$
is defined differently than in the main text. The classical field
$\Phi^{\mathrm{q}}$ plays the role of the canonical position and the quantum
field $\Phi^{\mathrm{q}}$, its conjugate momentum. The matrix
$\check{K}_{\alpha\beta}$ is given by:

$\check{K}=\begin{bmatrix}0&\check{H}^{\mathrm{T}}\\\
\check{H}&-i\check{\tau}^{3}\check{D}\check{\tau}^{3}\end{bmatrix}.$ (118)

Note that this is a symmetric matrix on the full $4N\times 4N$ Keldysh-Nambu
space.

The classical mechanics of the Keldysh action is a generalized Hamiltonian
mechanics on the $4N$-dimensional Keldysh phase space. The equations of motion
are given by Hamilton’s equations

$\partial_{t}\Phi^{\alpha}=-[\mathcal{K},\Phi^{\alpha}]_{\mathrm{PB}},$ (119)

where the Keldysh Poisson bracket is defined by:

$[\cdot,\cdot]_{\mathrm{PB}}=\check{\sigma}^{2}_{\alpha\beta}\reflectbox{$\vec{\reflectbox{$\partial$}}$}_{\Phi^{\alpha}}\vec{\partial}_{\Phi^{\beta}},$
(120a)
$[\Phi^{\mathrm{c}}_{s},\Phi^{\mathrm{q}}_{s^{\prime}}]_{\mathrm{PB}}=-i\delta_{ss^{\prime}}.$
(120b)

The matrix $\check{J}=i\check{\sigma}^{2}$ defines the symplectic form of the
Keldysh classical mechanics.

With the correct choice of coordinates, one can express $\mathcal{K}$ in a
diagonal form. Because $\check{K}$ is a symmetric matrix, $\check{J}\check{K}$
is an element of the complex symplectic algebra $\mathfrak{sp}(4N)$, defined
through the relation
$\check{J}(\check{J}\check{K})+(\check{J}\check{K})^{\mathrm{T}}\check{J}=0$.
It can be brought to Jordan canonical form by a (generically complex)
symplectic matrix $\check{V}\in\mathrm{Sp}(4N,\mathbb{C})$,

$\check{V}\check{J}\check{K}\check{V}^{-1}=\begin{bmatrix}\check{U}\check{H}\check{U}^{-1}&0\\\
0&-(\check{U}\check{H}\check{U}^{-1})^{\mathrm{T}}\end{bmatrix}.$ (121)

As a symplectic matrix, $\check{V}$ obeys
$\check{V}^{\mathrm{T}}\check{J}\check{V}=\check{J}$. The $2N\times 2N$ blocks
of this matrix $\check{V}_{\alpha\beta}$ are given by:

$\check{V}_{\mathrm{c}\mathrm{c}}=\check{U},\quad\check{V}_{\mathrm{c}\mathrm{q}}=-\check{F}_{\mathrm{st}}\check{\tau}^{1}\check{U}^{\mathrm{T}-1},\quad\check{V}_{\mathrm{q}\mathrm{c}}=0,\quad\check{V}_{\mathrm{q}\mathrm{q}}=\check{U}^{\mathrm{T}-1},$
(122)

where $\check{V}$ obeys the symplectic condition
$\check{V}^{\mathrm{T}}\check{J}\check{V}=\check{J}$. This matrix defines a
complex canonical transformation to a new set of coordinates,

$\begin{bmatrix}\zeta\\\
\bar{\zeta}\end{bmatrix}=\check{V}\begin{bmatrix}\Phi^{\mathrm{c}}\\\
\Phi^{\mathrm{q}}\end{bmatrix},$ (123)

Note that each of the new $2N$ fields not related by complex conjugation,
$\bar{\zeta}_{s}\neq\zeta_{s}^{*}$. They are however symplectic conjugates,
obeying the relation:

$[\zeta_{s},\bar{\zeta}_{s^{\prime}}]_{\mathrm{PB}}=\delta_{ss^{\prime}}.$
(124)

In these coordinates, the action is brought to a canonical form in terms of
decoupled fields. In the diagonalizable case, this is:

$S=\sum_{s}\int\text{d}t\bar{\zeta}_{s}(i\partial_{t}-\epsilon_{s})\zeta_{s}.$
(125)

Written in this form, one can see that there are $2N$ integrals of the
classical motion, $I_{s}=\bar{\zeta}_{s}\zeta_{s}$. They are related to the
Keldysh Hamiltonian through $\mathcal{K}=\sum_{s}\epsilon_{s}I_{s}$. Applying
the Bohr-Sommerfeld quantization rule, one puts $I_{s}=n_{s}$ with $n_{s}$ a
set of positive integers. The Lindbladian spectrum is given by the quantized
values taken by the Keldysh Hamiltonian,
$-i\mathcal{K}=-\sum_{s}n_{s}\epsilon_{s}$. In the non-diagonalizable, there
are less than $2N$ integrals $I_{s}$, corresponding to the degeneracy of
eigenvalues. In this situation, the action will have additional terms of the
form $\bar{\zeta}_{s}\zeta_{s+1}$ depending on the Jordan block structure of
$\check{H}$.

A similar argument can be made for fermions. Writing the fermion action in a
Hamiltonian form, one has:

$S=\int\text{d}t\Big{(}\Psi^{2}i\partial_{t}\Psi^{1}-\frac{1}{2}\Psi^{a}\check{K}_{ab}\Psi^{b}\Big{)},$
(126)

with $\Psi^{2}=[\bar{\psi}^{1}\ \psi^{2}]$ defined differently than in the
main text. The matrix $\check{K}$ is given by:

$\check{K}=\begin{bmatrix}0&-\check{H}^{\mathrm{T}}\\\
\check{H}&-i\check{D}\check{\tau}^{1}\end{bmatrix}.$ (127)

Note that this is an antisymmetric matrix on the Keldysh-Nambu space.

The pseudo-classical fermionic equations of motion can be defined in terms of
a fermionic Poisson bracket,

$\partial_{t}\Psi^{a}=-\\{\mathcal{K},\Psi^{a}\\}_{\mathrm{PB}}.$ (128)

The bracket is symmetric and defined through the Grassmann derivatives:

$\\{\cdot,\cdot\\}_{\mathrm{PB}}=-i\check{\sigma}^{1}_{ab}\reflectbox{$\vec{\reflectbox{$\partial$}}$}_{\Psi^{a}}\vec{\partial}_{\Psi^{b}},$
(129a) $\\{\Psi^{1}_{s},\Psi^{2}_{s^{\prime}}\\}=-i\delta_{ss^{\prime}}.$
(129b)

The matrix $\check{\sigma}^{1}$ defines the inner product on the Keldysh
Grassmann algebra. As such, the product $\check{\sigma}^{1}\check{K}$ is an
element of the complex orthogonal algebra $\mathfrak{so}(4N)$, obeying the
relation
$\check{\sigma}^{1}(\check{\sigma}^{1}\check{K})+(\check{\sigma}^{1}\check{K})^{\mathrm{T}}\check{\sigma}^{1}=0$.
It can be brought to Jordan canonical form by a complex orthogonal
transformation $\check{V}\in\mathrm{SO}(4N,\mathbb{C})$, which preserves the
inner product
$\check{V}^{\mathrm{T}}\check{\sigma}^{1}\check{V}=\check{\sigma}^{1}$. The
block components of $\check{V}$ are the same as in eq. (123), by replacing
index names $\mathrm{c}\to 1$ and $\mathrm{q}\to 2$. The new set of Grassmann
fields are given by:

$\begin{bmatrix}\xi\\\
\bar{\xi}\end{bmatrix}=\check{V}\begin{bmatrix}\Psi^{1}\\\
\Psi^{2}\end{bmatrix}.$ (130)

The semiclassical quantization of the moments $I_{s}=\bar{\xi}_{s}\xi_{s}$
gives the same result as for bosons, with the occupation numbers $n_{s}$
restricted to $0$ or $1$.

The usefulness of this representation beyond quadratic theory is limited.
Because the diagonal basis incorporates the stationary distribution matrix
$\check{F}_{\mathrm{st}}$ into its definition, it is not obvious how to use
this representation in an interacting theory or in the presence of time-
dependent perturbations. Conventionally $\check{F}$ is defined self-
consistently on the quadratic level to derive a renormalized quantum kinetic
equation that is non-linear in $\check{F}$. It is unclear how this procedure
could be performed, if at all, with $\check{F}$ absorbed into the definition
of the fields.

##  Appendix B Superoperator Quantization

In this appendix, the connection between the formalism presented here and the
third quantization superoperator formalism of [48, 49] is established. The
quantization of the Keldysh action reproduces the superoperator formulation of
the Lindbladian dynamics. The connection between the two approaches is
analogous to the relation between the Hilbert space and path integral
formulations of standard quantum theory, but has an added wrinkle: due to the
complex nature of the Keldysh classical mechanics, the quantization must be
“non-canonical” in that the generators of the boson and fermion superoperator
algebras are not related by adjoint.

For the bosonic theory, one can introduce the right/left boson superoperators
defined by:

$\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{a}}}_{+j}\rho=\hat{a}_{j}\rho;\quad\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{a}}}_{-j}\rho=\rho\hat{a}_{j}.$
(131)

These superoperators are bosonic in that they obey the bosonic commutation
algebra:

$[\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{a}}}_{+j},\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{a}}}_{+i}^{\dagger}]=\delta_{ij}=[\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{a}}}_{-j}^{\dagger},\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{a}}}_{-i}],\quad[\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{a}}}_{+j},\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{a}}}_{-i}]=0.$
(132)

With these, the Keldysh action is quantized by replacing the Keldysh fields
with the Keldysh-rotated superoperators,
$\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{a}}}_{\mathrm{c},\mathrm{q}j}=(\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{a}}}_{+j}\pm\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{a}}}_{-j})/\sqrt{2}$,
and the Keldysh Poisson structure from Appendix A eq. (120) with the quantum
commutator,

$[\cdot,\cdot]_{\mathrm{PB}}\to-i[\cdot,\cdot],$ (133a)
$\phi^{\alpha}\to\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{a}}}_{\alpha};\quad\bar{\phi}^{\alpha}\to\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{a}}}_{\alpha}^{\dagger}.$
(133b)

Using this procedure, the Keldysh Hamiltonian $-i\mathcal{K}$ of the action
eq. (11) is translated back into the Lindbladian defined by eq. (10).

The quantum and classical superoperators define the bosonic algebra:

$[\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{a}}}_{\mathrm{c}j},\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{a}}}^{\dagger}_{\mathrm{q}i}]=\delta_{ij},\quad[\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{a}}}_{\mathrm{c}j},\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{a}}}^{\dagger}_{\mathrm{c}i}]=0=[\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{a}}}_{\mathrm{c}j},\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{a}}}_{\mathrm{q}i}],$
(134)

The superoperator algebra in this basis reflects the complex nature of the
Keldysh classical mechanics. The canonical conjugate and complex conjugate
fields are note equivalent; the canonical conjugate of the field
$\phi^{\mathrm{c}}$ is $\bar{\phi}^{\mathrm{q}}$, not
$\bar{\phi}^{\mathrm{c}}$. As a consequence, the bosonic superoperator algebra
is non-conjugate. The classical superoperators
$\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{a}}}_{\mathrm{c}j}$ and
$\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{a}}}^{\dagger}_{\mathrm{c}j}$
act as creation operators, while the quantum superoperators
$\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{a}}}^{\dagger}_{\mathrm{q}j}$
and $-\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{a}}}_{\mathrm{q}j}$
are their corresponding annihilation operators respectively, despite not being
their adjoints.

The Lindbladian can be brought to Jordan canonical form by a generalized
Bogoliubov transformation of the bosonic superoperator algebra. This is
achieved by a change of basis implemented by $\check{V}$ from eq. (123). The
resulting set of bosonic superoperators
$(\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{b}}},\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{b}}}^{\prime})$
obey a non-conjugate version of the bosonic commutation algebra given by the
quantization of eq. (124),

$[\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{b}}}_{s},\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{b}}}^{\prime}_{s^{\prime}}]=\delta_{ss^{\prime}},\quad[\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{b}}}_{s},\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{b}}}_{s^{\prime}}]=0=[\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{b}}}^{\prime}_{s},\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{b}}}^{\prime}_{s^{\prime}}],$
(135)

where the prime denotes the creation superoperator that in general is not
related to the annihilation superoperator through adjoint,
$\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{b}}}_{s}^{\dagger}\neq\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{b}}}_{s}^{\prime}$.
A diagonalizable Lindbladian written in this basis adopts the familiar form of
a sum over number operators:

$\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{\mathcal{L}}}}=-i\sum_{s}\epsilon_{s}\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{b}}}_{s}^{\prime}\hat{\vphantom{\rule{1.0pt}{5.71527pt}}\smash{\hat{b}}}_{s}.$
(136)