Two-Time Correlation Functions in Dissipative and Interacting Bose–Hubbard Chains

Abstract

A method is presented for the systematic derivation of a hierarchy of coupled equations for the computation of two-time correlation functions of operators for open many-body quantum systems. We show how these systems of equations can be closed in mean-field and beyond approximations. Results for the specific example of the spectral weight functions are discussed. Our method allows one to access the full temporal evolution, not just the stationary solution, of non-equilibrium open quantum problems described by a Markovian master equation.

1. Introduction

A very controlled way of introducing non-trivial dynamics to a Bose–Einstein condensate is to deplete a narrow region of the condensate by localised loss and watch the subsequent evolution of the many-body quantum system. For the experimental situation, as realised, e.g., in Herwig Otts’s group at Kaiserslautern, atoms are ionised in a controlled way by an electronic beam and the produced ions and electrons are quickly extracted [1,2,3,4,5]. Consequently, there is scarcely any backaction onto the remaining atoms in the Bose condensate, provided that the filling factors (particle numbers per site) along the lattice are large. For such a setup, we can assume a Markovian coupling of the system to the environment. A corresponding Markovian master equation was used for such systems, taking into account localised loss in the lattice and phase noise, arising from interactions with the background gas or other experimental imperfections [6,7].

For small systems, typically two to four lattice sites, yet with reasonably large filling factors, we can unravel the Master equation exactly, using quantum jump Monte Carlo simulations [7,8,9]. For larger system sizes, approximative stochastic methods such as the truncated Wigner have been used to compute the single-particle density matrix (SPDM) [10,11] and normally ordered two-time correlation functions can be calculated in the Glauber–Sudarshan and positive-P representations [12,13]. Another beyond-mean-field method was successfully applied to propagate specifically chosen initial conditions, typically fully coherent Bose condensed states. This so-called Bogoliubov Back-reaction (BBR) method is based on a Bogoliubov–Born–Green–Kirkwood–Yvon (BBGKY) hierarchy expansion, consisting of dynamical equations for two-, four-, etc., point equal-time correlation functions, i.e., the expectation values of the two-, four-, etc., body reduced density matrices [6,14]. The interaction term in the Hubbard–Hamiltonian then induces the coupling between these dynamical equations toward higher orders. Typically, one truncates at the second order, approximating the six-point correlator by products of two- and four-point correlation functions, in the way the moments of Gaussian variables would exactly split. This allows one to arrive at a closed system of coupled equations, which we can subsequently solve [7].

In this paper, we present a method to compute non-equal time correlation functions of operators, much in the spirit of the BBGKY hierarchy truncation for equal-time observables, which allows one to take into account higher orders in the fluctuations in a systematic way. To do so, we adapt the quantum regression theorem, extensively used in quantum optics, see e.g., [15,16,17], to interacting ultracold atoms in an open system’s setting. Here, atom–atom interactions play an important role and result in hierarchies of dynamical equations that have to be truncated. More precisely, we want to compute two-time Green functions that are heavily used in solid-state physics, typically for fermonic transport problems, see e.g., [18,19,20]. In contrast to the latter applications, our open quantum systems are not time-translational invariant, which discards working in Fourier (frequency or energy) space and mapping the equations of motions into purely alegabric equations, as done e.g., in Ref. [21].

The paper is organised as follows: Section 2 introduces our open many-body boson system. Section 3 presents the computation scheme for two-time correlation functions, which is then applied in Section 4 in mean-field approximation. Section 5 discusses the next order beyond mean-field, with similar equations presented in Section 6 for the density–density correlation functions. The last Section 7 concludes the paper.

2. Dissipative Finite Bose–Hubbard Chain

We model ultracold bosons in sufficiently deep optical lattices by a tight-binding approximation, using a single-band Bose–Hubbard model. The geometry is assumed to be quasi one-dimensional, corresponding to a cigar-shaped confinement of the atoms (which is much stronger in the radial direction). Then, the dynamics of coherent interacting ultracold atoms tunnelling through an M-site lattice is described by the Bose–Hubbard Hamiltonian [22]:

$${\widehat{H}}_{\mathrm{BH}}=-J\sum _{j=1}^{M-1}({\widehat{a}}_{j+1}^{†}{\widehat{a}}_j^{\phantom{(}}+{\widehat{a}}_j^{†}{\widehat{a}}_{j+1}^{\phantom{(}})+\frac{U}{2}\sum _{j=1}^M{\widehat{a}}_j^{†}{\widehat{a}}_j^{†}{\widehat{a}}_j^{\phantom{(}}{\widehat{a}}_j^{\phantom{(}}.$$

(1)

Here, ${\widehat{a}}_j^{†}$ and ${\widehat{a}}_j^{\phantom{(}}$ denote respectively the creation and annihilation operators at site j, J is the tunnelling rate and U is the interaction strength. The reduced Planck’s constant ℏ is set to one, which corresponds to measuring all energies in frequency units. The natural unit of time is then $J^{-1}$.

Dissipation processes are accounted for by introducing the master equation in Lindblad form [17] together with a suitable Liouville superoperator:

$${\partial }_t\widehat{\rho }(t)=\mathcal{L}\widehat{\rho }(t)=-i\left[{\widehat{H}}_{\mathrm{BH}},\widehat{\rho }(t)\right]+\tilde{\mathcal{L}}\widehat{\rho }(t),$$

(2)

$$\tilde{\mathcal{L}}\widehat{\rho }=-\sum _{j=1}^M\frac{{\gamma }_j}{2}\left({\widehat{L}}_j^{†}{\widehat{L}}_j^{\phantom{(}}\widehat{\rho }+\widehat{\rho }{\widehat{L}}_j^{†}{\widehat{L}}_j^{\phantom{(}}-2{\widehat{L}}_j^{\phantom{(}}\widehat{\rho }{\widehat{L}}_j^{†}\right).$$

(3)

The Lindblad operators ${\widehat{L}}_j$ are chosen depending on the type of relaxation or decoherence process relevant for the system under study. In the following, we restrict to local single-body dissipation as motivated in the introduction, i.e., ${\widehat{L}}_j={\widehat{a}}_j$, and ${\gamma }_j$ is the dissipation rate at site j. This choice of the Lindblad operator, as shown in [7], leads to equations of motion for the SPDM equivalent to the heuristic non-Hermitian discrete nonlinear Schrödinger equation introduced in [23] and successfully applied to the description of localized single-body dissipation processes in Bose–Hubbard chains in good agreement with experimental realizations [3,5]. Further sets of Lindblad operators modelling other processes in Bose–Hubbard chains can be found in the review [24].

The expectation value of some system operator $\widehat{A}$ is then provided by the trace $\langle \widehat{A}(t)\rangle =Tr\{\widehat{A}\widehat{\rho }(t)\}=Tr\{{\widehat{A}}_H(t)\widehat{\rho }(0)\}$. From this, the equivalent Heisenberg representation is defined as ${\widehat{A}}_H(t)=V^{†}(t,0)\widehat{A}$, with the propagator $V(t,t_0)=exp((t-t_0)\mathcal{L})$, for a time-independent Liouvillian. Henceforth, the Heisenberg representation is not specified by an index but simply indicated by the presence of a time argument. In this representation, the time evolution of the Heisenberg operator $\widehat{A}(t)$ is carried out by the adjoint master equation:

$${\partial }_t\widehat{A}(t)=V^{†}(t,0)({\mathcal{L}}^{†}\widehat{A})=i\left[{\widehat{H}}_{\mathrm{BH}},\widehat{A}\right](t)+\left({\tilde{\mathcal{L}}}^{†}\widehat{A}\right)(t),$$

(4)

$${\tilde{\mathcal{L}}}^{†}\widehat{A}=-\sum _{j=1}^M\frac{{\gamma }_j}{2}\left({\widehat{L}}_j^{†}{\widehat{L}}_j^{\phantom{(}}\widehat{A}+\widehat{A}{\widehat{L}}_j^{†}{\widehat{L}}_j^{\phantom{(}}-2{\widehat{L}}_j^{†}\widehat{A}{\widehat{L}}_j^{\phantom{(}}\right)=-\sum _{j=1}^M\frac{{\gamma }_j}{2}\left({\widehat{L}}_j^{†}[{\widehat{L}}_j^{\phantom{(}},\widehat{A}]+[\widehat{A},{\widehat{L}}_j^{†}]{\widehat{L}}_j^{\phantom{(}}\right).$$

(5)

The main observable in the introduced system is the single-particle density matrix (SPDM) ${\sigma }_{j,k}=\langle {\widehat{a}}_j^{†}{\widehat{a}}_k^{\phantom{(}}\rangle =Tr\{{\widehat{a}}_j^{†}{\widehat{a}}_k^{\phantom{(}}\widehat{\rho }\}$. Its diagonal matrix elements give the local populations of the chain, while its off-diagonal elements provide information about the coherence of the state [7]. The quartic interaction term leads to a coupling of the equations of motion verified by the SPDM and higher-order moments. Consequently, computing the SPDM requires a prior truncation of the hierarchy into a closed set of differential equations. The mean-field (MF) approximation keeps only the first order of the hierarchy by neglecting the covariances ${\Delta }_{jmkn}=\langle {\widehat{a}}_j^{†}{\widehat{a}}_m^{\phantom{(}}{\widehat{a}}_k^{†}{\widehat{a}}_n^{\phantom{(}}\rangle -\langle {\widehat{a}}_j^{†}{\widehat{a}}_m^{\phantom{(}}\rangle \langle {\widehat{a}}_k^{†}{\widehat{a}}_n^{\phantom{(}}\rangle$, whereas the BBR close to mean-field method evolves simultaneously the SPDM and the covariance, truncating higher-order moments [6,7,14].

The aim of this article is to provide a method for obtaining multi-time correlation functions for quantum many-body systems from the knowledge of the equal-time correlation functions. Our method works for systems with relaxation channels according to the master Equation (3) and gives access to the transient dynamics of the functions, not just to the stationary-state solutions.

3. Computation Scheme for Two-Point Correlation Functions

Two-point Green functions (GF) provide useful information about many-body systems such as temporal and spatial correlations or information about the response of the system to an external perturbation [20]. As a first example, if one is interested in particular in evaluating the average probability of a particle propagating from some site k at time t to any other site j at time $t^{\prime }=t+\tau$, the retarded Green function provides relevant information. In this case, this correlation function reads:

$$G_{j,k}^R(t+\tau ;t)={\theta }_{\mathrm{H}}(\tau )(G_{j,k}^{>}(t+\tau ;t)-G_{j,k}^{<}(t+\tau ;t)),$$

(6)

where the lesser and greater bosonic GFs are defined as follows:

$$G_{j,k}^{<}(t^{\prime };t)=-i\langle {\widehat{a}}_k^{†}(t){\widehat{a}}_j^{\phantom{(}}(t^{\prime })\rangle ,G_{j,k}^{>}(t^{\prime };t)=-i\langle {\widehat{a}}_j^{\phantom{(}}(t^{\prime }){\widehat{a}}_k^{†}(t)\rangle ,$$

(7)

where the ladder operators are expressed in the above-recalled Heisenberg representation $V^{†}(t,0){\widehat{a}}_j^{(†)}$. In order to compute these last two GFs, we derive a closed set of equations of motion for a dissipative setting by means of a modified version of the quantum regression theorem.

Quantum Regression Hierarchy

Let $\{{\widehat{A}}_i\}$ be a set of arbitrary system operators, e.g., ${\widehat{n}}_i={\widehat{a}}_i^{†}{\widehat{a}}_i^{\phantom{(}}$, and $\mathcal{L}$ a time-independent Liouville superoperator, e.g., the one defined in Equations (4) and (5), the adjoint master equation reads:

$${\partial }_t{\widehat{A}}_i(t)=V^{†}(t,0)({\mathcal{L}}^{†}{\widehat{A}}_i).$$

(8)

Let us now make the general assumption that the adjoint Liouville operator acts on ${\widehat{A}}_i$ in such a way that its expectation value can be rewritten as:

$$Tr\{({\mathcal{L}}^{†}{\widehat{A}}_i)\widehat{\rho }(t)\}=Tr\left\{\left(\sum _{\ell }{T}_{i\ell }^{(1)}{\widehat{A}}_{\ell }+\sum _{\ell ,{\ell }^{\prime }}{K}_{i{\ell }^{\prime }\ell }^{(1)}{\widehat{B}}_{{\ell }^{\prime }}{\widehat{A}}_{\ell }\right)\widehat{\rho }(t)\right\},$$

(9)

where $\{{\widehat{B}}_{\ell }^{\prime }\}$ are operators composed of an even number of creation and annihilation operators, e.g., a power of the density operator. This is the typical relation one gets when dealing with non-quadratic Hamiltonians and/or nonlinear Lindbald operators, and in particular the previously defined Liouvillian, which result in a coupling between the equations of motion satisfied by n-point and higher-order correlation functions. This relation is assumed to hold for any initial density matrix $\widehat{\rho }$ so that one is able to identify both operators in parentheses:

$${\mathcal{L}}^{†}{\widehat{A}}_i=\sum _{\ell }{T}_{i\ell }^{(1)}{\widehat{A}}_{\ell }+\sum _{\ell ,{\ell }^{\prime }}{K}_{i{\ell }^{\prime }\ell }^{(1)}{\widehat{B}}_{{\ell }^{\prime }}{\widehat{A}}_{\ell }.$$

(10)

The K-term requires an extension of the quantum regression theorem (see [25,26] or the Section 5.2.3 in [16] and Section 3.2.3 in [17]). Combining Equations (8) and (10), the equation of motion of the two-point correlation function reads:

$$\begin{array}{cc}\hfill {\partial }_{\tau }\langle {\widehat{A}}_i(t+\tau ){\widehat{A}}_j(t)\rangle & =Tr\{({\mathcal{L}}^{†}{\widehat{A}}_i)V(t+\tau ,t){\widehat{A}}_{j}V(t,0)\widehat{\rho }(0)\}=\langle ({\mathcal{L}}^{†}{\widehat{A}}_i)V(t+\tau ,t){\widehat{A}}_{j}V(t,0)\rangle \hfill \\ & =\sum _{\ell }{T}_{i\ell }^{(1)}\langle {\widehat{A}}_{\ell }(t+\tau ){\widehat{A}}_j(t)\rangle +\sum _{\ell ,{\ell }^{\prime }}{K}_{i{\ell }^{\prime }\ell }^{(1)}\langle ({\widehat{B}}_{{\ell }^{\prime }}{\widehat{A}}_{\ell })(t+\tau ){\widehat{A}}_j(t)\rangle \hfill \\ & =\sum _{\ell }{T}_{i\ell }^{\prime (1)}(t+\tau )\langle {\widehat{A}}_{\ell }(t+\tau ){\widehat{A}}_j(t)\rangle +\sum _{\ell ,{\ell }^{\prime }}{K}_{i{\ell }^{\prime }\ell }^{(1)}\langle (\Delta {\widehat{B}}_{{\ell }^{\prime }}{\widehat{A}}_{\ell })(t+\tau ){\widehat{A}}_j(t)\rangle ,\hfill \end{array}$$

(11)

where $T_{i\ell }^{\prime (1)}(t)=T_{i\ell }^{(1)}+{\sum }_{{\ell }^{\prime }}{K}_{i{\ell }^{\prime }\ell }^{(1)}\langle {\widehat{B}}_{{\ell }^{\prime }}(t)\rangle$ and the central moment operator is defined as $\Delta \widehat{A}=\widehat{A}-\langle \widehat{A}\rangle$. Naturally, the same can be done for $\langle {\widehat{A}}_j(t){\widehat{A}}_i(t+\tau )\rangle$, with the only difference being that ${\widehat{A}}_j$ is then placed at the left end of the moments.

With this proper rewriting, one gets a hierarchy of coupled dynamical equations in the form of the BBGKY hierarchy:

$${\partial }_{\tau }\langle {\widehat{A}}_i(t+\tau ){\widehat{A}}_j(t)\rangle =\sum _{\ell }{T}_{i\ell }^{\prime (1)}(t+\tau )\langle {\widehat{A}}_{\ell }(t+\tau ){\widehat{A}}_j(t)\rangle +\sum _{\ell ,{\ell }^{\prime }}{K}_{i{\ell }^{\prime }\ell }^{(1)}\langle (\Delta {\widehat{B}}_{{\ell }^{\prime }}{\widehat{A}}_{\ell })(t+\tau ){\widehat{A}}_j(t)\rangle ,$$

(12)

$$\begin{array}{cc}\hfill {\partial }_{\tau }\langle (\Delta {\widehat{B}}_{i^{\prime }}{\widehat{A}}_i)(t+\tau ){\widehat{A}}_j(t)\rangle & =\sum _{\ell }{T}_{i\ell }^{\prime (2)}(t+\tau )\langle (\Delta {\widehat{A}}_{{\ell }^{\prime }}{\widehat{A}}_{\ell })(t+\tau ){\widehat{A}}_j(t)\rangle \hfill \\ & +\sum _{\ell ,{\ell }^{\prime },k}{K}_{i^{\prime }i{\ell }^{\prime }\ell k}^{(2)}\langle (\Delta {\widehat{B}}_k\Delta {\widehat{B}}_{{\ell }^{\prime }}{\widehat{A}}_{\ell })(t+\tau ){\widehat{A}}_j(t)\rangle .\hfill \\ & ⋮\hfill \end{array}$$

(13)

This system of differential equations can then be closed by truncating the hierarchy to some order in the fluctuations provided that the time-local averages $\langle {\widehat{B}}_{{\ell }^{\prime }}(t)\rangle$ can be calculated at any time t. In the simplest case, for which $K^{(1)}$ is a zero matrix, one naturally recovers the standard quantum regression theorem expression. Otherwise, the first approximation is to make a mean-field approximation and neglect the covariances of the operators ${\widehat{B}}_{{\ell }^{\prime }}$ and $({\widehat{A}}_{\ell }{\widehat{A}}_j)$ in Equation (12). To perform this approximation, it is convenient to first rewrite Equation (12) in such a way that the order in the fluctuations is explicit:

$$\begin{array}{cc}\hfill {\partial }_{\tau }\langle {\widehat{A}}_i(t+\tau ){\widehat{A}}_j(t)\rangle & =\sum _{\ell }{T}_{i\ell }^{\prime (1)}(t+\tau )\langle {\widehat{A}}_{\ell }(t+\tau ){\widehat{A}}_j(t)\rangle \hfill \end{array}$$

(14)

$$\begin{array}{c}+\sum _{\ell ,{\ell }^{\prime }}{K}_{i{\ell }^{\prime }\ell }^{(1)}\langle \Delta {\widehat{B}}_{{\ell }^{\prime }}(t+\tau )\Delta ({\widehat{A}}_{\ell }(t+\tau ){\widehat{A}}_j(t))\rangle .\hfill \end{array}$$

(15)

From this expression, one observes that here the mean-field approximation amounts to keeping only the line (14), which results in a quantum regression expression with a time-dependent coefficient matrix $T^{\prime (1)}$. The knowledge of $T^{\prime (1)}$ at any time puts then the system into a closed form. Moreover, if one is able to compute $\langle {\widehat{A}}_{{\ell }^{\prime }}\rangle$ and $\langle {\widehat{A}}_{\ell }{\widehat{A}}_j\rangle$ at any time at second order in the fluctuations, for instance by employing the BBR truncation [7], then one can obtain an approximation of the GFs at this order by including the line (15) and computing the Equation (13) of the hierarchy in the BBR approximation.

4. Mean-Field Approximation

We study now the Bose–Hubbard Hamiltonian of Equation (1) with single-body loss putting ${\widehat{L}}_j={\widehat{a}}_j$ in Equation (2). In the mean-field approximation, the equations of motion of the annihilation operators read in that case

$$\begin{array}{cc}\hfill {\partial }_t{\widehat{a}}_j^{\phantom{(}}(t)& =iJ\left({\widehat{a}}_{j+1}^{\phantom{(}}(t)+{\widehat{a}}_{j-1}^{\phantom{(}}(t)\right)-iU{\widehat{a}}_j^{†}(t){\widehat{a}}_j^{\phantom{(}}(t){\widehat{a}}_j^{\phantom{(}}(t)-\frac{{\gamma }_j}{2}{\widehat{a}}_j^{\phantom{(}}(t)\hfill \\ & =\sum _{\ell }\underset{T_{j\ell }}{\underbrace{\left(iJ({\delta }_{j+1,\ell }+{\delta }_{j-1,\ell })-\frac{{\gamma }_{\ell }}{2}{\delta }_{j,\ell }\right)}}{\widehat{a}}_{\ell }^{\phantom{(}}(t)+\sum _{\ell ,{\ell }^{\prime }}\underset{K_{j{\ell }^{\prime }\ell }}{\underbrace{\left(-iU{\delta }_{j,{\ell }^{\prime }}{\delta }_{{\ell }^{\prime },\ell }\right)}}{\widehat{n}}_{{\ell }^{\prime }}^{\phantom{(}}(t){\widehat{a}}_{\ell }(t)\hfill \\ & ={\sum }_{\ell }{T}_{j\ell }^{\prime }(t){\widehat{a}}_{\ell }^{\phantom{(}}(t)+{\sum }_{\ell ,{\ell }^{\prime }}{K}_{j{\ell }^{\prime }\ell }\Delta {\widehat{n}}_{{\ell }^{\prime }}^{\phantom{(}}(t){\widehat{a}}_{\ell }(t),\hfill \end{array}$$

(16)

where

$$T_{j\ell }^{\prime }(t)=iJ({\delta }_{j+1,\ell }+{\delta }_{j-1,\ell })-iU{\delta }_{j,\ell }{n}_{\ell }(t)-\frac{{\gamma }_{\ell }}{2}{\delta }_{j,\ell }.$$

(17)

The simplicity of the equation of motion satisfied by the annihilation operator allows one to get directly a relation of the form of Equation (10), without the need of considering its expectation value and taking advantage of the cyclic permutation invariance of the trace.

Then, the quantum regression theorem yields:

$${\partial }_t{\widehat{a}}_j^{\phantom{(}}(t)\stackrel{\mathrm{MF}}{\approx }{\sum }_{\ell }{T}_{j\ell }^{\prime }(t){\widehat{a}}_{\ell }^{\phantom{(}}(t)\Rightarrow \left\{\begin{array}{c}{\partial }_{\tau }\langle {\widehat{a}}_k^{†}(t){\widehat{a}}_j^{\phantom{(}}(t+\tau )\rangle \stackrel{\mathrm{MF}}{\approx }{\sum }_{\ell }{T}_{j\ell }^{\prime }(t+\tau )\langle {\widehat{a}}_k^{†}(t){\widehat{a}}_{\ell }^{\phantom{(}}(t+\tau )\rangle ,\hfill \\ {\partial }_{\tau }\langle {\widehat{a}}_j^{\phantom{(}}(t+\tau ){\widehat{a}}_k^{†}(t)\rangle \stackrel{\mathrm{MF}}{\approx }{\sum }_{\ell }{T}_{j\ell }^{\prime }(t+\tau )\langle {\widehat{a}}_{\ell }^{\phantom{(}}(t+\tau ){\widehat{a}}_k^{†}(t)\rangle .\hfill \end{array}\right.$$

(18)

Due to the dependence of $T^{\prime }$ on the local density $n_{\ell }=\langle {\widehat{a}}_{\ell }^{†}{\widehat{a}}_{\ell }^{\phantom{(}}\rangle ={\sigma }_{\ell \ell }$ these equations of motion have to be evolved along with those of the single particle density matrix. This leads to the following closed set of differential equations:

$$\begin{array}{cc}\hfill i{\partial }_t{\sigma }_{j,k}& \stackrel{\mathrm{MF}}{\approx }-J({\sigma }_{j,k+1}+{\sigma }_{j,k-1}-{\sigma }_{j+1,k}-{\sigma }_{j-1,k})+U(n_j-n_k){\sigma }_{j,k}-\frac{i}{2}({\gamma }_j+{\gamma }_k){\sigma }_{j,k},\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill i{\partial }_{\tau }{G}_{j,k}^{\lessgtr }(t+\tau ;t)& \stackrel{\mathrm{MF}}{\approx }-J\left(G_{j+1,k}^{\lessgtr }(t+\tau ;t)+G_{j-1,k}^{\lessgtr }(t+\tau ;t)\right)+\left(U{n}_j(t+\tau )-i\frac{{\gamma }_j}{2}\right)G_{j,k}^{\lessgtr }(t+\tau ;t).\hfill \end{array}$$

(20)

These can be evolved in $\tau \ge 0$ for each t. In this way, the first time argument is always later than the second one, but this is not restrictive as the opposite situation can be obtained from the relation:

$$G_{j,k}^{\lessgtr }(t;t+\tau )=-{(G_{k,j}^{\lessgtr }(t+\tau ;t))}^{*}.$$

(21)

Illustrative Results

As an example of the performance of the method presented here, we compute the correlation function defined as $A_{j,k}(t^{\prime },t)=\langle [{\widehat{a}}_j^{\phantom{(}}(t^{\prime }),{\widehat{a}}_k^{†}(t)]\rangle =i(G_{j,k}^{>}(t^{\prime };t)-G_{j,k}^{<}(t^{\prime };t))=i(G_{j,k}^R(t^{\prime };t)-G_{j,k}^A(t^{\prime };t))$ for two different settings. The Fourier transform of $A_{j,k}(t^{\prime },t)$ presents the spectral weight function [20], used e.g., for computing currents in transport setups. This particular correlation function is chosen because its equal-time values and its lower and upper bounds can readily be checked and because its magnitude is equal to that of the retarded Green’s function for $t^{\prime }>t,$ whereas the magnitude of the advanced Green’s function is given by the $t>t^{\prime }$ half-quadrant. In addition, in the noninteracting non-dissipative case, its profile simply consists of periodic oscillations extending to both sides away from the diagonal $t=t^{\prime }$. This harmonic behaviour of $A_{j,k}(t^{\prime },t)$ is related to the mentioned fact that its Fourier transform characterises the spectrum without perturbation in the closed system’s case, in the example shown below of a three-level system.

Figure 1 shows the magnitude of two matrix elements of the spectral weight function of a condensate initially made of $N_0=1000$ atoms loaded into a three-well Bose–Hubbard chain with an initial population imbalance between a low populated central well 2 and its two more populated neighbours 1 and 3 in the presence of strong dissipation occurring at the central site. The values taken at the diagonal $t=t^{\prime }$ correspond to the constant $\left|A_{j,k}\right|=\left|\langle [{\widehat{a}}_j^{\phantom{(}}(t),{\widehat{a}}_k^{†}(t)]\rangle \right|={\delta }_{j,k}$, as expected. In Figure 1a, which depicts the $\left|A_{2,1}\right|$ matrix element, the $t^{\prime }>t$ half-quadrant presents an elevation at early $t^{\prime }$ that drops as $t^{\prime }$ becomes larger than a few $J^{-1}$. This indicates that, on average, the flow of particles from the site 1 to the leaky site 2 is suppressed at times above this typical value.

This suppression is a signature of the so-called quantum Zeno effect, meaning that particles are blocked on average from flowing into the site with dissipation, see e.g., the results and descriptions in [3,24,27]. Dissipation can thus induce an effectively self-trapped regime, considerably lowering the average inter-well tunnelling, although the values of the population imbalance and the interaction strength do not suffice to reach this regime in the absence of dissipation. In Figure 1b, which depicts the $\left|A_{2,1}\right|$ matrix element, whereas, in the noninteracting non-dissipative case, this function shows periodic oscillations from either side of the diagonal around the value 1/2; in this case, the spectral weight function quickly attains a plateau at this value. The fact that it does not take values below one half after a few $J^{-1}$ is another signature of the quenching of the inter-well tunnelling induced by a quantum Zeno effect.

Figure 2 represents the same correlation functions for the same three-well setting but in the absence of dissipation. In this case, the profile of the spectral weight function is very rugged due to the interactions. The correlation function at $t\ne t^{\prime }$ fluctuates around $1/2$ regardless of the matrix element, which indicates that the system is not in a self-trapped regime, as the average inter-site flow of bosons is not suppressed.

From these two examples, one observes that correlations in Bose–Hubbard chains qualitatively differ depending on the presence or absence of dissipation. This as well as the observation of some particular structures in the profiles of the investigated correlation functions can be used to determine and differentiate the dynamical regimes in Bose–Hubbard chain set-ups, where the interplay of interactions and dissipation plays a major role.

5. Beyond Mean-Field Approximation

One of the assets that motivate seeking correlation functions as solutions to a hierarchy of coupled equations of motion is that the contribution of higher order fluctuations can be integrated in a systematic way—for instance, by taking into account the second equation, Equation (13), of the hierarchy and truncating sextic moments as follows [6,14]:

$$\begin{array}{cc}\hfill \langle {\widehat{a}}_j^{†}{\widehat{a}}_m^{\phantom{(}}{\widehat{a}}_k^{†}{\widehat{a}}_n^{\phantom{(}}{\widehat{a}}_r^{†}{\widehat{a}}_s^{\phantom{(}}\rangle \stackrel{\mathrm{BBR}}{\approx }& \langle {\widehat{a}}_j^{†}{\widehat{a}}_m^{\phantom{(}}{\widehat{a}}_k^{†}{\widehat{a}}_n^{\phantom{(}}\rangle \langle {\widehat{a}}_r^{†}{\widehat{a}}_s^{\phantom{(}}\rangle +\langle {\widehat{a}}_j^{†}{\widehat{a}}_m^{\phantom{(}}{\widehat{a}}_r^{†}{\widehat{a}}_s^{\phantom{(}}\rangle \langle {\widehat{a}}_k^{†}{\widehat{a}}_n^{\phantom{(}}\rangle +\langle {\widehat{a}}_k^{†}{\widehat{a}}_n^{\phantom{(}}{\widehat{a}}_r^{†}{\widehat{a}}_s^{\phantom{(}}\rangle \langle {\widehat{a}}_j^{†}{\widehat{a}}_m^{\phantom{(}}\rangle -2\langle {\widehat{a}}_j^{†}{\widehat{a}}_m^{\phantom{(}}\rangle \langle {\widehat{a}}_k^{†}{\widehat{a}}_n^{\phantom{(}}\rangle \langle {\widehat{a}}_r^{†}{\widehat{a}}_s^{\phantom{(}}\rangle .\hfill \end{array}$$

(22)

One gets a closed set of differential equations for the GFs in the Bogoliubov Back-reaction beyond mean-field approximation. For example, for the greater GF, this reads:

$$\begin{array}{cc}i{\partial }_{\tau }{G}_{j,k}^{>}(t+\tau ;t)\stackrel{\mathrm{BBR}}{\approx }\hfill & -J\left(G_{j+1,k}^{>}(t+\tau ;t)+G_{j-1,k}^{>}(t+\tau ;t)\right)+\left(U{n}_j(t+\tau )-i\frac{{\gamma }_j}{2}\right)G_{j,k}^{>}(t+\tau ;t)\hfill \\ & +U{F}_{j,j,j;k}^{>}(t+\tau ;t),\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill i{\partial }_{\tau }{F}_{j,m,n;k}^{>}(t+\tau ;t)\stackrel{\mathrm{BBR}}{\approx }& +J(\phantom{+}{F}_{j+1,m,n;k}^{>}(t+\tau ;t)+F_{j-1,m,n;k}^{>}(t+\tau ;t)-F_{j,m+1,n;k}^{>}(t+\tau ;t)\hfill \\ & \phantom{+J(}-F_{j,m-1,n;k}^{>}(t+\tau ;t)-F_{j,m,n+1;k}^{>}(t+\tau ;t)-F_{j,m,n-1;k}^{>}(t+\tau ;t))\hfill \\ & -U(\phantom{+}{F}_{j,j,n;k}^{>}(t+\tau ;t){\sigma }_{jm}(t+\tau )-({\Delta }_{j,m,m,m}+{\Delta }_{j,m,n,n})(t+\tau )i{G}_{n,k}^{>}(t+\tau ;t)\hfill \\ & \phantom{-U(}+F_{j,m,n;k}^{>}(t+\tau ;t)(n_j-n_m-n_n)(t+\tau )\hfill \\ & \phantom{-U(}+{\delta }_{m,n}(F_{j,m,n;k}^{>}(t+\tau ;t)+{\sigma }_{j,m}{G}_{n,k}^{>}(t+\tau ;t))).\hfill \end{array}$$

(24)

Here, we defined

$$i{F}_{j,m,n;k}^{>}(t^{\prime };t)=\langle \Delta ({\widehat{a}}_j^{†}{\widehat{a}}_m^{\phantom{(}})(t^{\prime })\Delta ({\widehat{a}}_n^{\phantom{(}}(t^{\prime }){\widehat{a}}_k^{†}(t))\rangle =\langle ({\widehat{a}}_j^{†}{\widehat{a}}_m^{\phantom{(}}{\widehat{a}}_n^{\phantom{(}})(t^{\prime }){\widehat{a}}_k^{†}(t)\rangle -{\sigma }_{j,m}(t^{\prime })i{G}_{n,k}^{>}(t^{\prime };t).$$

(25)

The SPDM and the covariances are given in the BBR approximation by (c.f. [7]):

$$\begin{array}{cc}\hfill i& {\partial }_t{\sigma }_{jk}\stackrel{\mathrm{BBR}}{\approx }-J\left({\sigma }_{j,k+1}+{\sigma }_{j,k-1}-{\sigma }_{j+1,k}-{\sigma }_{j-1,k}\right)+U\left({\Delta }_{jkkk}+{\sigma }_{jk}{\sigma }_{kk}-{\Delta }_{jjjk}-{\sigma }_{jj}{\sigma }_{jk}\right)-i\frac{{\gamma }_j+{\gamma }_k}{2}{\sigma }_{j,k},\hfill \end{array}$$

(26)

$$\begin{array}{c}i{\partial }_t{\Delta }_{jmkn}\stackrel{\mathrm{BBR}}{\approx }\hfill \\ -J\left({\Delta }_{j,m,k,n+1}+{\Delta }_{j,m,k,n-1}+{\Delta }_{j,m+1,k,n}+{\Delta }_{j,m-1,k,n}-{\Delta }_{j,m,k+1,n}-{\Delta }_{j,m,k-1,n}-{\Delta }_{j+1,m,k,n}-{\Delta }_{j-1,m,k,n}\right)\hfill \\ +U\left({\sigma }_{jm}({\Delta }_{mmkn}-{\Delta }_{jjkn})+{\sigma }_{kn}({\Delta }_{jmnn}-{\Delta }_{jmkk})+{\Delta }_{jmkn}(-{\sigma }_{jj}+{\sigma }_{mm}-{\sigma }_{kk}+{\sigma }_{nn})\right)\hfill \\ -i\frac{{\gamma }_j+{\gamma }_m+{\gamma }_k+{\gamma }_n}{2}{\Delta }_{jmkn}.\hfill \end{array}$$

(27)

However, although being more accurate in the low population regime, no qualitative difference is expected with respect to the mean-field version, at least in the limit of large filling factors $n_j\gg 1$.

6. Density–Density Correlation Function

Besides improving the mean-field precision, BBR also provides a way to compute non-trivial density–density correlation functions. The latter indicate, for instance, possible temporal bunching or anti-bunching effects [7,15]. Indeed, at the mean-field level, one just gets $\langle ({\widehat{a}}_j^{†}{\widehat{a}}_m^{\phantom{(}})(t^{\prime })({\widehat{a}}_k^{†}{\widehat{a}}_n^{\phantom{(}})(t)\rangle \stackrel{\mathrm{MF}}{\approx }{\sigma }_{jm}(t^{\prime }){\sigma }_{kn}(t)$. Instead, by defining the density–density correlation function as follows:

$$C_{j,m;k,n}(t^{\prime };t)=-i\langle \Delta ({\widehat{a}}_j^{†}{\widehat{a}}_m^{\phantom{(}})(t^{\prime })\Delta ({\widehat{a}}_k^{†}{\widehat{a}}_n^{\phantom{(}})(t)\rangle =\langle ({\widehat{a}}_j^{†}{\widehat{a}}_m^{\phantom{(}})(t^{\prime })({\widehat{a}}_k^{†}{\widehat{a}}_n^{\phantom{(}})(t)\rangle -{\sigma }_{jm}(t^{\prime }){\sigma }_{kn}(t)$$

(28)

and truncating sextic moments according to Equation (22) in the equation of motion similar to Equation (12) that it satisfies, one gets:

$$\begin{array}{c}i{\partial }_{\tau }{C}_{j,m;k,n}(t+\tau ;t)\stackrel{\mathrm{BBR}}{\approx }\hfill \\ +J\left(C_{j+1,m;k,n}(t+\tau ;t)+C_{j-1,m;k,n}(t+\tau ;t)-C_{j,m+1;k,n}(t+\tau ;t)-C_{j,m-1;k,n}(t+\tau ;t)\right)\hfill \\ -U\left((n_j-n_m)(t+\tau )C_{j,m;k,n}(t+\tau ;t)+{\sigma }_{jm}(C_{j,j;k,n}(t+\tau ;t)-C_{m,m;k,n}(t+\tau ;t))\right)\hfill \\ -i\frac{{\gamma }_j+{\gamma }_m}{2}{C}_{j,m;k,n}(t+\tau ;t),\hfill \end{array}$$

(29)

which has to be evolved together with Equations (26) and (27). Again, if one wants to compute the density–density correlation function with the first time argument evaluated at an earlier time than the second, one can simply use

$$C_{j,m;k,n}(t;t+\tau )=-{\left(C_{n,k;m,j}(t;t+\tau )\right)}^{*}.$$

(30)

This shows that the systematic expansions in order to approximate the temporal correlation functions of operators proposed here are indeed very useful to compute physical relevant quantities in open quantum many-body systems.

7. Conclusions

Usually for a quantum many-body problem, the Schrödinger equation, or in our case the master equation Equation (2), is not analytically solvable. One possible way to access the non-equilibrium dynamics of such systems is with the help of Green functions. The two-time correlators, such as the time-dependent first order coherences from Equations (6) and (7), or the density–density correlation functions Equation (28), allow one to characterize the out-of-equilibrium dynamics of many-body interacting quantum systems. We have presented a method on how to systematically compute a hierarchy of approximations for two-time correlation functions, in principle of arbitrary operators, which applies to a wide class of interacting systems whose dissipation is described by a master equation of the form of Equation (2). Based on our method, we can treat particle interactions on a mean-field level (see Section 4) and also one-order beyond (see Section 5). Higher-order expansions are possible, at the cost of very lengthy formulae, which at some point will become difficult to deal with even numerically (see e.g., [28,29] for just time-equal correlators).

With our method, we computed the spectral weight function in Section 4.1 as an example of a physical observable. In principle, we can compute the temporal evolution and the correlations of many interesting quantities for experiments, such as the onsite populations or the density–density correlations between two lattice sites. For about ten years, state-of-the-art experiments with ultracold atoms can access information on static correlation functions with high spatial resolution (see e.g., [30,31,32,33,34,35,36]). Time-dependent two-point correlation functions can be measured too, see e.g., [37], which would allow for direct applications of our developed theoretical approach.

Our problem is defined by a time-local master equation, which means that the coupling to the environment, in our case a zero-temperature sink, is supposed to be sufficiently weak for justifying the Markovian assumption underlying Equation (2). Strong coupling to the environment, such as that present in lead-to-lead transport across solid-state samples (see e.g., [18,19]) remains an open problem since then we may not apply the quantum regression theorem valid for time-local master equations. For tedious extensions of the theorem to non Markovian setups, the reader may consult e.g., [38,39,40,41]. An alternative approach would be then to fully include the leads into the treatment and use a diagrammatic expansion of the Green functions computed for the full lead-system-lead system, such as done for bosons in steady state in Refs. [42,43]. Needless to say, such an approach as just mentioned is very hard for many-body quantum systems whose main central part, without the leads, is itself non-integrable, such as our many-body Bose–Hubbard model for more than two wells [44,45].