Linear and non-linear response of quadratic Lindbladians

Quantum systems rarely exist in perfect isolation from their environment. While commonly regarded as detrimental for quantum information processing¹, the dynamics of open quantum systems have recently attracted substantial interest as a route towards harnessing engineered dissipation as a resource for creating exotic quantum states and information processing^2,3. Tailored environmental couplings in materials and quantum devices have long been appreciated as a possible ingredient for entanglement generation and quantum computation^4,5. At the same time, the steady states of open systems can offer rich new prospects for stabilizing unconventional topological^6,7 and many-body states^8,9,10 not readily found in thermal equilibrium, while posing new challenges in their theoretical and computational modeling¹¹. To understand signatures of such dissipative states of matter in experiment, a key question hence concerns how to model and interpret dynamical and spectroscopic responses in open quantum systems.

While open quantum systems are generically non-Markovian^12,13,14 and can be simulated effectively in only a limited number of cases¹¹, Markovian master equations and in particular the Lindblad master equation are prominent for their interpretibility, effectiveness as an approximation for steady-state behavior, and tractability for analytic work^15,16,17,18. Fundamentally, the Lindblad equation describes the time-local evolution of the density matrix of a small electronic or bosonic subsystem coupled to a larger Markovian reservoir, and has a natural field-theoretic representation via time-local self energies on the Keldysh contour^{19,20,21,22,23}. For fermions and bosons, the Lindblad master equation moreover can be usefully represented in terms of “third-quantized” superoperators which obey standard (anti)commutation relations^24,25,26. Remarkably, when such a Lindblad master equation is assembled from a non-interacting (quadratic) Hamiltonian and dissipative particle gain and loss processes, its eigenmodes can be computed exactly from diagonalizing a single particle matrix, which can then be used to construct the full many-body Lindbladian^24,25,26. This quadratic structure means that the Keldysh path integrals are Gaussian and hence Wick’s theorem can be used to reduce multi-point functions to two point functions as will be key in this work.

While Lindbladians originally originated from the desire to study individual quantum-optical emitters¹⁷, recently extended dissipative many-particle systems have attracted substantial attention. Particular model systems include boundary driven spin systems and quantum circuits that exhibit non-equilibrium phase transitions^{27,28,29,30,31}, and qubit systems where Lindbladians have been used to model dissipation³² and realize dissipative versions of the Ising model on quantum computers³³. Lindbladians have also been applied to reveal unusual dynamics in SYK systems^34,35, to realize dark states and dissipation-induced flat bands³⁶, study localization^37,38, and to stabilize entangled many-body states^39,40,41. In quantum circuits, Lindbladians have been used to complete logical operations⁴², and even as part of a feedback mechanism to mitigate errors in quantum computing experiments⁴³. Previous work on Lindbladian systems has focused on the structures of subspaces and symmetries^{32,36,44,45,46} as well as the response of these subsystems⁴⁷. Additionally, recent efforts have succeeded in classifying topological states in these systems⁴⁸ and in non-Hermitian systems more broadly⁴⁹. Building on these rapid advances of new and emergent behaviors in dissipative quantum systems, a natural question to address is how to elucidate and interpret their properties using experimentally-accessible spectroscopic tools.

In this paper, we develop a theory of linear and non-linear spectroscopy of the steady-state properties of Lindbladian quantum systems. For quadratic Lindbladians, we show that linear and non-linear spectroscopic responses have a succinct representation in terms of the biorthogonal eigenmodes of the third-quantized quadratic single-particle Lindbladian, and illustrate consequences of optical conductivities, second harmonic generation and shift currents of electronic Lindbladians, as well as dynamical spin correlation functions of boundary-driven spin chains. The paper is organized as follows: In the Lindblad Keldysh Green’s Function Formalism section, we start by reviewing Lindblad-Keldysh Green’s functions and third quantization, and compute the single-particle Green’s functions for quadratic Lindbladians. In the Spectroscopic Response Formalism section, we derive expressions for the linear and non-linear response of fermionic and bosonic open quantum systems governed by a quadratic Lindbladian; finally, in the Applications section, we apply our formalism to calculate response properties for (1) a 1D XY spin chain, (2) Bernal bilayer graphene, and (3) a bosonic optical lattice.

Lindblad Keldysh Green’s Function Formalism

We consider the time evolution of the density matrix ρ of a system of interest coupled to a reservoir (bath), generated by the Lindblad master equation (Liouvillian)

where the Lindblad superoperator is the generator of completely positive trace preserving (CPTP) time evolution and is given by^15,16

Here, H is the system Hamiltonian and J_m describe a set of quantum jump operators, weighted by a dissipation rate Γ (with rate differences between different jump operators absorbed in the definition of J_m). The dissipative component describes time evolution under a non-Hermitian Hamiltonian (ifrac{Gamma }{2}{sum }_{m}{J}_{m}^{dagger }{J}_{m}) subjected to quantum jumps to restore CPTP behavior—a key distinction from studies of classical non-Hermitian systems.

Quadratic Lindbladians comprise a rich class of open quantum systems for which the Hamiltonian is bilinear in fermions/bosons

and jump operators are linear in fermions/bosons

describing particle gain and loss processes. This quadratic structure ensures that Wick’s theorem can be used to obtain multi-point correlation functions in terms of two-point correlation functions. Additionally this quadratic structure leads to Gaussian path integrals with a single-particle matrix structure that can be leveraged to find these two-point functions: the single-particle Green’s functions.

The key object is the single-particle matrix (Xi ={mathcal{H}}+i{Sigma }^{R}), a combination of a Hermitian (coherent) term and an anti-Hermitian retarded self-energy contribution that describes Lindbladian dissipation. While Σ^R originates from gain and loss processes due to coupling with a bath, we note that similar dissipative processes can emerge in interacting closed systems. Lindbladians describe frequency-independent Σ^R, which naturally arise from Markovian baths with broad spectral functions. The formalism developed in this paper is therefore applicable for steady state or meta-stable responses in any regimes that permit approximating the Keldysh self energy via frequency-independent Lindbladian self energies.

Single-Particle Fermionic Lindbladians

To illustrate the analogy between quadratic Lindbladians and quadratic Hamiltonians, consider defining “left” and “right” superfermions^27,50, which act on the left side (({ell }_{alpha }rho ={c}_{alpha }rho {mathcal{P}})) and right side (({r}_{alpha }rho =rho {c}_{alpha }^{dagger }{mathcal{P}})) of the density matrix, with ({mathcal{P}}) the total fermionic parity operator. Here, α includes orbital/spin indices and continuous parameters other than the wave-vector. It is convenient to define these superoperators as ordinary operators acting on the vectorized density matrix, yielding

For finite-dimensional vector spaces, where matrix representations of the second quantized operator exist, the Lindblad master equation can be re-expressed in vectorized form as

where (overrightarrow{rho }) is a vector rather than a matrix and⁵¹

which decomposes in terms of jump operators J_L/R as³⁶

Here the operators are

Expanding the products of jump operators above now defines three single-particle matrices that comprise the dissipative action

where a_k,m = (a_k,m,1, …, a_k,m,N) and b_k,m = (b_k,m,1, …, b_k,m,N) are vectors of the jump operator coefficients and the matrices are constructed from the outer product of these coefficients. To make progress from here, we express (hat{{mathcal{L}}}) in terms of a generalization of Prosen’s “third-quantization” algebra^24,52. Defining ℓ = (ℓ₁, ℓ₂, …, ℓ_N), and r = (r₁, r₂, …, r_N) for N orbitals, the Lindbladian superoperator can be written as

with “third quantized” superfermions/superbosons ({Psi }_{{boldsymbol{k}}}^{dagger }=[{{boldsymbol{ell}}}_{{boldsymbol{k}}}^{dagger },{{boldsymbol{r}}}_{{boldsymbol{k}}}^{dagger},{{boldsymbol{ell}}}_{-{boldsymbol{k}}},{{boldsymbol{r}}}_{-{boldsymbol{k}}}])^24,36,52. The single-particle matrices ({L}_{{boldsymbol{k}}}^{{rm{coh}}}) and ({L}_{{boldsymbol{k}}}^{{rm{dis}}}) are

and

where the top and bottom of ± and ∓ refer to fermions and bosons respectively.

Here we note that the key to writing the Keldysh functional integral in the vectorized space is to reexpress operators in terms of their eigenvalues upon the application of complete sets of coherent states as detailed in^19,20; in this formalism the parity operators are naturally absorbed into the eigenvalues.

Single-Particle Fermionic Green’s Functions

To compute Green’s functions, we will be interested in the generating functional¹⁹

where vectors of left and right-contour fermion/boson fields are given by ℓ = (ℓ₁, ℓ₂, …, ℓ_N), and r = (r₁, r₂, …, r_N) for N orbitals. Here we use ℓ^†, r^† to indicate the Grassman conjugate field of ℓ, r. The left fermions correspond to operators acting on the left/contour moving to the left, and the right fermions correspond to operators acting on the right/contour moving to the right as illustrated in Fig. 1. The continuous-time Lindblad-Keldysh action S is formally given by¹⁹

where k is the wave-vector.

For quadratic Lindbladians, the expression for (hat{{mathcal{L}}}) given in Eq. (11) can be replaced by the quadratic superoperator of Eq. (19) written using left/right propagating fields and one obtains

Now, the physical responses are more transparent in a rotated basis obtained using the Keldysh–Larkin–Ovchinnikov transformation^53,54. For fermions the rotation matrix is

so that (tilde{{mathbf{Psi}}}=({{boldsymbol{psi}}}_{1},{{boldsymbol{psi}}}_{2},{bar{{boldsymbol{psi}}}}_{1},{bar{{boldsymbol{psi}}}}_{2}):= {U}_{L}({boldsymbol{ell}},{boldsymbol{r}},{{boldsymbol{ell}}}^{dagger},{{boldsymbol{r}}}^{dagger})={U}_{L}{mathbf{Psi }}). Note that ψ_i and ({bar{{boldsymbol{psi }}}}_{i}) are independent fields and while they are adjoints for bosons, they are unrelated for fermions²⁰. In this new basis,

for transformed single-particle matrices³⁶

and

These last expressions have recently also been obtained by a field-theoretic treatment conducted by Thompson and Kamenev in ref. ²¹, and are equivalent to the BdG form we obtained using third quantization in ref. ³⁶ and the form obtained by McDonald and Clerk in²². In this rotated basis one arrives at the inverse of the single particle Keldysh Green’s functions²⁰

for Lindbladians, which define spectroscopic response properties discussed below and take the usual Keldysh block structure

The Green’s functions can be succinctly written by introducing an effective non-Hermitian Hamiltonian

with Σ^R the retarded self energy, and the Lindblad-Keldysh self energy

Inverting Eq. (29), one obtains retarded and advanced Green’s functions

where ({G}^{A}={({G}^{R})}^{dagger }), since H = H^†, A = A^†, and B = B^†. The retarded/advanced Green’s functions encode the effective non-Hermitian single-particle spectrum that dictates spectroscopic responses. Conversely, the Keldysh Green’s function can be obtained via

which is true, only in the case of the steady state. For transient behavior Σ^K takes a different form and G^K = G^RF−FG^A for some Hermitian matrix F that specified the distribution function^20,21. Analogously to in closed systems one may consider the Fourier transform of these expressions to frequency space where i∂_t ↦ ω and we do not need to include an infinitesimal iη since there is naturally a finite self energy in these open systems. Doing so we obtain

which can be decomposed in terms of the right (leftvert urightrangle) and left (leftlangle bar{u}rightvert) eigenstates of Ξ as

and

Quadratic Bosonic Lindbladians

For bosonic systems, the generating functional Eq. (22) and action Eqs. (23)–(24) are the exact same and analogous arguments define “left” and “right” superbosons

where we note that Eq. (11) has the same form except that ({J}_{m,R}:= {sum }_{alpha }{a}_{m,alpha }{r}_{alpha }+{b}_{m,alpha }{r}_{alpha }^{dagger }) for bosons. The A, B, and C matrices have the same form as in the fermionic case and are given by Eq. (16).

Now we can perform the Keldysh–Larkin–Ovchinnikov rotation⁵³, where for bosons the rotation matrix is

so that (tilde{{mathbf{Psi }}}=({{boldsymbol{psi }}}_{1},{{boldsymbol{psi }}}_{2},{bar{{boldsymbol{psi }}}}_{1},{bar{{boldsymbol{psi }}}}_{2}):= {U}_{L}({boldsymbol{ell}},{boldsymbol{r}},{{boldsymbol{ell}}}^{dagger },{{boldsymbol{r}}}^{dagger })={U}_{L}{mathbf{Psi }}).

Expressed in the new basis, we have

and

which will be key to obtaining the Green’s functions and calculating response properties.

For bosons, the Keldysh block Green’s function takes the block structure

Introducing again an effective non-Hermitian Hamiltonian (Xi ={mathcal{H}}+i{Sigma }^{R}) for bosons

where contrary to the fermionic case the imaginary part of Σ^R can be negative corresponding to possible gain processes that are impossible in the fermionic case. This means that contrary to a “fully filled” steady state in fermionic systems there are bosonic systems that have physical choices for the dissipative jump operators where there is no steady state due to unbounded gain. Generically, the precise condition for there to be no steady state is that σ³Ξ (or if C = Δ = 0 then (h-ifrac{Gamma }{2}(A-B))) acquires an eigenvalue with positive imaginary part: in this case there will be a mode that grows exponentially with time. Expressions for densities and correlation functions below assume that all poles are in the negative half-plane. The retarded and advanced Green’s functions are

with ({G}^{A}={({G}^{R})}^{dagger }). The Pauli σ³ matrix acting in Nambu space (and as the identity in orbital space) arises since (bar{psi }) and ψ fields are related for bosons (but are independent for fermions). The Keldysh self energy term reads

and the Keldysh component of the Green’s function again becomes G^K = G^RΣ^KG^A for the steady state response.

Spectral Representation and Non-Hermitian Single-Particle Band Structure

Importantly, spectroscopic responses for quadratic Lindbladians generalize the response formalism of band insulators and metals in a two-fold manner: first, the effective single-particle band structure generalizes to a non-Hermitian problem, with the disparity of left and right eigenvectors having important consequences on, e.g., optical transition matrix elements. Second, the steady-state distribution of dissipative systems is generically distinct from thermal distributions in equilibrium.

To discuss ramifications of the non-Hermitian spectrum that enters into the definition of retarded and advanced Green’s functions, we start by describing its eigenmodes. We have right and left eigenvectors

where ({sigma }^{i}={mathbb{1}}) for fermions and σⁱ = σ³ for bosons which acts on the Nambu structure as σ³ and as the identity on the orbital structure. Here, (leftlangle {u}_{n}rightvert ={leftvert {u}_{n}rightrangle }^{dagger }ne leftlangle {bar{u}}_{n}rightvert), and

Left and right eigenvectors are biorthogonal for all k

and completeness (away from exceptional points)⁵⁵ dictates

This permits a spectral representation of the steady-state single particle Green’s functions

Finally, for the non-equilibrium steady-state, the Keldysh Green’s function can be written as

Green’s Functions

We will be interested in evaluating two-point correlation functions which follow naturally by considering the lesser and greater Green’s functions

where ϕ = π for fermions and 2π for bosons is the exchange angle. Expressed in terms of G^< and G^> the retarded, advanced, and Keldysh Green’s functions become⁵³

This naturally leads to the identity

The two point correlation functions read

where we used that 1/e^iϕ = e^iϕ for fermions and bosons, and reindexed in the first line. For quadratic Lindbladians, multi-point correlation functions can now be computed straightforwardly using Wick’s theorem.

Spectroscopic Response Formalism

Armed with a representation of Lindblad-Keldysh Green’s functions in terms of effective non-Hermitian spectra, we now transcribe linear and non-linear dynamical responses for dissipative system into the framework of quadratic Lindbladians. Detailed derivations of these frequency dependent response functions is given in the Supplemental Material (See Supplemental Material, Section 2 and 3 for derivations of the frequency dependent response functions.).

Spectral Function

To set the stage, we start by considering the single-particle spectral function that describes physical single-particle excitations and their density of states. The spectral function is given by

where ({sigma }^{i}={mathbb{1}}) for fermions and σⁱ = σ³ which acts on the Nambu structure for bosons. This can then be used to obtain the density of states

and the number of bands is then given by

where N_k is the same for all k.

Particle Density and Equal-Time Expectation Values

To obtain the steady-state particle density, consider the equal time expectation value

which follows from what we obtained above in Eq. (69). Using Eq. (61), one immediately finds the density using the frequency-space Keldysh Green’s function; we obtain

where we used the residue theorem to complete the ω integral and ({sigma }^{i}={mathbb{1}}) for fermions and σⁱ = σ³ which acts on the Nambu structure for bosons. The density is then the diagonal elements (langle {n}_{alpha }rangle =langle {c}_{alpha }^{dagger }{c}_{alpha }rangle).

By analogy, the steady state expectation value for an arbitrary single-particle operator

can be computed from the dissipative eigenbasis

where O_k is the 2N × 2N matrix with entries ({[{O}_{{boldsymbol{k}}}]}_{alpha beta }={O}_{{boldsymbol{k}},alpha beta }) and the trace runs over both state indices n and momentum k. Note that O_k is typically traceless in Nambu form, if time-reversal symmetry is respected.

Linear Response

We now turn to the steady state linear response of an arbitrary operator ({mathsf{O}}(t)) to a time-dependent perturbation ({{mathsf{O}}}^{{prime} }(t)), given by the correlation function

Generically, a Kubo formula can also be derived for dissipative systems. This is presented in the Supplemental Material (See Supplemental Material, Section 1 for a derivation of the Kubo formula for Lindbladian systems.) and remains generically applicable for interacting Lindbladians. For quadratic Lindbladians, Wick’s theorem permits expressing the correlation function in terms of single-particle Keldysh Green’s functions in frequency space

where the details are given in the Supplemental Material (See Supplemental Material, Section 2 for derivations of the frequency dependent linear response functions). Expressing in terms of a spectral representation and completing the integral over frequencies ω using the residue theorem we find

where ϕ = π and ({sigma }^{i}={mathbb{1}}) for fermions and ϕ = 2π and σⁱ = σ³ which acts on the Nambu structure for bosons. The response corresponds diagrammatically to two loops as illustrated in Fig. 2, and obeys standard Kramers-Kronig relations between real and imaginary parts.

Example: Optical Conductivity and f-Sum Rule

Consider now the steady state linear optical response of a dissipative electron systems. The linear optical conductivity in velocity gauge is given by

where Π_para is given by Eq. (81) with O_k = j^μ the paramagnetic current operator, and ({Pi }_{,text{dia},}^{mu nu }=-ilangle {j}^{mu nu }rangle) with j^μν the diamagnetic current operator.

In an open quantum system, a definition of a local current operator within the system is given by the time derivative of the electric dipole operator ex = ie ∇_k⁵⁶, which in a Lindbladian setting yields an equation of motion ({j}^{mu }=-iebar{{mathcal{L}}}{{x}^{mu }}) with (bar{{mathcal{L}}}) the adjoint Lindbladian for the Heisenberg equation of motion. We take e = 1 in the discussion that follows. For simplicity we will consider situations where the bath cannot carry a current, leading to momentum-independent jump operators. In this situation, the current operator recovers the usual closed-system form ({j}^{mu }=i[{mathcal{H}},{x}^{mu }]) with ({mathcal{H}}) the coherent system Hamiltonian.

In closed systems, the frequency-sum (f-sum) rule for optical responses⁵⁷ relates the integral over the paramagnetic (real part) of the optical response to the ground state expectation value of the diamagnetic current

a relation which is satisfied at each k-point individually and permits counting the number of carriers by measuring the frequency-dependent absorption. In quadratic Lindbladian systems we find that a f-sum rule holds for the k-integrated response for a bath that does not contribute to the current (with Lindbladian jump operators that are momentum-independent, i.e., a and b are momentum independent). Proving this relation and investigating other quantum geometric response properties of these systems, such as establishing a Thouless-Kohmoto-Nightingale-de Nijs (TKNN) formula⁵⁸ for open systems is an intriguing direction for future study.

Second-Order Non-Linear Response

For the second order response we have the correlation function

where the new term at this order is the triangle diagram illustrated in Fig. 3. In analogy to the linear response, we can use Wick’s theorem, express in terms of Green’s functions and Fourier transform to find

which can be diagrammatically expressed as in Fig. 3. As above, the trace includes an integral over k. Now, after substituting the expression for the spectral representation and completing the ω integral we find

where ϕ = π and ({sigma }^{i}={mathbb{1}}) for fermions and ϕ = 2π and σⁱ = σ³ which acts on the Nambu structure for bosons. This can now be immediately evaluated in terms of the normal modes of σⁱΞ. The derivation is given in Supplemental Material (See Supplemental Material, Section 3 for derivations of the frequency dependent nonlinear response function).

Example: Non-Linear Optical Conductivity

As an example we consider the non-linear response to light in velocity gauge, where as in the closed system the non-linear optical conductivity is given by⁵⁹

where

The terms including ({{mathsf{J}}}^{mu nu rho }) and ({{mathsf{J}}}^{mu nu }) can be evaluated using the linear response formalism of the previous subsection using Eq. (79) and Eq. (82) with ({j}^{mu nu rho }={partial }_{{k}_{mu }}{partial }_{{k}_{nu }}{partial }_{{k}_{rho }}{h}_{{boldsymbol{k}}}) and ({j}^{mu nu }={partial }_{{k}_{mu }}{partial }_{{k}_{nu }}{h}_{{boldsymbol{k}}}) respectively where ({mathsf{J}}) and j are related by Eq. (78). The final triangle term can be evaluated using Eq. (88) with (O={j}^{mu },{O}^{{prime} }={j}^{nu },{O}^{{primeprime} }={j}^{rho }).

Now, when the system is driven by one frequency of light there are two natural second order responses: second harmonic generation/frequency doubling, and shift current generation/frequency cancellation that are given by⁵⁹

and

which can be naturally evaluated using the formulae above.

Applications

Having established a response formalism for quadratic Lindbladians, we now apply it to diverse physical examples illustrating its applicability. In the first example we study an XY spin chain which through the Jordan-Wigner transformation can be expressed in terms of free fermions. This fermionic model has superconducting terms and leverages the full Nambu structure of the theory. In the second example, we consider the paradigmatic material Bernal bilayer graphene and its linear and non-linear optical responses. In the third example, we consider a bosonic optical lattice and consider the momentum space atomic occupation that results from realistic slight anisotropies in dissipation rates.

XY Spin Model

In equilibrium, the 1D J_x, J_y (XY) spin-1/2 model exhibits three phases: an oscillatory spin-density wave (SDW) phase at small fields and small J^x/J^y anisotropy, an ordered ferromagnetic (FM) phase at small fields and large anisotropy, and a disordered paramagnetic (PM) phase at large fields⁶⁰. Recently the Ising limit of this model has been realized on a quantum computer where simulated dissipation cooled the system towards its ground state³³. Remarkably, in the dissipatively boundary-driven XY model similar phases to equilibrium can emerge^25,27,28,31, however their signatures in physical spectroscopic responses and consequences on magnetic excitations are an open question. To address this, we recast the boundary driven model as a quadratic Lindbladian and compute the dynamical spin response at finite transverse field strengths. We find that, in contrast to the gapped modes of the equilibrium model, there are gapless modes whose coupling to the single site spin flip operator ({S}_{i}^{z}) decreases as the spin-density wavelength increases.

Model

Consider first the paradigmatic J_x, J_y spin chain model with N sites and open boundary conditions

Additionally we consider the effect of local potentials (transverse fields) h_n

so that we have H = H_XY + H_⊥. The standard Jordan-Wigner transformation recasts this model in terms of free fermions

with a single-particle Hamiltonian where we multiplied the J_x and J_y terms by 2 so that the critical transverse field in the Ising limit of the closed system is at h_⊥ = 1

Now, suppose that the model is subjected to boundary dissipators given by the jump operators

where the Jordan-Wigner string for J_N on the two contours cancels as the total parity of the superfermions is conserved—annihilation of one fermion on one contour corresponds to annihilation of one fermion on the other contour and so the total partiy is preserved. The sign is unimportant for the Lindblad evolution since each jump operator appears with its adjoint. These jump operators describe the coupling of the ends of the XY chain to two oppositely spin-polarized leads. These will preferentially (and incoherently) attempt to collapse the spins at the two ends of the chain to ↑ and ↓ respectively. The model is depicted schematically in Fig. 4(a) where a 1D chain is subjected to a transverse field and is sandwiched in between ferromagnetic reservoirs.

Magnetization Density

We are now in a position to calculate the steady-state z-polarization magnetization density, ({S}^{z}=n-frac{1}{2}), by using the fermionic response formalism given in Eq. (77) and above. We find a steady-state SDW phase with mean magnetization zero and long-range spin density wave order as illustrated in Fig. 4b, and a paramagnetic (PM) phase with mean magnetization zero and no long-range order: spin-spin correlations decay exponentially. The critical line between the SDW and the PM phases of the dissipative chain is at

Notably, the location of the critical line is independent of the dissipation rate Γ. While as illustrated in Fig. 4, Γ plays a key role in determining the amplitude of oscillations in the SDW phase, the oscillations only exactly vanish at and above the critical field value in the PM phase, independent of Γ. This critical field can be calculated from inspecting the bulk Jordan-Wigner spectrum of the XY spin chain Hamiltonian

For large transverse fields h in the PM phase, the Jordan-Wigner dispersion has a minimum at k = 0. Conversely, for small transverse fields h in the SDW phase, the dispersion exhibits two minimal at k = ±k_SDW with

The jump between k = 0 and k = ± k_SDW dispersion minima occurs at ({h}_{perp }^{c}({J}_{y})=4{J}_{x}{J}_{y}/({J}_{x}+{J}_{y})) and signifies the phase transition upon adding the boundary dissipation. As the phase boundary is crossed from PM to SDW, the dissipative steady state can now preferentially populate the ± k_SDW bulk modes, resulting in a steady state SDW order with momentum k_SDW. We note that this mechanism relies on an intrinsically non-equilibrium steady-state distribution, and stands in contrast to the ferromagnetic to paramagnetic transition in closed equilibrium systems which occurs when the spectrum becomes gapless, i.e., ϵ(0) = 0 at ({h}_{perp }^{c}({J}_{y})={J}_{y}+{J}_{x})⁶⁰.

We illustrate the phase diagram in Fig. 4c and note key features and symmetries of the phase diagram. First, at the Ising point J_x = 1, J_y = 0 there is no spin-density wave phase as there is only a non-magnetized phase. Second, at the finely tuned point J_y = J_x, the model is a critical XX spin chain and the SDW phase vanishes. Third, at J_y = − J_x, the SDW phase persists for all transverse fields: this can be understood through a seesaw mechanism where ferromagnetic correlations in one direction are precisely counterbalanced by correlations in the opposite direction leading to persistent oscillatory behavior, even at strong fields. Moreover the critical line of the model exhibits two symmetries: first, there is ({{mathbb{Z}}}_{2}) symmetry corresponding to spin-flip symmetry, and second there is a glide symmetry: upon reflecting the phase boundary lines across the seesaw point (purple line) in Fig. 4c they are the same as those on the other side of the seesaw point before reflection up to a shift in h_⊥ by ± 8J_x.

While the dissipation rate Γ does not affect the wavelength of the SDW oscillations, it has a drastic effect on the excitation spectrum, changing the amplitude of the oscillations as we illustrate in Fig. 4d, e. The amplitude of the SDW oscillations is maximized when Γ ~ J where J is the nearest neighbor spin-spin coupling. This can be understood using second order perturbation theory, where we consider the first spin in the system. This spin is coupled both to a ferromagnetic reservoir by Γ and to the rest of the system (more or less paramagnetic) by J. When Γ ≪ J, the coupling to the center of the system dominates and oscillations are weak. When Γ ≫ J there are no virtual transitions (via second order perturbation theory) from the first spin to the second spin and so oscillations are weak. Finally, when Γ ~ J, the amplitude of oscillations peaks since there is both a strong drive and a strong coupling to the next spin in the chain which means that the virtual transitions from perturbation theory are significant.

Near the critical point, the wavelength of the spin density wave diverges. We attribute a critical exponent to this divergence and find − 1/2 as expected for a free theory. This critical exponent matches the exponent for correlation functions defined by Eisert and Prosen in ref. ²⁸. We illustrate this power law divergence in Fig. 4f, with noise for long wavelengths arising when commensurability with the system size (N = 1000) becomes significant. Far from ({h}_{perp }^{c}), there are deviations from the power law behavior. The wavelengths λ are extracted by Fourier transforming the z-magnetization density of the central region of the spin chain and selecting the frequency of the largest peak.

Now, it is significant to note that while the wavelength of the SDW diverges at the critical point, its amplitude vanishes as we illustrate in Fig. 4g. Near the critical point for any l ≪ λ for spin-density wavelength λ we have

since the correlation function is over spins near each other on one wavelength and hence have asymptotically equal magnetizations. The average is over all sites n, and for numeric simulations we take l = 2. Studying this numerically, we find the amplitude vanishes with a critical exponent of roughly 5/2. The data is noisy over the whole range studied because the amplitude oscillates with the constructive and destructive interference of SDWs originating from the left and right ferromagnetic reservoirs. Again, far from the critical point we see deviations from power-law behavior.

Dynamic Spin Susceptibility

Above we considered the steady state magnetization, now we consider the dynamic response of the spin chain to a frequency dependent perturbation. This is analogous to calculations of time-dependent spin correlation functions in^61,62, except that in our framework frequency rather than time is the natural footing and so there is no need to Fourier transform over a finite time to obtain the dynamic spin susceptibility. Here we specialize to the case of an XY model with J_x = 1, J_y = 1/3 that exhibits both SDW and paramagnetic phases. As in the static case, the amplitudes peak when Γ~J, so we choose Γ = 1.

We proceed to calculate the dynamic spin-spin susceptibility using Eq. (81) for the one-loop finite frequency response. We consider the operators (O={S}_{i}^{z}) and ({O}^{{prime} }={S}_{j}^{z}) where i and j are site indices in the chain. We identify ({mathcal{S}}(i,j,Omega )={Pi }^{ij}(Omega )) as the spin-spin correlation function, and compute the specific correlations ({mathcal{S}}(2i,Omega )={mathcal{S}}(N/2-i,N/2+i,Omega )) which is the spin-spin correlation function centered around the middle of the chain, where i ranges from 0 to N/2 − ϵ where we choose N = 301 and ϵ = 25 to avoid finite size boundary effects. We then Fourier transform this to obtain ({mathcal{S}}(q,Omega )). Note that (langle [{S}_{j}^{z},{S}_{i}^{z}]rangle =langle [{n}_{j},{n}_{i}]rangle) so we are free to complete this computation using fermions.

Visualizing ({mathcal{S}}(q,Omega )) in Fig. 5, we see that there are gapless modes in contrast to the closed system where all modes are gapped for ({h}_{perp }le {h}_{perp }^{c}). Notably, while the real part of the spectrum of Ξ would suggest that the open system has a gap above the ({rm{Re}}(xi )=0) modes there are gapless excitations because the modes are partially filled as given by the distribution function. In fact, for our choice of jump operators all states are half-filled except four states corresponding to boundary modes acted on directly by the dissipators. In Fig. 5a(i) and Fig. 5c(i) a gap in the spectrum of ({mathcal{S}}(q,Omega )) appears between Ω = 1 and 1.2, corresponding to the gap in the spectrum of Ξ: while transitions are still allowed within the lower (and upper) manifold of states, the gap from the lower to the upper manifold exceeds the bandwidth of the lower (and upper) manifold and so no transitions are possible at these frequencies. In contrast in Fig. 5a(ii) and c(ii) no gap is visible since the bandwidth of the manifolds exceeds the gap between them.

In the weak and ultra-strong dissipation limits, Fig. 5a, c respectively, gapless modes disperse from (pm 2{q}_{{rm{SDW}}},{rm{mod}},2pi) for spin-density wavenumber q_SDW which exists for ({h}_{perp }le {h}_{perp }^{c}=1) here. Above ({h}_{perp }^{c}) closed systems similarly exhibit a gapless response. Additionally for ({h}_{perp }le {h}_{perp }^{c}) we note that there are excitations at q_SDW corresponding to exciting spin-density waves. The weak dissipation and ultra-strong dissipation cases have comparable spectra with a clear correspondence between the dominant modes because, as illustrated in Fig. 4e, for Γ not comparable to J coupling to the bulk is weak and the effect of dissipation is felt most strongly near the boundaries.

In the strong dissipation case, at commensurate Γ ~ J, inversion symmetry instead becomes strongly broken as the effects of dissipation persist deep in the bulk. At h_⊥ = J_x/3 and h_⊥ = 2J_x/3, as visualized in Fig. 5b(i–ii) the single spin flip operators (O={S}_{i}^{z}) and ({O}^{{prime} }={S}_{j}^{z}) couple strongly to a series of dispersive modes that are gapless. Meanwhile at larger field strengths, as illustrated in Fig. 5b(iii–iv), the operators couple more weakly to the gapless dispersive modes. This is also the regime in which the spin-density wave order of the steady state has a long wavelength. This may be because the ideal description of the excitations is as a more-extended object and the single-spin flip probe is too local to efficiently probe the excitations. Above ({h}_{perp }^{c}), Fig. 5b(iv), we see that at large energies there is only one mode with significant spectral weight, while in Fig. 5b(ii) there are many modes. We interpret this as modes that are scattered by multiples of q_SDW (and backfolded by 2π) in the spin-density wave phase, but that collapse onto a single mode above ({h}_{perp }^{c}) when the spin-density wave order vanishes and q_SDW goes to 0.

Bilayer Graphene

Over the past decade, twisted bilayer graphene has emerged as a platform of choice for realizing many unusual physical phenomena such as superconductivity^63,64, and many other correlated electronic behaviors⁶⁵ in flat bands near the Fermi energy^66,67. The simplest and most thermodynamically stable graphene bilayer is AB stacked “Bernal” bilayer. Recently, superconductivity⁶⁸ and the quantum anomalous Hall effect⁶⁹ have been observed in Bernal bilayer graphene. Here we take a new approach by strongly coupling bilayer graphene to metallic top and bottom leads and instead analyze its dissipative steady states due to potential differences between the leads. Importantly, dissipation will alter the single-particle electronic properties and lead to the formation of exceptional points. Given this fundamentally non-Hermitian phenomenon, an interesting question is how these properties are manifested in linear and non-linear optical response properties that can be probed in experiments. We illustrate the system geometry in Fig. 6a where the bottom layer is coupled to reservoir 1 and the top layer is coupled to reservoir 2.

Model

We start from a tight-binding model for the graphene bilayer; we consider the Bernal structure as illustrated in Fig. 6b. The hopping terms are t(r) = t_ppπ(r) + t_ppσ(r) where the terms are given by^66,70

which are the terms that emerge using the Slater-Koster method for the p_z orbitals. The energies that set the problem scale are V_ppπ = 2.7 eV and V_ppσ = −0.48 eV, and the key lengths are the intra-layer spacing a₀ = 0.142 nm, inter-layer spacing d = 0.335 nm, and exponential decay length δ = 0.0453 nm, finally e_z = (0, 0, 1) is the unit vector perpendicular to the bilayer. We implement a cutoff for the tight-binding hopping so that t(r) = 0 when ∣r∣ > 1 nm.

Next, we introduce spatially (sublattice) homogeneous, but layer dependent dissipation given by

where the first layer is subjected to particle gain with amplitude Γ(1 + γ) and the second layer is subjected to particle loss with amplitude Γ(1 − γ) as illustrated in Fig. 6(a). Γ sets the overall scale of dissipation, and γ sets an anisotropy in coupling between the lower layer to the lower reservoir and the upper layer to the upper reservoir. This corresponds to a bilayer sandwiched between two top and bottom reservoirs at different chemical potentials, where we neglect electric field effects induced by charge distributions in the reservoirs.

To emphasize the effect of dissipation, we will set Γ to 1 eV, to compete with energy scales of the (coherent) closed-system band structure; while perhaps exceedingly large for Bernal bilayers, similar physics occurs in large twist angle commensurate bilayer graphene but at a greatly reduced energy scales^71,72,73,74. This reduced energy scale means that in such large twist angle samples the relevant low energy system physics may occur at the same energy as a reduced system-reservoir coupling Γ.

Now, we observe some interesting features of the dissipation acting on the bilayer as exhibited in Fig. 6c. First, for zero anisotropy γ, dissipation acts solely to induce an overall broadening for electronic states. As a finite anistropy γ between tunneling into the high-bias and low-bias dissipative leads is induced, a non-equilibrium steady state forms, in which particles tunnel from top to bottom gate via the Bernal stack. Here, exceptional points form near the band maxima (where the bands are almost degenerate), and then migrate towards the K point as the anisotropy is increased. We can understand this by considering the bands away from the touching to be living in each layer and then hybridizing as they approach the crossing at K. Fig. 6d shows the steady state band occupation, which varies from completely empty at γ = − 1 corresponding to all loss processes (i.e., the high-bias gate is cut off from the sample) to completely filled at γ = 1 corresponding to all gain processes. In between these we see that the layer resolved occupation is quite anisotropic, where even at γ = 0 almost all the occupation is in layer 1 where the tunneling-in processes occur and layer 2 only becomes significantly occupied near γ = 1. We note that the difference in layer occupation corresponds to the rate V_ppσ at which the layers are coupled where at large V_ppσ, the layers have the same occupation.

Absence of Exceptional Points in the Spectral Function

A defining characteristic of non-Hermitian systems is possibility for the formation of exceptional points where eigenvalues and vectors coalesce and no longer form a complete basis. While extensively discussed in classical systems such as photonic crystals⁴⁹, signatures of exceptional points in open quantum systems are less clear.

To shed light on their spectroscopic signatures, we start from discussing their impact on the single-particle spectral function

While the spectral function exhibits a non-analyticity in ω at the exceptional point in each eigenstate n that coalesces, the full spectral function A(ω) = ∑_mA_m remains a smooth function of ω. We illustrate the eigenstate-resolved non-analyticity in the inset of Fig. 6f, meanwhile in the main panel of Fig. 6f we plot the full spectral function. The reason this non-analyticity is not inherited by the full spectral function is a fine cancellation due to contour-reversal symmetry where while differences of ξ terms lead to non-analyticities, these non-analyticities cancel when n is summed over. This property of exceptional points showing up in state-resolved responses but vanishing when all states are summed over is inherited by the other response properties as well.

Diamagnetic Optical Conductivity

We now turn to the diamagnetic optical response, which takes the form

where for the Bernal bilayer the current operator is traceless so we drop the trace term and ({j}_{{boldsymbol{k}}}^{mu nu }={partial }_{{k}_{mu }}{partial }_{{k}_{nu }}{h}_{{boldsymbol{k}}}), and ϕ = π and ({sigma }^{i}={mathbb{1}}) for fermions.

Notably as illustrated in Fig. 6e, the diamagnetic optical conductivity decreases above γ~0.79 which corresponds to a decreasing diamagnetic optical conductivity with increasing filling. This is notable since in closed systems the diamagnetic optical conductivity always increases monotonically with filling. So the decrease at large γ = 1 is an indication of the open nature of this system.

We note that in this case the optical conductivity is purely real (as in closed systems) which corresponds to the fact that the current operators here are Hermitian. Hence since Σ^K is Hermitian, Eq. (82) is equal to its conjugate. This is not the case when the current operators are non-Hermitian and then the response is both reactive and dissipative as opposed to just being reactive as in closed systems.

Paramagnetic Linear Optical Conductivity

We now study the paramagnetic contribution to optical conductivity, which again has both a real and an imaginary contribution in dissipative quantum systems. In Fig. 6h we see that the current-current correlation is relatively constant over a range of frequencies and anisotropies. Around this region of constant response there is an onset at a frequency that seems to correspond roughly to the energy separation to the remote bands at the K point, and this onset goes to zero frequency as γ → 1. On the other hand, the response vanishes above a threshold corresponding to high energy bands being extremely short lived (large imaginary part/broadening). Finally, we notice that around γ = 1 and Ω = 0.5 eV there is a slight crest in the optical conductivity which we attribute to some transitions near the K point.

Non-Linear Optical Conductivity

Now, we can use our non-linear response formalism to calculate the second harmonic and shift responses. For clarity we only include the contributions from the triangle diagram and not from the other paramagnetic and diamagnetic loops. Since the closed system is centrosymmetric, the second order response vanishes, as seen in Fig. 6i–j at γ = 0. For clarity, we have limited the non-linear response calculation to the triangle diagram term from Eq. (86), since the other contributions are qualitatively similar to the linear response. As γ increases, the dissipative couplings to the top and bottom reservoirs break inversion symmetry, leading to a non-centrosymmetric dissipative steady state. Consequently, the non-linear response increases but doesn’t peak until γ~0.8−1.0 where the layers are both nearly filled.

In the second harmonic response, Fig. 6i, we see two peaks, one corresponding to “narrow” interband transitions between neighboring pairs of bands around the Γ and M points, and one corresponding to “wide” interband transitions around K. To understand this response it is important to note that the distribution function of the non-equilibrium steady state is not that of a Fermi-Dirac distribution function but rather each of the bands is partially filled to some non-negligible extent.

Meanwhile in the shift response, Fig. 6j we see the complementary behavior where there is no response until an energy threshold corresponding to cycling at the K point. The onset of the shift response is rapid and it reaches a large magnitude, perhaps suggesting that using dissipation to break centrosymmetry is a promising direction for realizing devices with non-linear optical responses such as the shift current that enables solar cells.

Bosonic Optical Lattice

Optical lattices are a tunable platform to realize exotic physics^75,76, and can be used to realize arbitrary lattices in two dimensions^77,78. While optical lattices may approximate closed fermionic and bosonic systems, they are suspended in a near vacuum and are susceptible to loss and gain processes making dissipation fundamentally relevant for a full description of their physics^79,80. Here we consider non-interacting bosons in an optical lattice subject to gain and loss processes.

Here we consider a tight-binding model on a square lattice with nearest-neighbor hopping t and a sublattice dependent mass term m

where σ_i are the Pauli matrices. We consider dissipation that acts uniformly on the sublattices with dissipators

where Γ are scalars describing the amplitude of the processes and A and B refer to the two sublattices. The anisotropy between gain and loss rates and between sublattices could correspond to the physical scenario of a stronger confining field on the top side and a weaker confining field on the bottom side and a bilayer with A sites on the top layer and B sites on the bottom layer. We illustrate such a process in Fig. 7a where a cloud of atoms is incident from the upper left and leaves to the lower right.

Using this model we can obtain the steady state density using the expressions from the Particle Density and Equal-Time Expectation Values section. We find that at large m the gap is large and the density of bosons is relatively uniform across the Brillouin zone. In contrast at small m, the density varies across the zone and varies rapidly near the Dirac point/band gap minima. We illustrate this in Fig. 7b–d where the z-coordinate corresponds to the real part of the energy ξ and the coloration corresponds to the density of bosons at that momentum in the optical lattice.

We now study the finite-frequency correlation functions of the bosonic optical lattice. Using Eq. (82) with O = n_i on orbital i = 0, 1 for the interorbital correlation function which we plot in Fig. 8, where we have taken the magnitude to show the scale of the response, but not its phase. In the large mass limit, Fig. 8d, the bands are essentially the orbitals and the inter-orbital correlation is small. Meanwhile in the small mass limit, Fig. 8a, the bands are strong admixtures of the orbitals and the inter-orbital correlation is large.

Quadratic Lindbladian systems present a generalization of non-interacting closed quantum systems to include dissipation, leading to dissipative steady states with new non-equilibrium features. Understanding their spectral properties and physical responses is essential to interpret realizations in quantum devices or engineered solid state systems. Here we systematically developed a Lindblad-Keldysh formalism to compute linear and non-linear dynamical response properties. In quadratic Lindbladian systems, we showed that the response probes the spectrum of an effective non-Hermitian Hamiltonian for second-quantized fermion or boson superoperators. The formalism naturally applies to spin, electronic, and bosonic systems and their frequency dependent dynamic response properties, with generalizations to weakly interacting dissipative quantum systems an intriguing direction for future work.

We first focus on the 1D XY spin model in a transverse field with boundary dissipation realized via jumps to two oppositely polarized ferromagnetic reservoirs, which exhibits a non-equilibrium quantum phase transition between spin-density wave and paramagnetic phases^27,28. In this setting, we studied the dynamical spin susceptibility at and away from criticality, which exhibit markedly distinct regimes of response as a function the boundary dissipation strengths. In contrast to the closed-system transverse field XY model, the dissipative steady state exhibits a series of gapless dispersive modes that originate from the underlying spin density wave pattern and disperse from (pm 2{q}_{{rm{SDW}}},{rm{mod}},2pi). Additionally, spin-density wave amplitude excitations are present at ± q_SDW. The weak and ultra-strong dissipation limits exhibit strikingly similar spectra with an emergent inversion symmetry in the bulk and a correspondence of the dominant modes. This behavior originates from the fact that in both cases the effects of dissipation are confined to the boundary of the system and the bulk is relatively unaffected. In contrast, at commmensurate dissipation strength Γ ~ J, inversion symmetry is strongly broken and the asymmetry of dissipation at the two boundaries remains pronounced deep within the bulk. Finally, the response to local spin flips is most prominent at small fields, with excitations becoming more extended when the spin-density wavelength is large. With the dissipative Ising model recently realized in a quantum device³³, measuring the unconventional excitations predicted by this work are a promising direction for fingerprinting out-of-equilibrium ordered phases.

Turning to fermionic dissipative systems by example of bilayer graphene coupled to two reservoirs, we found that the diamagnetic optical conductivity can decrease with increasing occupation in stark contrast to closed systems; additionally for non-Hermitian current operators the diamagnetic response is dissipative rather than purely reactive. Moreover, second-order non-linear responses that are forbidden in the closed system by centrosymmetry emerge naturally as consequences of anisotropic coupling to a reservoir that lead to a non-centrosymmetric steady state. While we studied the AB stacked Bernal bilayer, the results are most reasonably obtained in large-twist angle twisted bilayer graphene where the energy of the quadratic band touching is drastically reduced^72,73.

Finally, in a bosonic optical lattice we found that sublattice-dependent gain and loss processes, for instance from sublattices living in different layers with different confining potentials, can lead to a steady state with a momentum-dependent occupation. When a sublattice dependent mass term is applied, the density in each band becomes uniform across the zone, in analogy to the equilibrium expectation of uniform occupation. Meanwhile in the inter-orbital dynamic response, the orbital-resolved response is small when the mass term is large and bands are well-approximated by atomic orbitals; in contrast the orbital-resolved response is large when the bands are mixtures of atomic orbitals.

This work has a number of direct applications and possible generalizations. First, the computation of dynamical and spectroscopic response properties is essential to interpret dissipative steady states in quantum circuit realizations of quadratic Lindbladians on NISQ devices³³. The developed formalism permits calculating response properties of arbitrary fermionic and bosonic systems that are described by a quadratic Lindbladian, making it immediately generalizable to weakly-interacting systems via perturbative expansions in Coulomb scattering. With rich phenomenologies established in closed interacting electron systems such as correlated metals^81,82 it will be interesting to study the ramifications of dissipation on their responses in an open system setting. Furthermore, charge or magnetic ordering transitions in dissipative interacting electron or boson systems are now simulable in cold atomic gases or quantum devices, necessitating dynamical probes of their properties. It will be interesting to generalize the presented formalism to mean field descriptions of dissipative ordered steady states, with applications to interacting photons and Rydberg atom arrays^11,83,84. Complementarily, driven-dissipative Floquet systems are known to exhibit interesting steady states and phase transitions with no thermal and closed-system counterpart^85,86,87. Probing their electromagnetic response would serve as a useful diagnostic of engineering non-equilibrium response properties. Furthermore, with the role of quantum geometry and topology recently under much scrutiny in resonant responses of free electron systems^88,89,90,91 and low-frequency responses of strongly-interacting electron systems⁹², it will be interesting to understand analogous signatures in dissipative responses in a non-Hermitian Lindbladian setting. Here, an interesting question is how to devise experimentally accessible signatures of Lindbladian topology⁴⁸. Additionally, following this work, key areas of interest include signatures of non-Hermitian exceptional points in Lindbladian systems⁴⁹.