Space-time first-order correlations of an open Bose Hubbard model with
incoherent pump and loss

Martina Zündel1*, Leonardo Mazza 2, Léonie Canet1,3, Anna Minguzzi1

1 Univ. Grenoble Alpes, CNRS, LPMMC, 38000 Grenoble, France

2 Université Paris-Saclay, CNRS, LPTMS, 91405, Orsay, France

3 Institut Universitaire de France, 5 rue Descartes, 75005 Paris

<EMAIL_ADDRESS>
## Abstract

We investigate the correlation properties in the steady state of driven-
dissipative interacting bosonic systems in the quantum regime, as for example
non-linear photonic cavities. Specifically, we consider the Bose-Hubbard model
on a periodic chain and with spatially homogeneous one-body loss and pump
within the Markovian approximation. The steady state corresponds to an
infinite temperature state at finite chemical potential with diagonal spatial
correlations. Nonetheless, we observe a nontrivial behaviour of the space-time
two-point correlation function in the steady state, obtained by exact
diagonalisation. In particular, we find that the decay width of the propagator
is not only renormalised at increasing interactions, as it is the case of a
single non-linear resonator, but also at increasing hopping strength. We then
compute the full spectral function, finding that it contains both a dispersive
free-particle like dispersion at low energy and a doublon branch at energy
corresponding to the on-site interactions. We compare with the corresponding
calculation for the ground state of a closed quantum system and show that the
driven-dissipative nature – determining both the steady state and the
dynamical evolution – changes the low-lying part of the spectrum, where
noticeably, the dispersion is quadratic instead of linear at small
wavevectors.

###### Contents

  1. 1 Introduction
  2. 2 Model and observables
    1. 2.1 Lindblad equation
    2. 2.2 Correlation functions
    3. 2.3 Keldysh field theory
  3. 3 Exact properties of the non-equilibrium steady state
    1. 3.1 Steady state density matrix
    2. 3.2 Equal-time correlations in the steady state
  4. 4 Correlation function in real time
    1. 4.1 Analysis of the parameter space
    2. 4.2 Time correlations
    3. 4.3 Renormalisation of the temporal decay
  5. 5 The spectrum
    1. 5.1 Spectrum of the open system
    2. 5.2 Position of the excitation spectral lines for small hopping
    3. 5.3 Spectral function for the ground state of the Bose-Hubbard model
    4. 5.4 Comparison of unitary evolution and Lindblad evolution of the NESS
  6. 6 Conclusions and perspectives
  7. A Detailed calculations
    1. A.1 Bosonic NESS density matrix
    2. A.2 Approach to the NESS
    3. A.3 Equal-time correlation functions
    4. A.4 Oscillation frequency
  8. B Normalisation with truncation of the local Hilbert space
  9. C Steady state solution for fermions
    1. C.1 NESS density matrix
    2. C.2 Equal-time two-point correlation functions
  10. D Steady-state solution for hard-core bosons
  11. E Decay rate at large interactions
  12. F Details on the numerical calculations
    1. F.1 Exact diagonalisation
    2. F.2 Fit of the decay
    3. F.3 Effect of the Hilbert space truncation
    4. F.4 System size dependence
    5. F.5 Exact diagonalisation for the unitary system
    6. F.6 Simulations with matrix product states

##  1 Introduction

Closed quantum systems under a unitary time evolution have been largely
studied in the last two decades, with particular impetus after the
experimental reach of the strongly interacting regime both in three-
dimensional optical lattices, as well as in lower-dimensional setups [1, 2, 3,
4]. One-dimensional interacting quantum systems have attracted much attention
since a wealth of theoretical techniques are available, including in
particular exact solutions by Bethe Ansatz for homogeneous systems with
contact interactions for bosonic or fermionic gases, or mixtures thereof [5,
6, 7, 8, 9], and for trapped systems in the limit of infinitely repulsive
interactions [10, 11, 12, 13, 14]. While lattice fermions are integrable [15],
lattice bosons are not, except for the case of the hard-core bosons (see e.g.
[16, 17, 18]), and the Bose-Hubbard model for two particles [19, 20]. A
manifold of attempts to effectively describe non-integrable or ‘weakly’ non-
integrable systems have been put forward, such as the generalised Gibbs
ensemble [21], generalised hydrodynamics [22, 23], macroscopic fluctuation
theory [24], and these fields are under active investigation.

While all these studies concern closed quantum systems, in several
experimental situations quantum systems are subjected to external pump and/or
losses, or put in contact with some type of bath. Examples of bosonic open
systems are Josephson junction arrays [25], superconducting microwave circuits
[26], polariton chains [27], atoms or ions in optical cavities [28, 29],
optomechanical systems [30], and lossy quantum gases [31, 32, 33, 34, 35, 36,
37, 38, 39, 40, 41]. In contrast to the equilibrium case, much less is known
for interacting open quantum systems, and several questions arise, ranging
from the properties of the steady state and its phase diagram to its
coherence, its excitations and its dynamical properties. It is also very
interesting to investigate which properties known for unitary systems are
still present in open quantum systems, and under which conditions or
parameters.

A way to model an open system is to describe it as a subsystem embedded in a
larger system acting as a bath and with which it interacts. If one requires
the time evolution of the subsystem to be completely-positive and trace-
preserving (CPTP map), this evolution has to be of the Lindblad form [42], and
described by a Markovian quantum master equation. This evolution follows the
composition law for universal dynamical maps which is the quantum analogue to
the classical differential Chapman-Kolmogorov equation [43]. Despite its
simplicity, the Lindblad master equation suffices to display rich physical
behaviour, as it has been demonstrated for the case of the steady state of
interacting driven-dissipative bosons [44, 45]. Quadratic Lindbladians,
meaning Lindbladians with quadratic Hamiltonian and linear jump operators in
the bosonic (fermionic) creation and annihilation operators, are known to be
exactly solvable by the method of third quantisation, presented by Prosen [46]
and Prosen and Seligman [47]. Within the Keldysh formalism these systems
correspond to a quadratic action. Also instances of quasi-free Lindbladians
for fermions and bosons with quadratic hermitean jump operators have been
solved [48, 49]. Recently, the connection between both formalisms and the
phase-space formulation has been pioneered [50]. Non-quadratic Lindbladians
require in most cases a numerical solution. We refer to [51] for a comparative
analysis of the state-of-art numerical methods, specifically tailored for non-
equal time correlation functions. This method was used to obtain the second-
order correlation function $g^{(2)}$ [52].

An especially important quantity to characterise a quantum system is its two-
point correlation function, corresponding to first-order correlations
$g^{(1)}$, see e.g. [53]. At equal times, it contains information about the
spatial coherence and the presence of (quasi)off-diagonal long-range order,
i.e. the existence of Bose-Einstein condensation according to the Penrose
Onsager criterion [54]. Its Fourier transform yields the momentum distribution
function, giving information on the velocity of the quantum particles in a
given quantum state. The full space-time correlation function describes the
evolution of the quantum system when removing or adding a particle to it. In
space-time, it allows one to study the spread of correlations, and to test the
Lieb-Robinson bound [55], and its generalisation to open quantum systems [56].
Its Fourier transform in frequency-wavevector space provides the spectral
function, whose poles give the dispersion of single-particle excitations. The
spectral function is also needed to describe the transport properties among
two reservoirs across a quantum channel within linear response [57].

In this work, we study the two-point correlation function of an interacting
bosonic open quantum system described by the Bose-Hubbard model subjected to
incoherent drive and losses. At non-zero interactions, similarly to its
closed-system counterpart, this model is not integrable. We choose the case of
spatially uniform pump and losses, where the non-equilibrium steady-state
(NESS) density matrix is known exactly [58], and it is closely related to an
infinite-temperature state at fixed chemical potential, independently of the
interaction strength. Despite the steady state taking a simple form, we find
that the two-point correlation function shows a non-trivial dependence on the
interaction strength, displaying in particular typical temporal oscillations
with frequency related to the interaction strength as well as an interaction-
dependent exponential decay. We compare these results to the exactly known
non-interacting limit as well as to the strongly interacting limit where the
dynamics of the different sites decouples and can be described using the exact
results for the dissipative Kerr resonator [59]. We then proceed in obtaining
the full spectral function, where in particular we find a doublon-like branch
at excitation energy $U$ as well as a cosine-like excitation branch in the
single-particle-hole sector. The features of the spectral function are a clear
signature of the non-equilibrium nature of the quantum system, further
supported by a Keldysh action analysis, and by the comparison with the closed-
quantum system spectral function [60].

The paper is organised as follows. We start by presenting the model and the
characterisation of the non-equilibrium steady state in Sec. 2, which includes
the Lindblad equation in Sec. 2.1 and the Keldysh action in Sec. 2.3. We first
provide our predictions for the particle density as a function of time in Sec.
3. We then proceed to present our results for the space-time correlations in
Sec. 4, where we display the time-resolved correlations and analyse their
decay. The spectral function is shown in Sec. 5. Finally, Sec. 6 presents our
concluding remarks.

##  2 Model and observables

Figure 1: Sketch of the Bose-Hubbard model with periodic boundary conditions,
nearest-neighbour hopping amplitude $J$ and on-site interaction strength $U$,
coupled to spatially homogeneous Markovian baths with pump rate $\gamma_{p}$
and loss rate $\gamma_{l}$. The illustrated couplings apply to all sites.

###  2.1 Lindblad equation

We consider a lattice of $L$ sites occupied by bosonic particles with on-site
repulsive interactions, where each site is homogeneously coupled to a bath
allowing for one-body losses and pump. The corresponding unitary evolution is
described by the Bose-Hubbard Hamiltonian:

$\displaystyle
H=\sum_{i=1}^{L}\Big{[}-J\big{(}b_{i}^{\dagger}b_{i+1}+\text{h.c.}\big{)}+\frac{U}{2}b_{i}^{\dagger}b_{i}^{\dagger}b_{i}b_{i}\Big{]}\;,$
(1)

where $J$ is the hopping amplitude and $U$ the interaction strength. The
bosonic operators $b_{i}$ and $b_{i}^{\dagger}$ fulfill the usual commutation
relations $[b_{i},b_{j}^{\dagger}]=\delta_{ij}$,
$[b_{i},b_{j}]=0=[b_{i}^{\dagger},b_{j}^{\dagger}]$, and we take periodic
boundary conditions, i.e. we set $L+l\equiv l$.

The full temporal evolution is described by the Lindblad – Gorini-Kossakowski-
Sudarsahan equation for the system density matrix $\rho$ [42, 61]:

$\displaystyle\partial_{t}\rho=\mathcal{L}[\rho]=-\text{i}[H,\rho]+\gamma_{l}\sum_{i=1}^{L}\big{(}b_{i}\rho
b_{i}^{\dagger}-\frac{1}{2}\\{b_{i}^{\dagger}b_{i},\rho\\}\big{)}+\gamma_{p}\sum_{i=1}^{L}\big{(}b_{i}^{\dagger}\rho
b_{i}-\frac{1}{2}\\{b_{i}b_{i}^{\dagger},\rho\\}\big{)}\;.$ (2)

The master equation (2) can be obtained from a microscopic model where the
system is weakly coupled to a bath, under the assumption that time scales of
system and bath are well separated and that the correlations built in the bath
do not affect the system (see e.g. [43] for a detailed derivation).

###  2.2 Correlation functions

The two-point correlation function of interest in this work is the retarded
Green’s function

$\displaystyle
G^{R}_{j\ell}(t,t^{\prime})=-\text{i}\theta(t-t^{\prime})\Big{\langle}[b_{j}(t),b_{\ell}^{\dagger}(t^{\prime})]\Big{\rangle},$
(3)

where $\langle...\rangle$ indicates the trace weighted with the respective
density matrix, $\theta(t)$ is the Heavyside step function, and $b_{j}(t)$
($b^{\dagger}_{j}(t)$) is the bosonic annihilation (creation) field operator
in the Heisenberg picture, time-evolved according to Lindbladian dynamics111In
the case of open quantum systems, a way of defining the time evolution of an
operator $A$ is given by the adjoint superoperator $\bar{\mathcal{L}}$,
according to $A(t)=e^{\bar{\mathcal{L}}t}[A]$ where
$\bar{\mathcal{L}}[A]=+\text{i}[H,A]+\gamma_{l}\sum_{i=1}^{L}\big{(}b_{i}^{\dagger}Ab_{i}-\frac{1}{2}\\{b_{i}^{\dagger}b_{i},A\\}\big{)}+\gamma_{p}\sum_{i=1}^{L}\big{(}b_{i}Ab_{i}^{\dagger}-\frac{1}{2}\\{b_{i}b_{i}^{\dagger},A\\}\big{)}$.
Notice that ${\rm Tr}\big{\\{}A(t)\rho(0)\big{\\}}={\rm
Tr}\big{\\{}A\rho(t)\big{\\}}$. (in the open case), or Hamiltonian one (in the
closed case).

In order to obtain the retarded Green’s function for the non-equilibrium
steady state (NESS) of the driven-dissipative system, we evaluate unequal time
correlation functions, such as

$\displaystyle\Big{\langle}b_{j}^{\dagger}(t)b_{0}\Big{\rangle}_{\text{NESS}}^{(\mathcal{L})}\equiv{\rm
Tr}\Big{\\{}b_{j}^{\dagger}\text{e}^{\mathcal{L}t}\,\Big{[}b_{0}\rho_{\text{NESS}}\Big{]}\Big{\\}}\;.$
(4)

It will be useful to compare it with the correlation function obtained
following a unitary evolution starting from the same NESS density matrix:

$\displaystyle\Big{\langle}b_{j}^{\dagger}(t)b_{0}\Big{\rangle}_{\text{NESS}}^{(H)}\equiv\text{Tr}\Big{\\{}\text{e}^{\text{i}Ht}\,b_{j}^{\dagger}\text{e}^{-\text{i}Ht}\,b_{0}\rho_{\text{NESS}}\Big{\\}}\;.$
(5)

For a closed system the corresponding correlation function evaluated on the
ground state $|\Psi_{0}\rangle$ reads

$\displaystyle\Big{\langle}b_{j}^{\dagger}(t)b_{0}\Big{\rangle}_{\text{GS}}\equiv\langle\Psi_{0}|\text{e}^{\text{i}Ht}\,b_{j}^{\dagger}\text{e}^{-\text{i}Ht}\,b_{0}|\Psi_{0}\rangle.$
(6)

We finally consider the spectral function, whose poles provide the excitation
branches of the system, and which is defined as

$\displaystyle A(\omega,k)=-\frac{1}{\pi}{\rm Im}\,G^{R}(\omega,k)\;$ (7)

for a homogeneous system, and with the Fourier conventions

$\displaystyle f(p,\omega)$
$\displaystyle=\sum_{j\in\mathbb{Z}}\int_{-\infty}^{\infty}dt\,f_{j}(t)\,e^{ipj}e^{-i\omega
t};$ $\displaystyle f_{j}(t)$
$\displaystyle=\dfrac{1}{L}\int_{-\infty}^{\infty}\dfrac{dp}{2\pi}\int_{-\infty}^{\infty}\dfrac{d\omega}{2\pi}\,f(p,\omega)\,e^{-ipj}e^{i\omega
t}\,.$ (8)

###  2.3 Keldysh field theory

It is useful to formulate the model in terms of the Keldysh formalism,
following the derivation of [45]. Starting from the Lindblad equation (2), and
expanding the system density matrix onto two coherent-state bases for each
space-time point identified by the fields $\Psi_{j}^{+}(t)$
($\Psi_{j}^{-}(t)$) for upper (lower) time contours, one obtains the time
evolution of the density matrix in terms of a path integral. Following the
usual notations, the Keldysh action in the rotated basis
$({\Psi_{c}}_{j}(t),{\Psi_{q}}_{j}(t))=\frac{1}{\sqrt{2}}(\Psi_{j}^{+}(t)+\Psi_{j}^{-}(t),\Psi_{j}^{+}(t)-\Psi_{j}^{-}(t))$
reads

$\displaystyle S=\int\text{d}t\sum_{j}\Bigg{[}$
$\displaystyle{\Psi_{c}}_{j}^{*}\,\text{i}\left(\partial_{t}-\frac{\gamma_{l}-\gamma_{p}}{2}\right){\Psi_{q}}_{j}+{\Psi_{q}}_{j}^{*}\,\text{i}\left(\partial_{t}+\frac{\gamma_{l}-\gamma_{p}}{2}\right){\Psi_{c}}_{j}+\text{i}(\gamma_{l}+\gamma_{p})|{\Psi_{q}}_{j}|^{2}$
$\displaystyle+J\Big{(}{\Psi_{c}}_{j+1}^{*}{\Psi_{q}}_{j}+{\Psi_{q}}_{j+1}^{*}{\Psi_{c}}_{j}+{\Psi_{c}}_{j}^{*}{\Psi_{q}}_{j+1}+{\Psi_{q}}_{j}^{*}{\Psi_{c}}_{j+1}\Big{)}$
$\displaystyle-\frac{U}{2}\Big{(}({\Psi_{c}}_{j}^{2}+{\Psi_{q}}_{j}^{2}){\Psi_{c}}_{j}^{*}{\Psi_{q}}_{j}^{*}+({{\Psi_{c}}_{j}^{*}}^{2}+{{\Psi_{q}}_{j}^{*}}^{2}){\Psi_{c}}_{j}{\Psi_{q}}_{j}\Big{)}\Bigg{]}\;.$
(9)

In the same notation, the retarded Green’s function (3) reads

$\displaystyle
G^{R}_{j\ell}(t,t^{\prime})=-\text{i}\Big{\langle}\Psi_{cj}(t)\Psi_{q\ell}^{*}(t^{\prime})\Big{\rangle}_{c}\,.$
(10)

The saddle-point approximation in the path integral with the action (2.3)
allows one to study the existence and stability of a condensate. For a static
homogeneous field solution, if $\gamma_{p}<\gamma_{l}$, the imaginary part of
the Keldysh potential has a global minimum for $\Psi_{cj}=\Psi_{qj}=0$, which
implies that the system is in the symmetric phase without breaking of the U(1)
symmetry. Hence, there is no finite order parameter. As we shall see below,
this is in agreement with the exact NESS solution, which rather corresponds to
an infinite-temperature state for the bosons. Notice that this is at variance
with the case where both one- and two-body losses are present, where a U(1)
symmetry-broken phase emerges for $\gamma_{p}>\gamma_{l}$ and a
(quasi-)condensate forms at weak interactions [62].

Expanding the action around this solution to second order in the fluctuations,
we obtain the inverse of the retarded sector of the propagator as

$\displaystyle[G_{0}^{R}]^{-1}_{ml}(t^{\prime},t)$
$\displaystyle\equiv\begin{pmatrix}\frac{\delta^{2}S}{\delta{\Psi^{*}_{qm}}(t^{\prime}){\delta\Psi_{cl}}(t)}&&\frac{\delta^{2}S}{\delta{\Psi^{*}_{qm}}(t^{\prime}){\delta\Psi^{*}_{cl}}(t)}\\\
\frac{\delta^{2}S}{\delta{\Psi_{qm}}(t^{\prime}){\delta\Psi_{cl}}(t)}&&\frac{\delta^{2}S}{\delta{\Psi_{qm}}(t^{\prime}){\delta\Psi^{*}_{cl}}(t)}\end{pmatrix}\Bigg{|}_{\Psi_{c}=\Psi_{q}=0}$
$\displaystyle=\begin{pmatrix}(\text{i}\partial_{t}+\text{i}\kappa_{0})\delta_{m,l}+J(\delta_{m,l+1}+\delta_{m,l-1})&&0\\\
0&&-(\text{i}\partial_{t}+\text{i}\kappa_{0})\delta_{m,l}+J(\delta_{m,l+1}+\delta_{m,l-1})\end{pmatrix}$
$\displaystyle\times\delta(t-t^{\prime})\;,$ (11)

where we defined

$\kappa_{0}=(\gamma_{l}-\gamma_{p})/2,$ (12)

and where the 0 index indicates that this is the ‘bare’ propagator, obtained
within the quadratic (non-interacting) approximation. In order to calculate
the excitation branches $\Omega(k)$ in this approximation, we search for the
solutions implicitly given by [45, 63]

$\displaystyle[G_{0}^{R}]^{-1}(k,\Omega(k))\begin{pmatrix}\Psi_{c}\\\
\Psi^{*}_{c}\end{pmatrix}=0\;.$ (13)

By Fourier transforming Eq. (11), we find

$\displaystyle\Omega_{\pm}(k)=-\text{i}\kappa_{0}\mp 2J\cos(k)\;.$ (14)

This approach allows us to identify $\kappa_{0}$ with the decay rate at weak
interactions corresponding to the broadening of the spectral line, as well as
the excitation dispersion $\Omega_{+}$, corresponding to the free-particle
dispersion in a lattice. Inverting the matrix (11) in Fourier space and
selecting the relevant element (10), we obtain the retarded Green’s function

$G_{0}^{R}(k,\omega)=\frac{1}{\omega+2J\cos k+\text{i}\kappa_{0}}\;,$ (15)

which has poles for $\omega=\Omega_{+}(k)$. Notice that the second branch
$\Omega_{-}(k)$ has zero spectral weight in this approximation. Hence, in the
open system under investigation, there is no ghost branch predicted at
quadratic order in the fluctuations, due to the absence of a quasi Bose-
Einstein condensate. To include the effect of interactions, we turn to
numerical simulations in Secs. 4 and 5 below. Prior to this, we present the
steady-state density matrix.

##  3 Exact properties of the non-equilibrium steady state

###  3.1 Steady state density matrix

It has been shown in Ref. [58] that the ansatz

$\displaystyle\rho_{\text{NESS}}=\frac{1}{\mathcal{N}}\sum_{N=0}^{\infty}z^{N}\mathbb{1}_{N}\;,\qquad
z=\frac{\gamma_{p}}{\gamma_{l}}\;,$ (16)

is a NESS of the Lindblad equation (2) for $z<1$. Notice that this state
corresponds to an infinite-temperature grand-canonical partition function with
fugacity $z$. The normalisation $\mathcal{N}$ in Eq. (16) is given by

$\displaystyle\mathcal{N}=\sum_{N=0}^{\infty}D^{B}_{L,N}z^{N}\;,\qquad
D^{B}_{L,N}=\operatorname{Tr}\mathbb{1}_{N}=\begin{pmatrix}\begin{pmatrix}L\\\
N\end{pmatrix}\end{pmatrix}\equiv\begin{pmatrix}L-1+N\\\ N\end{pmatrix}\;,$
(17)

where we used the fact that the dimension of the Hilbert space for $N$ bosons
in a system of length $L$ is given by the number of $L$-tuples of total length
$N$. This is mathematically equivalent to the ‘star and bar’ problem of
identifying in how many ways one can place $L-1$ bars (separating sites)
between $N$ stars (occupied sites). The normalisation can be further
simplified to $\mathcal{N}=(1-z)^{-L}$. A generalised formula for the
normalisation in a truncated Hilbert space can be found in Appendix B.

We remark that, in this idealised model, there is no high-energy cutoff on the
pump intensity. We refer to Ref. [58] for a model with frequency dependent
effective jump operators. The infinite dimensional Hilbert space of the model
as written above renders the Hamiltonian unbounded. We point out that in the
numerical calculations, we have to introduce a cutoff for the local
occupation, which gives an upper bound for the Hamiltonian and ensures that
the time evolution is a CPT map as well as the existence of a NESS [43]. We
have generalised the solution (16) for $\rho_{\text{NESS}}$ to the case of
driven-dissipative free fermions and hard-core bosons on the lattice. The
results are given in Appendix C.

###  3.2 Equal-time correlations in the steady state

The simple form of the steady state solution (16) allows one to exactly
calculate the equal-time correlation functions of the bosonic field operators
in the NESS. We here present the solution for the one-body density matrix
$\langle b_{i}^{\dagger}b_{j}\rangle$ and for the variance of the local
particle occupation number $\text{VAR}(n_{i})$, as they give insights in the
occupation of the system and fluctuations around it. The higher moments can be
easily obtained following the same method. The equal-time correlation function
in the NESS is given by

$\displaystyle\langle b_{i}^{\dagger}b_{j}\rangle_{\text{NESS}}$
$\displaystyle=\operatorname{Tr}\\{b_{i}^{\dagger}b_{j}\rho_{\text{NESS}}\\}$
$\displaystyle=(1-z)^{L}\sum_{N=0}^{\infty}z^{N}\sum_{\\{m_{r}\\}_{r=1,...,L}}^{{}^{\prime}}\langle{\\{m_{r}\\}}|b_{i}^{\dagger}b_{j}|{\\{m_{r}\\}}\rangle\;,$
(18)

where we inserted a complete basis in the Fock space $|\\{m_{r}\\}\rangle$ for
$N$ particles on $L$ sites, with
$b_{r}^{\dagger}b_{r}|\\{m_{r}\\}\rangle=m_{r}|\\{m_{r}\\}\rangle$, and where
$\sum^{\prime}_{\\{m_{r}\\}_{r=1,...,L}}$ runs over the sets $\\{m_{r}\\}$
such that $\sum_{r}m_{r}=N$. We note that, due to the diagonal form of
$\rho_{\text{NESS}}$, the matrix elements
$\langle{\\{m_{r}\\}}|b_{i}^{\dagger}b_{j}|{\\{m_{r}\\}}\rangle$ are non-
vanishing only if $i=j$. After some combinatorial manipulations (see analogous
derivation in Appendix A.3), we find that

$\displaystyle\langle b_{i}^{\dagger}b_{j}\rangle_{\text{NESS}}$
$\displaystyle=\delta_{ij}\frac{z}{1-z}\;,$ (19)

where one must assume that $z<1$, i.e. $\gamma_{l}>\gamma_{p}$. The loss rate
must be larger than the pump rate, in agreement with the Keldysh analysis of
Sec. 2.3. Note that since $\langle b_{i}^{(\dagger)}\rangle_{\text{NESS}}=0$,
the correlation function (19) corresponds also to the connected correlation
function.

With a similar method, we next calculate the fluctuations of the local
particle number density. This quantity is useful in order to assess the effect
of the truncation of the local Hilbert space in the numerical calculations.
The variance is given by

$\displaystyle\text{VAR}(n_{i})\coloneqq\langle n_{i}^{2}\rangle-\langle
n_{i}\rangle^{2}=\frac{z}{(1-z)^{2}}\;.$ (20)

We note that at vanishing pump rate $z\rightarrow 0$ the local number
fluctuations vanish, while they diverge for $z\rightarrow 1$. In the latter
regime the system is close to the instability point given by the change of
curvature of the Keldysh potential. The ratio of the number fluctuation
$\Delta(n_{i})=\sqrt{\text{VAR}(n_{i})}$ to the average occupancy $\langle
n_{i}\rangle$ is given by $\frac{1}{\sqrt{z}}$ and decreases for large
occupation. This calculation provides a good estimate on the error made by
truncating the Hilbert space to $N_{s}$ states per lattice site. For example,
choosing $z=1/10$ ($z=1/5$) gives an average occupation of $1/9$ ($1/4$) with
the number fluctuation of $\sqrt{10}/9$ ($\sqrt{5}/4$). Hence, for these
parameters, a truncation of the Hilbert space to $N_{s}=3$ allowing only for
the states $\\{|0\rangle,|1\rangle,|2\rangle\\}$ on each site is justified.
This will be confirmed by the numerics in the following.

For the total particle number
$N=\langle\sum_{j=1}^{L}b_{j}^{\dagger}b_{j}\rangle$ in the system, one can
solve the full time evolution analytically. The result reads

$\displaystyle\langle N(t)\rangle=(N_{0}-N_{\rm
NESS})\,\text{e}^{-(\gamma_{l}-\gamma_{p})t}\,+N_{\rm NESS}\;,$ (21)

where $N_{0}=N(t=0)$ and $N_{\rm NESS}=\frac{z}{1-z}L$ is the total steady
state occupation. The details can be found in Appendix A.2. Notice that there
exists a NESS if and only if the pump rate is chosen strictly smaller than the
loss rate, i.e. $z<1$, consistently with the previous considerations.

##  4 Correlation function in real time

We report in this section our results for the non-equal time two-point
correlation function.

###  4.1 Analysis of the parameter space

The original model (2) depends on four couplings:
$\\{J,U,\gamma_{l},\gamma_{p}\\}$. By rescaling time by the typical loss time
$\tau_{l}=1/\gamma_{l}$, we are left with the three independent parameters
$\\{J/\gamma_{l},U/\gamma_{l},z\\}$, with $z=\gamma_{p}/\gamma_{l}$ for which
$J\in\mathbb{R}^{+}_{0}$, $U\in\mathbb{R}^{+}_{0}$ and $z\in[0,1)$ for a NESS
to exist. The parameter space and known limiting cases are illustrated in Fig.
2. The plane $z=0$ corresponds to a purely lossy system with empty steady
state. The plane $J=0$ corresponds to the case of absence of tunneling among
the sites – in this case the problem maps to a set of single-site dissipative
Kerr resonator which has an exact solution [59]. The plane $U=0$ corresponds
to the non-interacting limit which can be solved exactly using the non-
interacting (Gaussian) Keldysh action. In the following, we fix a value for
$z$ and explore numerically the two-point correlation function in the
interaction/tunnel energy plane.

Figure 2: Sketch of the parameter space under investigation:
$J/\gamma_{l}\in\mathbb{R}^{+}_{0}$, $U/\gamma_{l}\in\mathbb{R}^{+}_{0}$ and
$z\in[0,1)$; the plane $z=0$ (light green) describes a purely lossy system
with empty steady state; the plane $J=0$ (yellow) corresponds to a local
problem: the solution to $L$ independent dissipative Kerr resonators is known
exactly; the plane $U=0$ is the non-interacting limit which can also be solved
exactly. The inserted plane for finite $\tilde{z}$ indicates the numerically
investigated phase space.

Figure 3: Left panel: space-time evolution of the real part of the retarded
Green’s function (dimensionless). Right panel: cuts for various spatial
differences as indicated on the legend as a function of time (in units of
$\tau_{l}=1/\gamma_{l}$) and extracted exponential decay rate (solid black
line). All results are obtained by exact diagonalisation with parameters
$U=50$, $J=10$, $z=0.2$, $L=8,\,N_{s}=3$.

###  4.2 Time correlations

Using the analytical expression for the NESS density matrix, we have performed
numerical calculations using exact diagonalisation (see details in Appendix
F.1) to obtain the retarded Green’s function in space-time. Fig. 3 displays an
example of such a Green’s function, which displays the emergence of two main
features. First, the Green’s function decays exponentially in time, and
second, it oscillates at well-defined frequencies.

Figure 4: Left panel: decay rate $\kappa$ of the retarded Green’s function for
$z=0.1$ in the $\\{U,J\\}$ plane with $U$, $J$ in units of $\gamma_{l}$. Right
panel: cut of the decay rate as a function of $U/\gamma_{l}$ at fixed
$J/\gamma_{l}=10$. The analytical solution for a single site Eq. (58) is also
shown, as well as the two limits $\kappa_{0}=\kappa(U=0)$ from Eq. (12) and
$\kappa_{\infty}=\kappa(U\rightarrow\infty)$ from Eq. (23). Calculations are
done by exact diagonalisation with $L=4$, $N_{s}=4$.

Quite remarkably, both the decay rate and oscillation frequencies are affected
by the interactions. We next focus on the properties of the decay rate
$\kappa$. The frequencies are discussed in Sec. 5 and Appendix A.4.

###  4.3 Renormalisation of the temporal decay

In order to capture the main features of the non-static quantities in the
NESS, we determine the inverse life time, i.e. the decay $\kappa$ of the
equal-space retarded Green’s function. We investigate how the decay
$\kappa(U,J)$ changes for fixed ratio of pump-to-loss rate $z$. As a benchmark
of our numerical procedure, we recover the known exact results, both in the
limit of weak interaction and in the limit of weak hopping.

To extract the decay rate from the numerical data, we assume that the retarded
Green’s function $G^{R}_{00}(t)$ endows at coincident spatial points a single-
pole form corresponding to

$\displaystyle G^{R}_{\rm sp}(t)=-\text{i}\theta(t)\text{e}^{-\text{i}\Omega
t-\kappa t}\,\;.$ (22)

This leads us to define the decay rate as
$\kappa=-\Re\Big{[}\ln\text{i}G^{R}_{\rm{sp}}(t)\Big{]}/t$ for $t>0$. We thus
represent the numerical data as ${\cal
G}(t)=-\Re\Big{[}\ln\text{i}G^{R}_{00}(t)\Big{]}$ and determine the slope of
${\cal G}$ to obtain an estimate of the decay rate (see Appendix F.2 for
details). The resulting decay rate is depicted in Fig. 4 in the $\\{U,J\\}$
plane. At weak interactions, i.e. for $U/\gamma_{l}\lesssim 1$, we recover the
prediction from the free, quadratic theory of Sec. 2.3, given by Eq. (12). At
strong interactions, when tunneling is negligible with respect to both
interaction energy $U/J\gg 1$ and pump/loss rates
$J\ll\gamma_{l},\,\gamma_{p}$, the decay rate tends to a plateau independent
of the interaction strength which is to the best of our knowledge unknown for
extended systems. For the single-site case, the value of this plateau can be
extracted from the exact solution of Ref. [59] for a single dissipative Kerr
resonator as

$\displaystyle\kappa_{\infty}=\frac{\gamma_{l}+3\gamma_{p}}{2}=\kappa_{0}+2\gamma_{p}\;$
(23)

(see Appendix E for details). In between these two limits, we observe a
renormalisation of the decay smoothly connecting weak and strong interactions,
and weak and large decay rate, similarly to the predictions for the single-
site case [59]. We have checked that our numerical simulation for $L=1$ agrees
well with the exact solution all through the crossover from weak to strong
interactions.

We notice that the decay rate increases both when the interaction strength
increases at fixed $J$, as well as when the tunnel amplitude increases at
fixed (large) $U$. We understand this latter effect as being due to increased
decay possibilities upon allowing for tunneling among particles. For the
parameter regime we could access numerically, we find that the value of the
decay rate for large tunneling amplitude and large interaction strengths
depends on the system size (see Appendix F.4).

##  5 The spectrum

We now focus on the properties of the single-particle excitations on top of
the NESS by analysing the spectral function $A(k,\omega)$ defined in Eq. (7).
For interacting 1D bosons in closed systems, this function displays several
noticeable features as broad spectra with power-law singularities in
correspondence of the Lieb-I and Lieb-II excitation branches [64]. At low
energy, the dispersion of the Goldstone branch is linear, confirming the $z=1$
dynamical critical exponent of the equilibrium model. On the lattice, the
spectral function also displays a third excitation branch associated with the
change of curvature of the single-particle dispersion [17], as well as a
doublon branch at energy close to $U$ [65]. These results are used as a
reference to analyse and discuss the spectral function of the open system.

###  5.1 Spectrum of the open system

Our exact diagonalisation results for the spectral function are shown in Fig.
5. At weak interaction, we find the free spectral line with the Lorenzian
shape in frequency predicted from the field theory calculation in Sec. 2.3.
Two main features emerge at strong interactions. First, differently from the
closed-system case, the low-energy branch is quadratic rather than linear,
with an emerging dispersion branch corresponding to the free particle one. The
background of excitations behind this main branch is due to the $N>1$ particle
sectors in $\rho_{\text{NESS}}$. They provide additional contributions to the
spectral function whose precise shape is difficult to resolve given the system
size. Second, we identify an additional excitation branch centered around
energy $U$, which is the analogue of the doublon branch of the closed-system
case. This branch provides a clear signature of interactions. This is
remarkable as $\rho_{\text{NESS}}$ does not contain information on the
interaction itself and does not depend explicitly on $U$. We provide an
estimate of the position of this excitation branch using a strong-coupling
approach in Sec. 5.2.

The spectral function here obtained has limited resolution due to the size of
the system accessible in numerical diagonalisation, but displays nevertheless
emergent features characteristic also of larger system sizes. This is
illustrated in Sec. 5.3, where we present the calculation of the spectral
function for the closed-system case for two choices of system size.

Figure 5: Spectral function $A(\omega,k)$ for the NESS of a driven-dissipative
Bose-Hubbard model in the weakly (left panel,
$U/\gamma_{l}=1,\,J/\gamma_{l}=10$) and strongly (right panel
$U/\gamma_{l}=500,\,J/\gamma_{l}=50$) interacting regimes. The other
parameters are: $z=0.2,\,L=8,\,N_{s}=3$. The horizontal white lines are set at
$\omega_{b}\coloneqq kU\;,k\in\\{0,...N_{s}-1\\}$ (positions of the branches
for closed system and $J=0$, compare with (27)), and dotted lines at
$\omega=\omega_{b}\pm 2J$.

###  5.2 Position of the excitation spectral lines for small hopping

In order to estimate the position of the excitation branches in the
interacting case, we analytically calculate the oscillation frequency of the
doublon branch of the NESS under unitary time evolution. As shown in Sec. 5.4,
the unitary time evolution well accounts for the position of the peaks of the
spectral function.

We use the one-site model, i.e. we set $J=0$. In this case, the Hamiltonian
becomes local in position space

$\displaystyle H=\frac{U}{2}\sum_{j}n_{j}(n_{j}-1)\;.$ (24)

Upon evaluating the two-point correlation function in the NESS,

$\displaystyle\Big{\langle}b_{j}^{\dagger}(t)b_{0}\Big{\rangle}_{\text{NESS}}$
$\displaystyle=(1-z)^{L}\sum_{N}z^{N}\sum_{\\{m\\}}\langle\\{m\\}|\text{e}^{\text{i}Ht}\,b_{j}^{\dagger}\text{e}^{-\text{i}Ht}\,b_{0}|\\{m\\}\rangle\;,$
(25)

we obtain that the only non-vanishing term is the one with $j=0$ due to the
locality of the strong-coupling Hamiltonian. This term reads

$\displaystyle\langle
m_{0}|\text{e}^{\text{i}\frac{U}{2}n_{0}(n_{0}-1)t}\,b_{0}^{\dagger}\text{e}^{\text{i}\frac{U}{2}n_{0}(n_{0}-1)t}\,b_{0}|m_{0}\rangle$
$\displaystyle=m_{0}\text{e}^{\text{i}U(m_{0}-1)t}\,\;.$ (26)

This shows that there is an oscillation with period depending on the
interaction strength $U$. Altogether, we find

$\displaystyle\Big{\langle}b_{j}^{\dagger}(t)b_{0}\Big{\rangle}_{\text{NESS}}$
$\displaystyle=\delta_{0j}z(1-z)\sum_{m=0}^{\infty}z^{m}(m+1)\text{e}^{\text{i}Umt}\,\;,$
(27)

with $m+1$ corresponding to the occupation number of the site $j=0$. Further
details of the calculation can be found in Appendix A.4. We notice from Eq.
(27) that the position of the first excitation peak for vanishing $J$ and
large $U$ corresponds to the contribution with $m=1$, i.e.
$\omega_{d}(m=1)=U$. This is at the origin of the doublon branch in Fig. 5.
Furthermore, we observe that the number of dominant peaks in the excitation
spectrum at $k=0$ is exactly $N_{s}-1$, as predicted by (27).

###  5.3 Spectral function for the ground state of the Bose-Hubbard model

In order to assess the role of drive and dissipation on the spectral function,
we show here the results for the spectral function of a closed Bose-Hubbard
system at small filling for the same values of parameters, namely interactions
and hopping strength, as well as the same system size. In the closed system,
numerical calculations allow one to tackle larger system sizes, and to explore
the interplay between interaction and finite-size effects in the spectral
function. As shown in Fig. 6, features clearly identified for large system
size are also found at smaller size and serve as guideline to interpret our
results for the open case.

Figure 6: Spectral function $A(k,\omega)$ for the ground state of the Bose-
Hubbard model under unitary evolution (negative values in grey). Left panel:
Exact solution without truncation at finite interaction ($L=12$, $N=2$
particles, $U/J=10$). Right panel: larger system size ($L=24$, $N=6$,
$U/J=10$, with $N_{s}=3$ obtained with DMRG).

Let us comment on the results for the closed system. In the lowest-energy
excitation manifold, we clearly identify the Lieb-I and Lieb-II branches, as
well as the third branch predicted for the lattice case[17]. The predictions
for the dispersion branches in the infinite interaction limit is also shown in
Fig. 6. Notice that in this limit all the branches are slightly shifted
upwards with respect to the finite-$U$ spectrum, since they are centered at
$\omega=\mu-2J$ and the value of the chemical potential $\mu$ at infinite
interactions is larger than the one at finite interactions. At frequencies
$\omega\simeq\mu+U+2J$, we also see in the figure a dispersive doublon branch.

We stress that this result is very different from the one obtained by the
Lindblad evolution presented in the previous section. In the latter case, the
spectral function displays a quadratic behaviour at low momenta, while in the
current case the dispersion at low momenta is linear. This is an illustration
of the fact that the nature of the low-energy excitations involved is
different depending on the state of the system considered. For the ground
state, low-energy excitations are particle-hole excitations on top of an
effective Fermi sphere, while for the NESS state, all excitations are
possible, and single particle excitations dominate the intensity of the
spectral function.

###  5.4 Comparison of unitary evolution and Lindblad evolution of the NESS

Finally, to get further insights on the spectral function of the driven-
dissipative system, we compare the results obtained with the Lindblad
evolution (4) to the ones obtained starting from the same
$\rho_{\text{NESS}}$, but with unitary evolution (5), in the regime where the
coherent parameters are much larger then the incoherent drive and dissipation
$J,\,U\gg 1,z$. The latter can be interpreted as a quenched system: The system
is initially open and in the steady state and then, drive and dissipation are
turned off and the temporal evolution considered is the unitary one. As shown
in Fig. 7 for strong interactions, the spectral functions in both cases
closely resemble each other, as one would expect in the limit of vanishing
incoherent parameters. The broadening of the lines is very small comparatively
to the typical features (depth, width) of the spectrum which are of the order
of the coherent parameters.

Figure 7: Spectral function $A(\omega,k)$ of the Markovian open system
dynamics in the NESS. Left: Under unitary dynamics; Right: Under Lindbladian
dynamics. The parameters are: $U=500,\,J=50,\,z=0.2,\,L=8,\,N_{s}=3$. The
horizontal white lines are as in Fig. 5.

##  6 Conclusions and perspectives

The study of the bosonic quantum gas coupled to an environment is a present
problem of particular interest. Open quantum systems play a crucial role in
the current noisy intermediate-scale quantum (NISQ) era [66], with quantum
science and quantum engineering bringing forth a variety of platforms for
quantum simulations of emergent collective many-body phenomena.
Notwithstanding the exceptional technological advances, photonic platforms
remain intrinsically dissipative and we foresee as particularly important for
the near future the characterisation of the interplay of the open-system
dynamics of these setups with the coherent Hamiltonian dynamics, which will be
characterised by increasingly strong Kerr non-linearities.

In this article we have studied the space-time behaviour of two-point
correlation functions and the spectral function for a driven-dissipative
system of interacting bosons on a lattice, subjected to uniform incoherent
pump and one-body losses. Our starting point has been the analytical
expression for the NESS density matrix of the system derived in Ref. [58]. We
have used several techniques, namely numerical exact diagonalisation, Keldysh
formalism and strong coupling approaches, to provide a complete understanding
of the excitation spectrum on top of the non-equilibrium steady state. We have
found that the decay of the retarded Green’s function is exponential in time,
and the decay rate is renormalised by interactions and increases at both
increasing interactions and tunnel amplitude. The spectral function displays
at low frequency a quadratic-like branch, very different from the linear
branch found for the excitations on top of the ground state of the closed
system, for which we have calculated the spectral function at finite
interactions. At large interactions, we have identified, both in the open and
in the closed system, a dispersive doublon branch at energy of order $U$. By
performing calculations for various system sizes in the closed system case, we
have shown that specific features of the spectrum clearly emerging at large
system sizes are also recognisable in smaller systems, of typical size
accessible via numerical diagonalisation in the open case.

In outlook, on a technical level, it would be interesting to push the
calculations for the spectral function of the open system to larger system
sizes using techniques based on matrix-product states such as those detailed
in Refs. [67, 68]. This would provide conclusive evidence about the presence
of excitations at negative energies (the ‘ghost’ branches) which, according to
Bogoliubov theory, are predicted to be more elusive in the driven-dissipative
case than in the closed one [69, 70, 71] and could not be resolved in our
current numerical calculations. On a fundamental level, our work opens very
interesting perspectives for the possibility of observing the effects of
interactions in out-of-equilibrium systems within actual experiments, where
short one-dimensional lattices of small photonic resonators can be engineered.

## Acknowledgements

We warmly thank I. Carusotto, L. Garbe and S. Diehl for fruitful discussions,
and L. Rosso for significant help in the early stages of the project. MZ
thanks for visiting the LPTMS and using their computational infrastructure.

#### Funding information

AM acknowledges funding from the Quantum-SOPHA project ANR-21-CE47-0009; This
work is part of HQI (www.hqi.fr) initiative and is supported by France 2030
under the French National Research Agency grant number ANR-22-PNCQ-0002; LM
acknowledges funding from the ANR project LOQUST ANR-23-CE47-0006-02; LC
acknowledges support from Institut Universitaire de France (IUF).

##  Appendix A Detailed calculations

###  A.1 Bosonic NESS density matrix

We show in this section that the ansatz (16) is a NESS of the Lindblad
evolution (2). The Hamiltonian conserves the particle number and hence
commutes with each identity operator in each particle number subspace.
Furthermore, the non-hermitian Hamiltonian terms in (2) can be written as

$\displaystyle\rho_{\text{NESS}}\sum_{i}b_{i}^{\dagger}b_{i}=\rho_{\text{NESS}}N=N\rho_{\text{NESS}}=\sum_{i}b_{i}^{\dagger}b_{i}\rho_{\text{NESS}}\;,$
(28)
$\displaystyle\rho_{\text{NESS}}\sum_{i}b_{i}b_{i}^{\dagger}=\rho_{\text{NESS}}(N+L)=\sum_{i}b_{i}b_{i}^{\dagger}\rho_{\text{NESS}}\;,$
(29)

such that we only have to compute the remaining terms

$\displaystyle\sum_{i}b_{i}^{\dagger}\rho_{\text{NESS}}b_{i}$
$\displaystyle=\sum_{i}\sum_{N}\omega_{N-1,L}\sum_{\begin{subarray}{c}\\{l_{r}\\}_{r=1,...,L;N}\\\
\\{n_{r}\\}_{r=1,...,L;N}\end{subarray}}|\\{l_{p}\\}\rangle\langle\\{n_{s}\\}|\sum_{\\{m_{r}\\}_{r=1,...,L;N-1}}\langle\\{l_{p}\\}|b_{i}^{\dagger}|\\{m_{r}\\}\rangle\langle\\{m_{r}\\}|b_{i}|\\{n_{s}\\}\rangle$
$\displaystyle=\sum_{N}\omega_{N-1,L}\sum_{\begin{subarray}{c}\\{l_{r}\\}_{r=1,...,L;N}\\\
\\{n_{r}\\}_{r=1,...,L;N}\end{subarray}}|\\{l_{p}\\}\rangle\langle\\{n_{s}\\}|\langle\\{l_{p}\\}|\sum_{i}b_{i}^{\dagger}b_{i}|\\{n_{s}\\}\rangle$
$\displaystyle=\sum_{N}N\omega_{N-1,L}\mathbb{1}_{N}\;,$ (30)

and similarly,

$\displaystyle\sum_{i}b_{i}\rho_{\text{NESS}}b_{i}^{\dagger}$
$\displaystyle=\sum_{N}(N+L)\omega_{N+1,L}\mathbb{1}_{N}\;.$ (31)

The sum of all these contributions vanishes

$\displaystyle\mathcal{L}_{\text{loss}}[\rho_{\text{NESS}}]+\mathcal{L}_{\text{gain}}[\rho_{\text{NESS}}]$
$\displaystyle=\sum_{N}\mathbb{1}_{N}\Big{[}\gamma_{l}\Big{\\{}(N+L)\omega_{N+1,L}-N\omega_{N,L}\Big{\\}}$
$\displaystyle\qquad\qquad+\gamma_{p}\Big{\\{}N\omega_{N-1,L}-(N+L)\omega_{N,L}\Big{\\}}\Big{]}$
$\displaystyle=\sum_{N}\mathbb{1}_{N}\omega_{N-1,L}\Big{[}\gamma_{l}\Big{\\{}z^{2}(N+L)-zN\Big{\\}}+\gamma_{p}\Big{\\{}N-z(N+L)\Big{\\}}\Big{]}$
$\displaystyle=0\;,$ (32)

which shows that the ansatz (16) is indeed a steady state.

###  A.2 Approach to the NESS

The time evolution of the total particle number
$N=\sum_{j=0}^{L}b_{j}^{\dagger}b_{j}$ is given by the Lindblad master
equation (2). Since the total particle number commutes with the Hamiltonian
$[H,N]=0$, its evolution is purely generated by the dissipative and pump
parts, which yield

$\displaystyle\dot{\langle N\rangle}=(-\gamma_{l}+\gamma_{p})\langle
N\rangle+\gamma_{p}L\;.$ (33)

Starting with a total occupation of $N_{0}\coloneqq N(t=0)$ and depending on
the value of $\gamma_{l}$ w.r.t $\gamma_{p}$, we find the following solutions:

  * •

$\gamma_{p}=\gamma_{l}$: The particle number grows linearly, there is no
steady state. It can be interpreted as a bosonic enhancement, following

$\displaystyle\langle N\rangle(t)=\gamma_{p}Lt+N_{0}\;.$ (34)

  * •

$\gamma_{p}>\gamma_{l}$: The particle number grows exponentially, there is no
steady state.

  * •

$\gamma_{p}<\gamma_{l}$: The particle number changes according to

$\displaystyle\langle N\rangle(t)=(N_{0}-N_{\rm
NESS})\text{e}^{(-\gamma_{l}+\gamma_{p})t}\,+N_{\rm NESS}\;,$ (35)

where the steady state is given by

$\displaystyle N_{\rm NESS}=\frac{z}{1-z}L\;.$ (36)

###  A.3 Equal-time correlation functions

We start from equation (18). If the number of particles in position $j$ is
$k$, then all other $N-k$ particles are distributed over the remaining $L-1$
sites. Hence, we have

$\displaystyle\sum_{\\{m_{r}\\}_{r=1,...,N}}\langle{\\{m_{r}\\}}|b_{i}^{\dagger}b_{j}|{\\{m_{r}\\}}\rangle$
$\displaystyle=\delta_{ij}\sum_{k=1}^{N}k\begin{pmatrix}\begin{pmatrix}L-1\\\
N-k\end{pmatrix}\end{pmatrix}$
$\displaystyle=\delta_{ij}\sum_{k=0}^{N-1}(k+1)\begin{pmatrix}\begin{pmatrix}L-1\\\
N-1-k\end{pmatrix}\end{pmatrix}\;.$ (37)

Summing over $N$, we get

$\displaystyle\sum_{N=1}^{\infty}z^{N}\sum_{k=0}^{N-1}(k+1)\begin{pmatrix}\begin{pmatrix}L-1\\\
N-1-k\end{pmatrix}\end{pmatrix}$
$\displaystyle=z\sum_{N=0}^{\infty}z^{N}\sum_{k=0}^{N}(k+1)\begin{pmatrix}\begin{pmatrix}L-1\\\
N-k\end{pmatrix}\end{pmatrix}$
$\displaystyle=z\sum_{l=0}^{\infty}(l+1)z^{l}\sum_{M=0}^{\infty}z^{M}\begin{pmatrix}\begin{pmatrix}L-1\\\
M\end{pmatrix}\end{pmatrix}$ $\displaystyle=z\frac{1}{(1-z)^{2}}(1-z)^{-L+1}$
$\displaystyle=\frac{z}{(1-z)^{L+1}}\;,$ (38)

where from the first to the second line, we used

$\displaystyle\sum_{p=0}^{\infty}\sum_{q=0}^{p}f(q,p-q)=\sum_{m=0}^{\infty}\sum_{n=0}^{\infty}f(m,n)\;.$
(39)

The geometrical series in the second line only converges if $|z|<1$. All
together, we find the result (19) as expected from the calculation of the
total particle number in the NESS.

For the variance, we need to calculate the average over quadratic terms in the
local particle number density. Performing similar steps, we obtain

$\displaystyle\langle n_{i}^{2}\rangle_{\text{NESS}}$
$\displaystyle=\sum_{N}\frac{z^{N}}{(1-z)^{-L}}\sum_{\\{m_{r}\\}_{r=1,...,N}}\langle{\\{m_{r}\\}}|b_{i}^{\dagger}b_{i}b_{i}^{\dagger}b_{i}|{\\{m_{r}\\}}\rangle$
$\displaystyle=\frac{z}{(1-z)^{-L}}\sum_{l=0}^{\infty}(l+1)^{2}z^{l}\sum_{M=0}^{\infty}z^{M}\begin{pmatrix}\begin{pmatrix}L-1\\\
M\end{pmatrix}\end{pmatrix}$ $\displaystyle=\frac{z(1+z)}{(1-z)^{2}}\;.$ (40)

###  A.4 Oscillation frequency

The dependence of the oscillation frequency on the interaction $U$ can be best
understood for the single-particle Green’s function under the unitary dynamics
in the limit of strong interactions $U\gg J$. We hence consider the evolution
under

$\displaystyle H=\frac{U}{2}\sum_{j}n_{j}(n_{j}-1)\;.$ (41)

We evaluate the two-point correlation function in the NESS as

$\displaystyle\Big{\langle}b_{j}^{\dagger}(t)b_{0}\Big{\rangle}_{\text{NESS}}$
$\displaystyle=(1-z)^{L}\sum_{N}z^{N}\sum_{\\{m\\}}\langle\\{m\\}|\text{e}^{\text{i}Ht}\,b_{j}^{\dagger}\text{e}^{-\text{i}Ht}\,b_{0}|\\{m\\}\rangle$
(42)
$\displaystyle=(1-z)^{L}\sum_{N}z^{N}\sum_{\\{m\\}}\langle\\{m\\}|\Big{(}\prod_{k}\text{e}^{\text{i}\frac{U}{2}n_{k}(n_{k}-1)t}\,\Big{)}b_{j}^{\dagger}\Big{(}\prod_{l}\text{e}^{-\text{i}\frac{U}{2}n_{l}(n_{l}-1)t}\,\Big{)}b_{0}|\\{m\\}\rangle\;.$

The only non-vanishing contribution is the $j=0$ one, which is just a
consequence of the density matrix being diagonal in the Fock basis. (We can
show the same oscillating behaviour for $j\neq 0$ for a more generic density
matrix). The matrix element factorises and we find

$\displaystyle\langle
m_{0}|\text{e}^{\text{i}\frac{U}{2}n_{0}(n_{0}-1)t}\,b_{0}^{\dagger}\text{e}^{\text{i}\frac{U}{2}n_{0}(n_{0}-1)t}\,b_{0}|m_{0}\rangle$
$\displaystyle=m_{0}\text{e}^{\text{i}\frac{U}{2}m_{0}(m_{0}-1)t}\,\text{e}^{-\text{i}\frac{U}{2}(m_{0}-1)(m_{0}-2)t}\,$
$\displaystyle=m_{0}\text{e}^{\text{i}U(m_{0}-1)t}\,\,,$ (43)

which indicates an overall oscillation with a period depending on the
interaction strength $U$. We finally obtain

$\displaystyle\Big{\langle}b_{0}^{\dagger}(t)b_{0}\Big{\rangle}_{\text{NESS}}$
$\displaystyle=(1-z)^{L}\sum_{N=0}^{\infty}z^{N}\sum_{m_{0}=0}^{N}m_{0}\text{e}^{\text{i}U(m_{0}-1)t}\,\begin{pmatrix}\begin{pmatrix}L-1\\\
N-m_{0}\end{pmatrix}\end{pmatrix}$
$\displaystyle=z(1-z)^{L}\sum_{N=0}z^{N}\sum_{m_{0}=0}^{N}(m_{0}+1)\text{e}^{\text{i}Um_{0}t}\,\begin{pmatrix}\begin{pmatrix}L-1\\\
N-m_{0}\end{pmatrix}\end{pmatrix}$
$\displaystyle=z(1-z)^{L}\sum_{k=0}^{\infty}\sum_{l=0}^{\infty}z^{l}z^{k}(k+1)\text{e}^{\text{i}Ukt}\,\begin{pmatrix}\begin{pmatrix}L-1\\\
l\end{pmatrix}\end{pmatrix}$
$\displaystyle=z(1-z)\sum_{k=0}^{\infty}z^{k}(k+1)\text{e}^{\text{i}Ukt}\,\;,$
(44)

with $k$ being the occupation minus one of site zero (relabeled), where we
used the infinite double sum rule (39).

##  Appendix B Normalisation with truncation of the local Hilbert space

This section gives the normalisation of the NESS density matrix in each
$N$-particle sector depending on the system size $L$ and on the local Hilbert
space dimension $d_{l}$. For this, one has to solve a restricted star
(balls/particles) and bar (box separations) problem which is well known in
combinatorics

$\displaystyle\sum_{i=1}^{L}x_{i}=N\,,$ (45)

and $x_{i}\in{0,...,N}$. The latter condition is equivalent to the existence
of an upper bound for the integers $x_{i}<d_{l}$ $\forall i$. This problem can
be solved by using the inclusion-exclusion principle on the lower-bound
integer sum. The latter is given by the previous sum with the constraint
$x_{i}\geq a_{i}$. We substitute $x_{i}^{\prime}\coloneqq x_{i}-a_{i}$, such
that the modified problem is the unbounded one with $x_{i}^{\prime}\geq 0$:

$\displaystyle\sum_{i}x_{i}^{\prime}=N-\sum_{i}a_{i}\;.$ (46)

The number of all possible solutions (as for the bosons) is given by

$\displaystyle\begin{pmatrix}L+N-1\\\ L-1\end{pmatrix}\;,$ (47)

then we subtract all the cases where there is at least one box with $x_{i}\geq
d_{l}$ (solving the lower-bound integer sum):

$\displaystyle\begin{pmatrix}L\\\ 1\end{pmatrix}\begin{pmatrix}L+N-1-d_{l}\\\
L-1\end{pmatrix}\;,$ (48)

and finally add the two set intersections, such that

$\displaystyle
D^{d_{l}}_{L,N}\coloneqq\\#\text{states}(d_{l},L,N)=\sum_{k=0}^{k(d_{l}-1)\leq
N}(-1)^{k}\begin{pmatrix}L\\\ k\end{pmatrix}\begin{pmatrix}L+N-1-kd_{l}\\\
N-1\end{pmatrix}\;.$ (49)

##  Appendix C Steady state solution for fermions

We consider the Lindblad master equation (2) where now the annihilation and
creation operators are fermionic ones $b_{i}\coloneqq
c_{i},b_{i}^{\dagger}\coloneqq c_{i}^{\dagger}$, obeying the anti-commutation
relations $\\{c_{i},c_{j}^{\dagger}\\}=\delta_{ij}$ and
$\\{c_{i},c_{j}\\}=0=\\{c_{i}^{\dagger},c_{j}^{\dagger}\\}$. The unitary
evolution is given by the Fermi-Hubbard model.

###  C.1 NESS density matrix

In finite dimensional Hilbert spaces, there always exists a steady state
solution of the Lindbladian superoperator [43]. Starting from the same ansatz
(16) as in the main text, we obtain from the normalisation of the state that

$\displaystyle\mathcal{N}^{F}_{L}=\frac{1}{(1+\gamma)^{L}}\;,$ (50)

where we used that the dimension of the Hilbert space for $N$ fermions in a
system of length $L$ corresponds to $D^{F}_{L,N}=\binom{L}{N}$. This is
exactly the size of the unit matrix
$\operatorname{Tr}\mathbb{1}_{N}=D^{F}_{L,N}$. This provides the NESS density
matrix for fermions.

###  C.2 Equal-time two-point correlation functions

We assume that the sums for fermions run only over non-equal elements. The
correlation function is given by

$\displaystyle\langle c_{i}^{\dagger}c_{j}\rangle_{\text{NESS}}$
$\displaystyle=\sum_{N=0}^{L}\frac{z^{N}}{(1+z)^{L}}\sum_{\\{m_{r}\\}_{r=1,...,N}}\langle\\{m_{r}\\}|c_{i}^{\dagger}c_{j}|\\{m_{r}\\}\rangle$
$\displaystyle=\delta_{i,j}\sum_{N=0}^{L}\frac{z^{N}}{(1+z)^{L}}\binom{L-1}{N-1}$
$\displaystyle=\delta_{i,j}\frac{z}{1+z}\;,$ (51)

where we introduced a complete basis in position space $\\{m_{r}\\}$ for $N$
particles and evaluated the sum as

$\displaystyle\sum_{\\{m_{r}\\}_{r=1,...,N}}\langle{\\{m_{r}\\}}|c_{i}^{\dagger}c_{j}|{\\{m_{r}\\}}\rangle$
$\displaystyle=\sum_{i,\\{m_{r}\\}_{r=1,...,N-1}\in\\{1,...,L\\}/\\{i\\}}\langle{\\{m_{r}\\}}|c_{i}^{\dagger}c_{j}|{\\{m_{r}\\}}\rangle\delta_{i,j}$
$\displaystyle=\binom{L-1}{N-1}\delta_{i,j}\;.$ (52)

The variance is given by

$\displaystyle\langle n_{i}^{2}\rangle_{\text{NESS}}$
$\displaystyle=\sum_{N}\frac{z^{N}}{(1+z)^{L}}\sum_{\\{m_{r}\\}_{r=1,...,N}}\langle{\\{m_{r}\\}}|c_{i}^{\dagger}c_{i}c_{i}^{\dagger}c_{i}|{\\{m_{r}\\}}\rangle$
$\displaystyle=\frac{z}{1+z}\;.$ (53)

Since the number of fermions in state $i$ is the same as its squared value,
the result of this operator is the same as the first moment. This leads to

$\displaystyle\text{VAR}(n_{i})=\langle n_{i}^{2}\rangle-\langle
n_{i}\rangle^{2}=\frac{z}{(1+z)^{2}}\;.$ (54)

This function has a maximum variance of $1/4$ for $\gamma_{p}=\gamma_{l}$. For
either limit of pump much bigger/smaller than the loss, the variance tends to
zero.

##  Appendix D Steady-state solution for hard-core bosons

We consider a periodic chain of hard-core bosons, which can be mapped by the
transformation due to Holstein and Primakoff [72] to the XX
spin-$\frac{1}{2}$-chain given by the Hamiltonian

$\displaystyle
H=-2J\sum_{j=1}^{L}(S^{x}_{j}S^{x}_{j+1}+S^{y}_{j}S^{y}_{j+1})\;,$ (55)

with $S^{\mu}=\frac{1}{2}\sigma^{\mu}$ where $\\{\sigma^{\mu}\\}_{\mu=x,y,z}$
are the Pauli matrices. We define the first state as the spin-up state
$|\uparrow\rangle=|+\rangle=|0\rangle$ and note that the operator
$S^{+}_{j}S^{-}_{j}$ acts as a local density operator for the spin-up state at
site $j$, where $S^{+}=S^{x}+\text{i}S^{y}$ and $S^{-}=S^{x}-\text{i}S^{y}$,
which act as spin-flip operators.

The pump and loss of hard-core bosons hence correspond to spin flips in the
XX-model, implying that the resulting Lindbladian is of the form

$\displaystyle\mathcal{L}[\rho]=-\text{i}[H,\rho]+\gamma_{l}\sum_{i=1}^{L}\big{(}{S^{-}}_{i}\rho{{S^{-}}_{i}}^{\dagger}-\frac{1}{2}\\{{{S^{-}}_{i}}^{\dagger}{S^{-}}_{i},\rho\\}\big{)}+\gamma_{p}\sum_{i=1}^{L}\big{(}{S^{+}_{i}}\rho{S^{+}}_{i}^{\dagger}-\frac{1}{2}\\{{S^{+}}_{i}^{\dagger}{{S^{+}}_{i}},\rho\\}\big{)}\,.$
(56)

We make the following ansatz for the non-equilibrium steady state

$\displaystyle\rho_{\text{NESS}}=\frac{1}{\mathcal{N}_{F}}\sum_{N=0}^{L}\Big{(}\frac{\gamma_{p}}{\gamma_{l}}\Big{)}^{N}\mathbb{1}_{N}\;,$
(57)

where the identity operator is the sum over the complete Fock basis for the
$N$-particle sector. The number of states in the system of hard-core bosons is
identical to the one of the fermionic system, which explains the identical
normalisation $\mathcal{N}_{F}$ and structure of the solution. One can show
that the Hamiltonian conserves the particle number, or equivalently the number
of up-spins, and hence it commutes with the steady state density matrix. The
dissipative part of the Lindbladian evolution vanishes with the same reasoning
as for the bosons and fermions.

##  Appendix E Decay rate at large interactions

For $J=0$, each site in the Keldysh action (2.3) decouples. Hence, in the
strong interaction limit $U\rightarrow\infty$, one is left with solving the
problem of $L$ independent dissipative Kerr resonators. This is a well-known
exactly solvable problem [59]. One can exploit a weak symmetry to block-
diagonalise the Lindbladian and then introduce an operator which renders the
Lindbladian quadratic in the superoperators. Equivalently, one can introduce a
time-dependent gauge transformation to obtain a quadratic action as shown in
[50]. The Kerr non-linearity then transforms into a fluctuating frequency
depending on the number density which in turn leads to dephasing. The result
for the retarded Green’s function is given by [59]

$\displaystyle
G^{R}(t)=-\text{i}\theta(t)\frac{\text{e}^{\text{i}Ut+\frac{\gamma_{l}-\gamma_{p}}{2}t}\,}{\Big{(}\cosh{\frac{\Gamma}{2}t}+R_{1}\sinh{\frac{\Gamma}{2}t}\Big{)}^{2}}\;,$
(58)

where

$\displaystyle\Gamma\coloneqq\sqrt{(\gamma_{l}-\gamma_{p})^{2}-U^{2}+2\text{i}U(\gamma_{l}+\gamma_{p})}\;,\qquad
R_{1}\coloneqq\frac{1}{\Gamma}\Big{[}(\gamma_{l}-\gamma_{p})+\text{i}\frac{U(\gamma_{l}+\gamma_{p})}{\gamma_{l}-\gamma_{p}}\Big{]}\;.$
(59)

We can expand the term (assuming $t>0$)

$\displaystyle\kappa_{\infty}t$
$\displaystyle=-\Re\Big{[}\ln\text{i}G^{R}(t)\Big{]}$
$\displaystyle=-\Re\Big{[}\text{i}Ut+\frac{\gamma_{l}-\gamma_{p}}{2}t\Big{]}+2\Re\ln\Big{[}\cosh{\frac{\Gamma}{2}t}+R_{1}\sinh{\frac{\Gamma}{2}t}\Big{]}$
$\displaystyle=-\frac{\gamma_{l}-\gamma_{p}}{2}t+t\Re\Gamma+2\Re\ln\Big{[}1+R_{1}+\text{e}^{-\Gamma
t}\,(1-R_{1})\Big{]}\;.$ (60)

The last term is suppressed in the large $t$ limit. Expanding the second term
for $U\gg J$, we find that

$\displaystyle\Re\,\Gamma=\gamma_{l}+\gamma_{p}+\mathcal{O}{(U^{-2})}\;.$ (61)

such that the renormalised decay for $U\gg J$ reads

$\displaystyle\kappa_{\infty}=\frac{\gamma_{l}+3\gamma_{p}}{2}=\kappa_{0}+2\gamma_{p}\;.$
(62)

##  Appendix F Details on the numerical calculations

###  F.1 Exact diagonalisation

Numerical results were obtained by exact diagonalisation of the Hamiltonian or
Lindbladian in the Fock basis of particle occupation numbers per site and the
time evolution implemented in the qutip library [73, 74]. In the case of the
closed-system evolution, we additionally use the conservation of particle
number, as explained in section F.5.

###  F.2 Fit of the decay

In order to extract the decay of the retarded Green’s function (3) in the
NESS, we assume that it is of the form given by the single-pole ansatz (22)
and define

$\displaystyle{\cal G}(t)\coloneqq-\Re\ln\text{i}G^{R}_{00}(t)\;.$ (63)

In the non-interacting regime, the single-pole ansatz is exact and one has
${\cal G}(t)=\kappa t$. In the interacting case, the decay is renormalised and
we determine it from the linear trend of the oscillating curve ${\cal G}(t)$.

Restricting the dimension of the local Hilbert space, we obtain results for
all interaction strength by exact diagonalisation. Using this method, we
extract the slope of ${\cal G}(t)$ by fitting a family of linear curves over
different intervals. The multiple scale dependence of the oscillation of the
two-point function makes this fitting procedure necessary. As the result
depends on the chosen interval of the fit, we fit a family of 20 curves in the
interval accessible in order to reach sufficient numerical precision. The
average and variance are shown as an example in Fig. 8.

Figure 8: Logarithm of the retarded Green’s function ${\cal
G}(t)=-\Re\Big{[}\ln\text{i}G_{00}^{R}(0t)\Big{]}$ (dimensionless) obtained
from exact diagonalisation (red) together with the single-pole approximation
(black, solid) as a function of time (in units of $\tau_{l}$) for the site
$j=0$. The corresponding fitted decay rate is indicated by the black solid
line. The parameters are $U=5.623$, $J=1.0$, $z=0.1$, $L=4$, $N_{s}=4$.

###  F.3 Effect of the Hilbert space truncation

We showed in the main text that we recover the analytically predicted
renormalisation of the decay $\kappa(U)$ for small hopping $J$. In order to
benchmark the numerical truncation over the local Hilbert space, we here
investigate the dependence of $\kappa(U)$ for large $J=10$. We show in Fig. 9
the analytical (one-sited) prediction for small $J$ and the numerical results
for $N_{s}\in\\{3,4,5\\}$. For large $U$, we expect a deviation from the one-
sited result and we observe that the reached plateau is independent of the
local Hilbert space truncation. In the weak interaction limit, we interpret
the closeness of the results for $N_{s}=4$, $N_{s}=5$ with the exact
prediction as an indication of numerical convergence.

Figure 9: Decay $\kappa(U)$ for $J=10$, $z=0.1$ and $L=4$, and for different
local truncation of the Hilbert space $N_{s}$ in the numerics. We expect to
recover the non-interacting solution for small $U$ and a deviation from the
one-sited solution in the limit of large interaction.

###  F.4 System size dependence

We further investigate the dependence of the reached plateau value of the
decay for large $U$ and $J$ on the system size. From Fig. 9, we conclude that
a truncation of $N_{s}=3$ is sufficient to reach the converged plateau value.
We hence fix $N_{s}=3$ and vary the system size $L$. The dependence of the
renormalised decay on $L$ is shown in Fig. 10, in a limit of large hopping and
interaction ($J=10$ and $U=100$). We observe that the effective decay
increases monotonously with the system size. Further analytical insight could
be obtained from the perturbative calculation of corrections in $J$ in the
strong interaction limit around the exact result (58).

Figure 10: Decay $\kappa(U)$ for different system sizes $L$. The parameters
are $J=10$, $z=0.1$ and $N_{s}=3$.

###  F.5 Exact diagonalisation for the unitary system

We calculate the Green’s function as the overlap of two time-evolved states.
Since we are interested in the unitary time evolution, the dimensionality of
the problem greatly reduces, and is further diminished by the symmetries
present in the problem (particle number conservation, momentum conservation
and parity conservation). We calculate the Green’s function as

$\displaystyle G^{R}_{j0}(t,0)$
$\displaystyle=-\text{i}\Theta(t)\text{Tr}\Big{\\{}\Big{[}b_{j}(t),b_{0}^{\dagger}(0)\Big{]}\rho_{\text{NESS}}\Big{\\}}$
(64)
$\displaystyle=-\text{i}\Theta(t)\Big{(}G_{j}^{(1)}(t)-G_{j}^{(2)}(t)\Big{)}$
(65)

with

$\displaystyle G_{j}^{(1)}(t)$
$\displaystyle\coloneqq\text{Tr}\Bigg{\\{}\text{e}^{\text{i}Ht}\,b_{j}\text{e}^{-\text{i}Ht}\,b_{0}^{\dagger}\rho_{\text{NESS}}\Bigg{\\}}$
$\displaystyle=\mathcal{N}^{-1}\sum_{N=0}z^{N}\sum_{\\{m_{i}\\}_{N}}\text{Tr}\Bigg{\\{}b_{j}\Big{|}\text{e}^{-\text{i}H_{N+1}t}\,b_{0}^{\dagger}\\{m_{i}\\}_{N}\Big{\rangle}\otimes\Big{|}\text{e}^{-\text{i}H_{N}t}\,\\{m_{i}\\}_{N}\Big{\rangle}^{\dagger}\Bigg{\\}}\;,$
(66) $\displaystyle G_{j}^{(2)}(t)$
$\displaystyle\coloneqq\text{Tr}\Bigg{\\{}b_{0}^{\dagger}\text{e}^{\text{i}Ht}\,b_{j}\text{e}^{-\text{i}Ht}\,\rho_{\text{NESS}}\Bigg{\\}}$
$\displaystyle=\mathcal{N}^{-1}\sum_{N=0}z^{N}\sum_{\\{m_{i}\\}_{N}}\text{Tr}\Bigg{\\{}b_{j}\Big{|}\text{e}^{-\text{i}H_{N}t}\,\\{m_{i}\\}_{N}\Big{\rangle}\otimes\Big{|}\text{e}^{-\text{i}H_{N-1}t}\,b_{0}\\{m_{i}\\}_{N}\Big{\rangle}^{\dagger}\Bigg{\\}}\;.$
(67)

In order to reduce the Hilbert space dimension further, we construct the
Hamiltonian in each particle number sector, since $[N,H]=0$. Then, we draw
states from the distribution $\rho_{\text{NESS}}$, which constitute a sequence
of numbers allowed by the local Hilbert space cutoff, and evolve them with the
block Hamiltonian constructed in the reduced basis. This is performed using
qutip [73, 74], which optimises the temporal evolution. We add and subtract
particles (before/ after the evolution) by action on the sequence directly.

###  F.6 Simulations with matrix product states

In order to produce the plots in Fig. 6 we employ a representation of the
many-body state as a matrix-product state [75] and code our algorithm using
the Itensors library [76, 77]. We prepare the ground state considering a
maximum number of bosons per site equal to $2$ on a lattice with periodic
boundary condition and at density $N/L=1/4$. We subsequently compute the
lesser and greater Green’s functions $G^{<}_{0j}(t)=-i\langle
b_{j}^{\dagger}(t)b_{0}\rangle$ and $G^{>}_{0j}(t)=-i\langle
b_{0}b_{j}^{\dagger}(t)\rangle$ and using these data we reconstruct the
spectral functions plotted in the figure. For the time evolution, we fix $J=1$
and use a time step of $\delta t=0.01$.

## References

  * [1] M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch and I. Bloch,  _Quantum phase transition from a superfluid to a Mott insulator in a gas of ultracold atoms_ ,  Nature 415(6867), 39 (2002),  10.1038/415039a.
  * [2] T. Kinoshita, T. Wenger and D. S. Weiss,  _Observation of a one-dimensional Tonks-Girardeau gas_ ,  Science 305(5687), 1125 (2004),  10.1126/science.1100700,  https://www.science.org/doi/pdf/10.1126/science.1100700.
  * [3] B. Paredes, A. Widera, V. Murg, O. Mandel, S. Fölling, I. Cirac, G. V. Shlyapnikov, T. W. Hänsch and I. Bloch,  _Tonks–Girardeau gas of ultracold atoms in an optical lattice_ ,  Nature 429(6989), 277 (2004),  10.1038/nature02530.
  * [4] I. Bloch, J. Dalibard and W. Zwerger,  _Many-body physics with ultracold gases_ ,  Rev. Mod. Phys. 80, 885 (2008),  10.1103/RevModPhys.80.885.
  * [5] E. H. Lieb and W. Liniger,  _Exact analysis of an interacting Bose gas. i. the general solution and the ground state_ ,  Phys. Rev. 130, 1605 (1963),  10.1103/PhysRev.130.1605.
  * [6] C. N. Yang,  _Some exact results for the many-body problem in one dimension with repulsive delta-function interaction_ ,  Phys. Rev. Lett. 19, 1312 (1967),  10.1103/PhysRevLett.19.1312.
  * [7] M. Gaudin,  _Un systeme a une dimension de fermions en interaction_ ,  Physics Letters A 24(1), 55 (1967),  https://doi.org/10.1016/0375-9601(67)90193-4.
  * [8] B. Sutherland,  _Further results for the many-body problem in one dimension_ ,  Phys. Rev. Lett. 20, 98 (1968),  10.1103/PhysRevLett.20.98.
  * [9] S. G. Jackson, G. E. Astrakharchik and M. Olshanii,  _Asymmetric Bethe Ansatz_ (2023), 2311.15155.
  * [10] M. Girardeau,  _Relationship between Systems of Impenetrable Bosons and Fermions in One Dimension_ ,  Journal of Mathematical Physics 1(6), 516 (1960),  10.1063/1.1703687,  %****␣Markovian_Bose_Hubbard_SUBMIT_arxiv.bbl␣Line␣75␣****https://pubs.aip.org/aip/jmp/article-pdf/1/6/516/19055341/516_1_online.pdf.
  * [11] M. D. Girardeau and A. Minguzzi,  _Soluble models of strongly interacting ultracold gas mixtures in tight waveguides_ ,  Phys. Rev. Lett. 99, 230402 (2007),  10.1103/PhysRevLett.99.230402.
  * [12] F. Deuretzbacher, K. Fredenhagen, D. Becker, K. Bongs, K. Sengstock and D. Pfannkuche,  _Exact solution of strongly interacting quasi-one-dimensional spinor Bose gases_ ,  Phys. Rev. Lett. 100, 160405 (2008),  10.1103/PhysRevLett.100.160405.
  * [13] A. G. Volosniev, D. V. Fedorov, A. S. Jensen, M. Valiente and N. T. Zinner,  _Strongly interacting confined quantum systems in one dimension_ ,  Nature Communications 5(1), 5300 (2014),  10.1038/ncomms6300.
  * [14] G. Aupetit-Diallo, G. Pecci, C. Pignol, F. Hébert, A. Minguzzi, M. Albert and P. Vignolo,  _Exact solution for SU(2)-symmetry-breaking bosonic mixtures at strong interactions_ ,  Phys. Rev. A 106, 033312 (2022),  10.1103/PhysRevA.106.033312.
  * [15] E. H. Lieb and F. Y. Wu,  _Absence of Mott transition in an exact solution of the short-range, one-band model in one dimension_ ,  Phys. Rev. Lett. 20, 1445 (1968),  10.1103/PhysRevLett.20.1445.
  * [16] M. Rigol and A. Muramatsu,  _Emergence of quasicondensates of hard-core bosons at finite momentum_ ,  Phys. Rev. Lett. 93, 230404 (2004),  10.1103/PhysRevLett.93.230404.
  * [17] J. Settino, N. Lo Gullo, F. Plastina and A. Minguzzi,  _Exact spectral function of a Tonks-Girardeau gas in a lattice_ ,  Phys. Rev. Lett. 126, 065301 (2021),  10.1103/PhysRevLett.126.065301.
  * [18] O. I. Pâţu,  _Exact spectral function of the Tonks-Girardeau gas at finite temperature_ ,  Phys. Rev. A 106, 053306 (2022),  10.1103/PhysRevA.106.053306.
  * [19] C. D. E. Boschi, E. Ercolessi, L. Ferrari, P. Naldesi, F. Ortolani and L. Taddia,  _Bound states and expansion dynamics of interacting bosons on a one-dimensional lattice_ ,  Phys. Rev. A 90, 043606 (2014),  10.1103/PhysRevA.90.043606.
  * [20] J. Polo, P. Naldesi, A. Minguzzi and L. Amico,  _Exact results for persistent currents of two bosons in a ring lattice_ ,  Phys. Rev. A 101, 043418 (2020),  10.1103/PhysRevA.101.043418.
  * [21] M. Rigol, V. Dunjko and M. Olshanii,  _Thermalization and its mechanism for generic isolated quantum systems_ ,  Nature 452(7189), 854 (2008),  10.1038/nature06838.
  * [22] O. A. Castro-Alvaredo, B. Doyon and T. Yoshimura,  _Emergent hydrodynamics in integrable quantum systems out of equilibrium_ ,  Phys. Rev. X 6, 041065 (2016),  10.1103/PhysRevX.6.041065.
  * [23] B. Bertini, M. Collura, J. De Nardis and M. Fagotti,  _Transport in out-of-equilibrium $XXZ$ chains: Exact profiles of charges and currents_,  Phys. Rev. Lett. 117, 207201 (2016),  10.1103/PhysRevLett.117.207201.
  * [24] L. Bertini, A. De Sole, D. Gabrielli, G. Jona-Lasinio and C. Landim,  _Macroscopic fluctuation theory_ ,  Rev. Mod. Phys. 87, 593 (2015),  10.1103/RevModPhys.87.593.
  * [25] R. Fazio and H. van der Zant,  _Quantum phase transitions and vortex dynamics in superconducting networks_ ,  Physics Reports 355(4), 235 (2001),  https://doi.org/10.1016/S0370-1573(01)00022-9.
  * [26] B. Saxberg, A. Vrajitoarea, G. Roberts, M. G. Panetta, J. Simon and D. I. Schuster,  _Disorder-assisted assembly of strongly correlated fluids of light_ ,  Nature 612(7940), 435 (2022),  10.1038/s41586-022-05357-x.
  * [27] J. Bloch, I. Carusotto and M. Wouters,  _Non-equilibrium Bose–Einstein condensation in photonic systems_ ,  Nature Reviews Physics 4(7), 470 (2022),  10.1038/s42254-022-00464-0.
  * [28] K. Toyoda, Y. Matsuno, A. Noguchi, S. Haze and S. Urabe,  _Experimental realization of a quantum phase transition of polaritonic excitations_ ,  Phys. Rev. Lett. 111, 160501 (2013),  10.1103/PhysRevLett.111.160501.
  * [29] J. Raftery, D. Sadri, S. Schmidt, H. E. Türeci and A. A. Houck,  _Observation of a dissipation-induced classical to quantum transition_ ,  Phys. Rev. X 4, 031043 (2014),  10.1103/PhysRevX.4.031043.
  * [30] M. Ludwig and F. Marquardt,  _Quantum many-body dynamics in optomechanical arrays_ ,  Phys. Rev. Lett. 111, 073603 (2013),  10.1103/PhysRevLett.111.073603.
  * [31] N. Syassen, D. M. Bauer, M. Lettner, T. Volz, D. Dietze, J. J. García-Ripoll, J. I. Cirac, G. Rempe and S. Dürr,  _Strong dissipation inhibits losses and induces correlations in cold molecular gases_ ,  Science 320(5881), 1329 (2008),  10.1126/science.1155309.
  * [32] J. J. García-Ripoll, S. Dürr, N. Syassen, D. M. Bauer, M. Lettner, G. Rempe and J. I. Cirac,  _Dissipation-induced hard-core Boson gas in an optical lattice_ ,  New J. Phys. 11(1), 013053 (2009),  10.1088/1367-2630/11/1/013053.
  * [33] T. Tomita, S. Nakajima, I. Danshita, Y. Takasu and Y. Takahashi,  _Observation of the Mott insulator to superfluid crossover of a driven-dissipative Bose-Hubbard system_ ,  Sci. Adv. 3(12), e1701513 (2017),  10.1126/sciadv.1701513.
  * [34] T. Tomita, S. Nakajima, Y. Takasu and Y. Takahashi,  _Dissipative Bose-Hubbard system with intrinsic two-body loss_ ,  Phys. Rev. A 99, 031601 (2019),  10.1103/PhysRevA.99.031601.
  * [35] D. Rossini, A. Ghermaoui, M. B. Aguilera, R. Vatré, R. Bouganne, J. Beugnon, F. Gerbier and L. Mazza,  _Strong correlations in lossy one-dimensional quantum gases: From the quantum Zeno effect to the generalized Gibbs ensemble_ ,  Phys. Rev. A 103, L060201 (2021),  10.1103/PhysRevA.103.L060201.
  * [36] I. Bouchoule, B. Doyon and J. Dubail,  _The effect of atom losses on the distribution of rapidities in the one-dimensional Bose gas_ ,  SciPost Phys. 9, 44 (2020),  10.21468/SciPostPhys.9.4.044.
  * [37] I. Bouchoule and J. Dubail,  _Breakdown of Tan’s relation in lossy one-dimensional Bose gases_ ,  Phys. Rev. Lett. 126, 160603 (2021),  10.1103/PhysRevLett.126.160603.
  * [38] L. Rosso, A. Biella and L. Mazza,  _The one-dimensional Bose gas with strong two-body losses: the effect of the harmonic confinement_ ,  SciPost Phys. 12, 44 (2022),  10.21468/SciPostPhys.12.1.044.
  * [39] C.-H. Huang, T. Giamarchi and M. A. Cazalilla,  _Modeling particle loss in open systems using Keldysh path integral and second order cumulant expansion_ ,  Phys. Rev. Res. 5, 043192 (2023),  10.1103/PhysRevResearch.5.043192.
  * [40] F. Riggio, L. Rosso, D. Karevski and J. Dubail,  _Effects of atom losses on a one-dimensional lattice gas of hard-core Bosons_ ,  Phys. Rev. A 109, 023311 (2024),  10.1103/PhysRevA.109.023311.
  * [41] J. Maki, L. Rosso, L. Mazza and A. Biella,  _Loss induced collective mode in one-dimensional Bose gases_ (2024), 2402.05824.
  * [42] G. Lindblad,  _On the generators of quantum dynamical semigroups_ ,  Communications in Mathematical Physics 48(2), 119 (1976),  10.1007/BF01608499.
  * [43] A. Rivas and S. F. Huelga,  _Open Quantum Systems: An Introduction_ ,  SpringerBriefs in Physics. Springer,  ISBN 978-3-642-23353-1 978-3-642-23354-8,  10.1007/978-3-642-23354-8 (2012).
  * [44] S. Diehl, A. Micheli, A. Kantian, B. Kraus, H. P. Büchler and P. Zoller,  _Quantum states and phases in driven open quantum systems with cold atoms_ ,  Nature Physics 4(11), 878 (2008),  10.1038/nphys1073.
  * [45] L. M. Sieberer, M. Buchhold and S. Diehl,  _Keldysh field theory for driven open quantum systems_ ,  Reports on Progress in Physics 79(9), 096001 (2016),  10.1088/0034-4885/79/9/096001.
  * [46] T. Prosen,  _Third quantization: a general method to solve master equations for quadratic open Fermi systems_ ,  New Journal of Physics 10(4), 043026 (2008),  10.1088/1367-2630/10/4/043026.
  * [47] T. Prosen and T. H. Seligman,  _Quantization over Boson operator spaces_ ,  Journal of Physics A: Mathematical and Theoretical 43(39), 392004 (2010),  10.1088/1751-8113/43/39/392004.
  * [48] B. Horstmann, J. I. Cirac and G. Giedke,  _Noise-driven dynamics and phase transitions in fermionic systems_ ,  Phys. Rev. A 87, 012108 (2013),  10.1103/PhysRevA.87.012108.
  * [49] T. Barthel and Y. Zhang,  _Solving quasi-free and quadratic Lindblad master equations for open fermionic and bosonic systems_ ,  Journal of Statistical Mechanics: Theory and Experiment 2022(11), 113101 (2022),  10.1088/1742-5468/ac8e5c.
  * [50] A. McDonald and A. A. Clerk,  _Third quantization of open quantum systems: Dissipative symmetries and connections to phase-space and Keldysh field-theory formulations_ ,  Phys. Rev. Res. 5, 033107 (2023),  10.1103/PhysRevResearch.5.033107.
  * [51] S. Wolff, A. Sheikhan and C. Kollath,  _Numerical evaluation of two-time correlation functions in open quantum systems with matrix product state methods: a comparison_ ,  SciPost Phys. Core 3, 010 (2020),  10.21468/SciPostPhysCore.3.2.010.
  * [52] B. Sciolla, D. Poletti and C. Kollath,  _Two-time correlations probing the dynamics of dissipative many-body quantum systems: Aging and fast relaxation_ ,  Phys. Rev. Lett. 114, 170401 (2015),  10.1103/PhysRevLett.114.170401.
  * [53] D. Walls and G. Milburn,  _Quantum Optics_ ,  Springer Berlin Heidelberg,  ISBN 9783540285731 (2008).
  * [54] O. Penrose and L. Onsager,  _Bose-Einstein condensation and liquid Helium_ ,  Phys. Rev. 104, 576 (1956),  10.1103/PhysRev.104.576.
  * [55] E. H. Lieb and D. W. Robinson,  _The finite group velocity of quantum spin systems_ ,  Communications in Mathematical Physics 28(3), 251 (1972),  10.1007/BF01645779.
  * [56] D. Poulin,  _Lieb-Robinson bound and locality for general Markovian quantum dynamics_ ,  Phys. Rev. Lett. 104, 190401 (2010),  10.1103/PhysRevLett.104.190401.
  * [57] G. Stefanucci and R. van Leeuwen,  _Nonequilibrium Many-Body Theory of Quantum Systems: A Modern Introduction_ ,  Cambridge University Press (2013).
  * [58] J. Lebreuilly, M. Wouters and I. Carusotto,  _Towards strongly correlated photons in arrays of dissipative nonlinear cavities under a frequency-dependent incoherent pumping_ ,  Comptes Rendus. Physique 17(8), 836 (2016),  10.1016/j.crhy.2016.07.001.
  * [59] A. McDonald and A. A. Clerk,  _Exact solutions of interacting dissipative systems via weak symmetries_ ,  Phys. Rev. Lett. 128, 033602 (2022),  10.1103/PhysRevLett.128.033602.
  * [60] L. M. Sieberer, A. Chiocchetta, A. Gambassi, U. C. Täuber and S. Diehl,  _Thermodynamic equilibrium as a symmetry of the Schwinger-Keldysh action_ ,  Phys. Rev. B 92, 134307 (2015),  10.1103/PhysRevB.92.134307.
  * [61] V. Gorini, A. Kossakowski and E. C. G. Sudarshan,  _Completely positive dynamical semigroups of n‐level systems_ ,  Journal of Mathematical Physics 17(5), 821 (1976),  10.1063/1.522979,  https://pubs.aip.org/aip/jmp/article-pdf/17/5/821/19090720/821_1_online.pdf.
  * [62] L. M. Sieberer, S. D. Huber, E. Altman and S. Diehl,  _Non-equilibrium functional renormalization for driven-dissipative Bose-Einstein condensation_ ,  Phys. Rev. B 89(13), 134310 (2014),  10.1103/PhysRevB.89.134310,  1309.7027[cond-mat].
  * [63] M. Wouters and I. Carusotto,  _Excitations in a nonequilibrium Bose-Einstein condensate of exciton polaritons_ ,  Phys. Rev. Lett. 99, 140402 (2007),  10.1103/PhysRevLett.99.140402.
  * [64] A. Imambekov, T. L. Schmidt and L. I. Glazman,  _One-dimensional quantum liquids: Beyond the Luttinger liquid paradigm_ ,  Rev. Mod. Phys. 84, 1253 (2012),  10.1103/RevModPhys.84.1253.
  * [65] T. Guaita, L. Hackl, T. Shi, C. Hubig, E. Demler and J. I. Cirac,  _Gaussian time-dependent variational principle for the Bose-Hubbard model_ ,  Phys. Rev. B 100, 094529 (2019),  10.1103/PhysRevB.100.094529.
  * [66] J. Preskill,  _Quantum Computing in the NISQ era and beyond_ ,  Quantum 2, 79 (2018),  10.22331/q-2018-08-06-79.
  * [67] S. R. White,  _Minimally entangled typical quantum states at finite temperature_ ,  Phys. Rev. Lett. 102, 190601 (2009),  10.1103/PhysRevLett.102.190601.
  * [68] A. J. Daley,  _Quantum trajectories and open many-body quantum systems_ ,  Advances in Physics 63(2), 77 (2014),  10.1080/00018732.2014.933502,  https://doi.org/10.1080/00018732.2014.933502.
  * [69] P. Stepanov, I. Amelio, J.-G. Rousset, J. Bloch, A. Lemaître, A. Amo, A. Minguzzi, I. Carusotto and M. Richard,  _Dispersion relation of the collective excitations in a resonantly driven polariton fluid_ ,  Nature Communications 10(1), 3869 (2019),  10.1038/s41467-019-11886-3.
  * [70] I. Frérot, A. Vashisht, M. Morassi, A. Lemaître, S. Ravets, J. Bloch, A. Minguzzi and M. Richard,  _Bogoliubov excitations driven by thermal lattice phonons in a quantum fluid of light_ ,  Phys. Rev. X 13, 041058 (2023),  10.1103/PhysRevX.13.041058.
  * [71] F. Claude, M. J. Jacquet, I. Carusotto, Q. Glorieux, E. Giacobino and A. Bramati,  _Spectrum of collective excitations of a quantum fluid of polaritons_ ,  Phys. Rev. B 107, 174507 (2023),  10.1103/PhysRevB.107.174507.
  * [72] T. Holstein and H. Primakoff,  _Field dependence of the intrinsic domain magnetization of a ferromagnet_ ,  Phys. Rev. 58, 1098 (1940),  10.1103/PhysRev.58.1098.
  * [73] J. Johansson, P. Nation and F. Nori,  _Qutip: An open-source python framework for the dynamics of open quantum systems_ ,  Computer Physics Communications 183(8), 1760 (2012),  https://doi.org/10.1016/j.cpc.2012.02.021.
  * [74] J. Johansson, P. Nation and F. Nori,  _QuTiP 2: A Python framework for the dynamics of open quantum systems_ ,  Computer Physics Communications 184(4), 1234 (2013),  https://doi.org/10.1016/j.cpc.2012.11.019.
  * [75] J. I. Cirac, D. Pérez-García, N. Schuch and F. Verstraete,  _Matrix product states and projected entangled pair states: Concepts, symmetries, theorems_ ,  Rev. Mod. Phys. 93, 045003 (2021),  10.1103/RevModPhys.93.045003.
  * [76] M. Fishman, S. R. White and E. M. Stoudenmire,  _The ITensor Software Library for Tensor Network Calculations_ ,  SciPost Phys. p. 4 (2022),  10.21468/SciPostPhysCodeb.4.
  * [77] M. Fishman, S. R. White and E. M. Stoudenmire,  _Codebase release 0.3 for ITensor_ ,  SciPost Phys. Codebases pp. 4–r0.3 (2022),  10.21468/SciPostPhysCodeb.4-r0.3.