Modeling Particle Loss in Open Systems using Keldysh Path Integral and Second Order Cumulant Expansion

For open quantum systems, integration of the bath degrees of freedom using the second order cumulant expansion in the Keldysh path integral provides an alternative derivation of the effective action for systems coupled to general baths. The baths can be interacting and not necessarily Markovian. Using this method in the Markovian limit, we compute the particle loss dynamics in various models of ultra-cold atomic gases including a one-dimensional Bose-Hubbard model with two-particle losses and a multi-component Fermi gas with interactions tuned by an optical Feshbach resonance. We explicitly demonstrate that the limit of strong two-body losses can be treated by formulating an indirect loss scheme to describe the bath-system coupling. The particle-loss dynamics thus obtained is valid at all temperatures. For the one-dimensional Bose-Hubbard model, we compare it to solutions of the phenomenological rate equations. The latter are shown to be accurate at high temperatures.


I. INTRODUCTION
Open quantum systems constitute one of the most interesting challenges of the field, both from a fundamental point of view and of course in connection with applications or experimental realizations.This class of problems offers an interesting interplay between coherent Hamiltonian dynamics and incoherent, dissipative dynamics emerging from coupling to the environment [1][2][3].From a fundamental point of view, after being considered as a nuisance due to the destruction of coherence, dissipative processes are now regarded as a resource allowing to engineer novel quantum states of matter [4][5][6][7] or to assist for quantum transport properties [8][9][10].Among the various types of dissipative processes, particle losses have played a special role, and recent experiments in cold atomic systems have allowed to control them in an exquisite manner.This is for example the case of losses in weakly interacting Bose gases [11,12] or fermionic systems where local losses are realized [13,14].More recently, cavities have also provided an interesting playground for this kind of physics [15][16][17].
Theoretically studying dissipative phenomena is a considerable challenge, and several approaches have been used in order to deal with the coupling to the environment, which is often modelled as a bath.One common approach resorts to non-Hermitian Hamiltonians [18][19][20][21].The resulting (non-Hermitian) models can be analyzed using e.g. the powerful tools of integrability [22,23] as well as various field-theoretic and numerical techniques [18,19] which include bosonization and the renormalization group [20,23,24].However, because such treatments often neglect or at most treat in an approximate way the so-called quantum-jump term of the Lindblad master equation, it is difficult to assess how reliably they can describe the dynamics of the system observables.Other approaches deal with the full Lindblad master equation [25][26][27][28][29] assuming the bath is Markovian [1,30].In this work, we are concerned with another generic method to deal with such out of equilibrium systems, namely the Keldysh formalism and its path integral formulation [31,32].
Approaches based on the path integral have a long history beginning with the pioneering work of Feynman and Vernon [33].They have been often used in the context of quantum dissipation where the coupling to ohmic, subohmic or superohmic baths generates an effective longrange interaction in imaginary time [34][35][36][37].In a more general context, they have also been used to describe non-equilibrium dynamics in the presence of coupling to baths [31,38], time dependent environments [10,39] and single particle losses [32,[40][41][42].
Despite these important developments, several phenomena linked to the coupling to the environment still remain elusive.In particular, in ultracold atomic systems the dynamics in the presence particle losses involve a priori processes beyond single-particle losses and include two-or three-particle losses as well.Quite generally, the description of particle losses has assumed the validity of phenomenological rate equations of the form: where n is the particle density and γ 1 and γ 2 are the oneand two-particle loss rates, respectively.However, the correct functional form and range of validity of these phenomenological equations and how to account for manyparticle effects are still not fully understood.For ultracold atoms near Feshbach resonances, a theoretical description of the two-body loss rate, γ 2 , was developed based on S-matrix calculation in Ref. [43] and some properties of three-body losses could be determined from Bethe ansatz solutions [44].Braaten et al. [45] combined effective field theory with the Lindblad master equation to obtain a universal relation for the two-atom inelastic loss rate for two-species Fermi gases.
In this work, in order to go beyond the phenomenological description of particle losses provided by Eq. ( 1) and be able to account for many-particle effects, we develop a path integral formalism that relies on the second order cumulant expansion.The formalism can in principle also deal with interacting as well as non-Markovian baths.However, here we demonstrate it by studying the Markovian dynamics of many-particle systems in the presence of one-and two-particle losses and briefly discuss how to go beyond the Markovian limit in the Appendix C. In the absence of interactions in the bath, the resulting approach is exact for one-particle losses and approximately valid for two-particle loss problems at weak coupling.In order to deal with the strong coupling regime, we also introduce an indirect two-body loss scheme.We discuss two experimentally relevant examples of this indirect-loss scheme (a lossy one-dimensional Bose-Hubbard model and a multi-component Fermi gas near an optical Feshbach resonance) and explicitly derive the loss rate equations.In both cases, we find their functional form deviates from Eq. ( 1).In the case of the lossy Bose-Hubbard model, we explicitly show that, at low-temperatures where quantum coherence is important, predictions of the microscopic theory deviate substantially from those obtained using the phenomenological equation (1).
The rest of this article is organized as follows.In Sec.II we describe the derivation of the Keldysh path integral using the 2nd order cumulant expansion.We first illustrate the method with a model consisting of a singlemode (fermonic or bosonic) oscillator coupled to a bath to which particles can be lost (cf.Sec.II A).The general formalism is discussed in Sec.II B. However, since only the results of subsection II A will be used in the rest of the article, Sec.II B can be skipped on a first reading.Sec.III briefly describes how the results of Sec.II A are applied to describe a quantum gas with one-body losses.Sections IV and V deal with the applications of the formalism to two models for relevant for the physics of ultracold atomic gases with two-body losses: Sec.IV is concerned with a lossy Bose-Hubbard model in one dimension, whereas Sec.V deals with the model of a multicomponent Fermi gas near optical Feshbach resonance.Some technical details and useful results are described to the Appendices.

A. Single-mode case
To begin with, let us illustrate the method by considering a single bosonic or fermionic mode (system A) coupled to a bath (B) by means of a quadratic Hamilto-nian: Here a, a † describe a bosonic (fermionic) mode in the A and b α , b † α a set of bosonic (fermionic) modes in B labelled by a continuum index α.We shall assume that A and bath B are in equilibrium (not necessarily with each other) at t = −∞ and the interaction between them H AB (t) is switched according to a protocol described by the function f (t).Both A and B can be interacting systems with e.g.
The Keldysh generating functional [31] for the system introduced above can be written as follows: The Keldysh action is S = S A + S B + S V + S AB , where In the above expressions, C is the Keldysh contour, which runs from t = −∞ to t = +∞ and back to t = −∞ [31]; V , V are sources that couple the system degrees of freedom: Any n-point correlation of the system A can be obtained by conveniently taking functional derivatives of Z[ V , V ] with respect to V and V .The generating functional is normalized so that Z[ V = 0, V = 0] = 1 since in the absence of external sources it merely describes the unitary evolution of the initial state from t = −∞ to +∞ and back to t = −∞.Starting from the above functional integral representation, we can define the Feynman-Vernon influence functional [33] F[ā, a] as the result of formally integrating out the bath degrees of freedom, i.e.
It is often not possible to obtain F in a closed form and therefore we have to resort to approximations.For a general (e.g.interacting) bath that is weakly coupled to the system the 2nd order cumulant expansion provides a good starting point to capture the dissipative dynamics induced by the bath.Furthermore, it is exact if the bath Hamiltonian H B and the coupling H AB are quadratic and linear in the fields bα , b α , respectively.Using the cumulant expansion to second order (see Appendix B), we approximate: where ⟨• • • ⟩ B stands for average over the bath degrees of freedom.Taking ⟨S AB ⟩ B = 0 after assuming that the bath conserves the number of bath excitations in its initial state, we obtain the following dissipative effective action: In the last expression, we have split the integrals over the Keldysh contour C.After rewriting them as a single integral where t 1 , t 2 run from −∞ to +∞, we have introduced the subindices m, n = ± to denote on which branch of C the time argument of the fields a, b, etc lies.We have also introduced the functions v mn (t 1 , t 2 ) = f (t 1 )f (t 2 )g mn (t 1 − t 2 ), where In the above expression s m=± = ±1 and z = +1 (z = −1) for bosons (fermions).In the Appendix, within Markovian limit, we show that g mn (t 1 − t 2 ) ∼ δ(t 1 − t 2 ).Thus, we arrive at the following expression for the effective couplings (see Appendix C for more detail), where ⟨g⟩ is an average system-bath coupling constant.Upon denoting γ(t) = ν 0 |⟨g⟩f (t)| 2 for the loss rate of particles to the bath, we obtain the following result: This is the dissipative part of the action characteristic of a Markovian bath.It can also be obtained from the evolution of the density matrix according to the Lindblad master equation (see Appendix A, which is based on Ref. [32] and references therein).Finally, let us compute the particle loss in the system A caused by switching on the coupling to the bath B at t = 0 for an infinitesimal time δt.This calculation can be carried out by the path integral version of timedependent perturbation theory, i.e. by perturbatively expanding the effective dissipative action to leading order in L, which yields: = where we have set γ(t) = γθ(t) and denoted n a = ⟨a † a⟩, which is the occupation in the initial state (i.e. for t < 0).Assuming the system A is non-interacting, we have Hence, Setting t = δt ≫ D −1 (where D −1 is the characteristic response time of the bath, see Appendix C), we arrive at the following one-body loss rate equation: In the right hand-side, with accuracy O(γδt), we have replaced n a by n a (δt).

B. General case
In this section, we generalize the above results.Since only the results obtained in previous section will be necessary the applications discussed in this work, this section can be entirely skipped on a first reading.The starting point is again the Hamiltonian describing the unitary dynamics of a system (A) coupled to a bath (B).
where H A (t), H B , and H AB (t) denote the Hamiltonians of the system, bath and their coupling, respectively.We have assumed that, in general, the system Hamiltonian and its coupling to the bath can be explicitly timedependent.We shall consider a rather general form of the coupling between system and bath: where are products of arbitrary number operators (with N ̸ = N ′ and M ̸ = M ′ in general) acting either on the system or the bath; are the quantum numbers carried by those operators; g qp (t) are the set of systembath couplings.The system-bath coupling is switched on according to a certain protocol that determines the explicit time dependence of the g qp (t).
In general, we are not able to obtain the Feynman-Vernon influence functional F exactly and here we will resort to a second order cumulant expansion: Explicitly, for the system-bath coupling introduced in Eq. ( 27) the first order correction takes the form: Note that this term does not describe dissipation and only modifies the unitary evolution of the system A.
Thus, it can be conveniently absorbed into H S by defining the operators B p and B † p in Eq. ( 27) to have zero averages, i.e. ⟨B p ⟩ B = ⟨B † p ⟩ B = 0.In addition, this is automatically fulfilled if B and B † change the number of particle/excitations in the bath but the bath Hamiltonian H B and its initial state conserve this number.Therefore, in what follows, we set ⟨S AB ⟩ B = 0 and do not discuss it any further.
The second order correction can be brought to the following form in terms of bath correlators: where In the above equation, the following notation has been introduced (X, Y = B, B): for the bath two-point correlation functions of operators B, B † .In the context of the Keldysh path integral each one of the above correlation functions becomes a 2 × 2 matrix in the superindices m, n = ± after expanding the integrals over C so that t 1 , t 2 run from −∞ to +∞.The superindices are inherited by the couplings u, v, ū, v introduced in Eq. (36), which also become the 2 × 2 matrices u mn q1q2 (t 1 , t 2 ), v mn q1q2 (t 1 , t 2 ), etc.
Explicit consideration of the time-dependence of the coupling matrices u mn q1q2 (t 1 , t 2 ), v mn q1q2 (t 1 , t 2 ), etc, for a given bath may allow to identify regimes where the Markovian approximation applies.This is typically the case when the response of the bath is much faster than the characteristic time scale of the system dynamics.Thus, in the Markovian regime we can assume that u mn q1q2 (t 1 , t 2 ) ∝ δ(t 1 − t 2 ), etc (although some of them may also vanish as it was the case in the previous example).This provides an additional simplification of the second order term of the cumulant expansion and leads to the following result: where S N is the "normal" part: and L A is the "anomalous" part In the above expressions, we have introduced the following system-bath coupling matrices: where z = −1 (z = +1) if the operator products A q , A † q in H SB have fermionic (bosonic) statistics, and s + = +1 (s − = −1) .We can classify the different terms according to the type of time arguments of the system degrees of freedom A q (t), Āq (t).The first two terms on the right-hand site contain A qm (t), Āqn (t) whose time arguments lie on different branches of the Keldysh contour C and therefore m ̸ = n.These terms correspond to the so-called quantum jump terms of the Lindblad master equation (see Sec. III and Appendix A).The remaining terms contain A qm (t), Āqn (t) with time arguments lying on the same branch of C, i.e. with m = n.Such terms contribute to the anti-commutator terms of the Lindblad master equation in the operator language (see Appendix A) and give rise to the anti-Hermitian part of the effective non-hermitian Hamiltonian in the operator language [1].
Finally, although the last few expressions above have been derived under the assumptions of Markovianity of the bath, we want to emphasize that the approach used here is not limited to the Markovian regime and it can be used as a starting point to include effects beyond Markovianity.The latter are outside the range of applicability of the Lindblad master equation or its path integral formulation as introduced in Ref. [32].In the following sections, we shall consider a number of applications to particle loss and show that, although derived using the cumulant expansion up to second order which may appear to be only valid for a weak system-bath coupling, it is possible reformulate the models to describe the limit of very strong two-body losses exactly.

III. ONE-BODY LOSS
Before considering systems with two-body losses, it is interesting to generalize the results of Sec.II A to describe an ultracold gas coupled to a bath to which it can loose one particle at a time.This section largely relies on the results obtained in Sec.II A. We generalize the Hamiltonian to a uniform gas and therefore the fields carry a momentum index k.The Hamiltonian reads: where the fields a k and b k,α are either fermionic or bosonic.In ultra-cold atomic gases, the couplings g α , g α are often not known from first principles.Instead, what is measured is the loss rate of particles.Following the same steps as in Sect.II A while keeping track the momentum index k, we arrive at the following effective action in the Markovian limit (m, n = ±): where we have introduced the following short-hand notations: , and zero otherwise.
Like in Sect.II A, we can obtain the rate of particle loss by assuming the coupling to the bath is switched on at t = 0 (i.e.f (t) = θ(t)) and computing the change in the total particle number density using perturbation theory in L: Following the same steps as in Sec.II A, we compute the leading order change of n k (t) = ⟨a † k a k ⟩.Thus, we obtain rate equations for the momentum distribution dn k (t)/dt = −γn k (t), and hence the rate of change of the particle density follows: Note that the rate of change is proportional to the density, which is characteristic of the one-body loss process.

IV. LOSSY 1D BOSE HUBBARD MODEL
Next, we discuss how to describe two-body losses by coupling a system to a bath.We first consider an ultracold gas in a deep optical lattice which can be described by the one-dimensional Bose-Hubbard model [46].A laser is applied that photo-associates atoms in doubly occupied lattice sites (doublons) into molecules.The molecules are quickly lost from the trap, resulting in a loss of two bosons (see Fig. 1).In the limit where the loss of the photoassociated molecule is very fast, the double occupancy is strongly suppressed.In other words, in the presence of this coupling, the states containing doublons are rapidly projected out.However, virtual transitions to states containing doublons can still have an effect on the system dynamics.This kind of loss dynamics was experimentally studied in [47,48] and theoretically described using an approach based on an effective master equation derived in [49], which yields a phenomenological rate equation like (1).In the following, we provide a microscopic theory for the loss dynamics and compare it to the phenomenological loss equation.
A convenient way to describe the dynamics of the Bose-Hubbard model in a subspace containing no doublons relies on the Jordan-Wigner transformation [46] that relates hard-core bosons to fermions: where n j = a † j a j = c † j c j = 0, 1 measures the occupation of site j.The transformation holds true provided (a † j ) 2 = a 2 j = 0 (i.e.no doublons) in the relevant Fock subspace.In this subspace, the kinetic energy of the J U FIG. 1. Scheme of the lossy 1D Bose-Hubbard model.J is the nearest-neighbor hopping and U is the onsite interaction.The loss is parametrized by γ is the one-body loss rate of the doublons, i.e. doubly occupied sites.
hardcore bosons can be written in terms of the Jordan-Wigner fermions (see Fig. 1): In addition, we note that the original hopping operator of bosons −J i a † j a j+1 + h.c. also allows for transitions that create virtual doublons.In order to allow for such processes, we introduce a doublon field on each site d † i that is coupled to the (hardcore) bosons by means of −J j (d j + d j+1 ) a † j a † j+1 + h.c.. Following the Jordan-Wigner transformation, this coupling becomes: Finally, we note that the doublon has an excitation energy equal to U , and therefore, This rather heuristic derivation of the Hamiltonian in the limit where the doublons are suppressed is confirmed below by showing that it is a convenient Hubbard-Stranovich decoupling of the effective interaction generated in the strongly interacting limit of the onedimensional Bose-Hubbard model i.e. for U ≫ J [50].The generating functional for the above model is: where As mentioned above, the doublon field will be treated as Hubbard-Stranotovich field and therefore its "action" does not contain a time-derivative term i d∂ t d: The calculations to be described below can be also carried out including such derivative term, which is not important in the limit where U ≫ J (see next section for a discussion where a similar situation is encountered).We further assume that the system is coupled to a bath that removes the doublons.This coupling is described by the following term in the Keldysh action: where the bath modes have the following quadratic action: We integrate out the bath following the same steps as in Sect.II A, which in the Markovian limit yields: Combining this term with the action for the doublon, the following effective action is obtained: Finally, we integrate out the doublon field by making the following change of integration variables in the functional integral: 65) where we have denoted h j,m = −J(c j+1,m c j,m + c j,m c j−1,m ), hj,m = −J(c j,m cj+1,m + cj−1,m cj,m ), and The resulting integral over d′ j , d ′ j is gaussian, and yields a constant prefactor to the generating functional in Eq. ( 55).In addition, there is an exponential factor with the following effective action in the exponent: dt hj,m (t) σ 3 G(t)σ 3 mn h j,n (t).(68) Note that in the limit where the coupling to the bath vanishes and γ(t) = 0, G mn (t) = −σ 3 mn /U and we obtain which is the Keldysh action for the effective Hamiltonian obtained using strong coupling perturbation theory in Ref. [50] for the 1D Bose-Hubbard model in the limit where U ≫ J and in the subspace with no doublons.Next, we switch on the coupling to the bath so that γ(t) = γθ(t) and compute the particle-loss rate.It is convenient to work in the Bloch wave basis where c k = j e −ikxj c j / √ M , M being the number of lattice sites, x j = j, k = 2πl/M , l = −M/2 + 1, . . ., M/2 (assuming periodic boundary conditions and M to be even).Thus, the full effective action S eff = S c + S ′ eff reads: where ϵ k = −J cos k is the single-particle dispersion, Ũ (p, k, q) = U pkq (t) − iσ 3 mn Γ pkq (t) and Γ pkq (t) are given by: In the limit where J ≪ max{U, γ(t)}, both U pkq (t) and Γ pkq (t) are perturbatively small.Using perturbation theory to leading order in U pkq (t) and Γ pkq (t) we obtain the following loss rate equation for the distribution function of Jordan-Wigner fermions (see Appendix D for further details), where C kp = sin p−k 2 cos p+k 2 and the effective loss rate is Integrating over p yields the following rate equation for the lattice filling n A (t) = N A (t)/M : The form factor C 2 kp was not present in the approach used in Ref. 49, which neglected inter-site correlations.The expression containing the form factor was later obtained by solving the Lindblad master equation using an approximation termed as time-depependent generalized Gibbs ensemble [51,52].As illustrated below, the form factor turns out very important when the filling of the lattice is low and at low initial temperatures.The phenomenological loss rate equation and the two-body loss coefficient γ 2 (cf.Eq. 1) [47][48][49] can be obtained in the high-temperature limit where the Fermi-Dirac distribution of the Jordan-Wigner fermions approaches n A (t) and is independent of the momentum.Therefore, with two-body loss coefficient γ 2 = γ T = γ eff /4.Fig. 2 compares the solution of the phenomenological rate equation with the numerical solution of Eq. ( 75) obtained from the microscopic Keldysh action at different initial temperatures determined by the initial momentum distribution of the Jordan-Wigner fermions.Indeed, since after t = 0 the system is out of equilibrium, a global temperature can no longer be defined.At a small values of the initial lattice filling (e.g.n A (t = 0) = 0.25), the results from the phenomenological rate equation strongly deviate from the predictions of the microscopic theory, Eq. ( 75).The difference between the two rate equations becomes smaller as temperature increases becasue, as mentioned above, the phenomenological rate equation is recovered from Eq. (75) in the high temperature limit where T ≫ J.However, at a higher initial lattice filling (n A (t = 0) = 0.75), the error incurred by using the phenomenological rate equation becomes smaller at all the studied temperatures.Thus, the phenomenological rate equation only applies at high temperatures or high lattice fillings (i.e.n A (t = 0) ≃ 1).The latter are indeed the conditions of the previous experiments [47,48].Our formalism thus allows to access other regimes of lattice fillings and temperature, which can be explored in future experiments.
In addition, a closer examination of the numerical solution of Eq. (75) shows that the situation is indeed more complex than what can be naively inferred from the above discussion of the validity of the phenomenological rate equation.In order to see from where the complexity emerges, we have plotted the evolution of the distribution function of Jordan-Wigner fermions in Fig. 3 for an initially half-filled lattice and different values of the effective loss rate γ eff (right panels).In the same figure, we also show the evolution of the lattice filling and kinetic energy for different initial temperatures (left panels).Note that for the lattice filling shown on the top left panel, the particle loss dynamics is largely independent of the initial temperature (note the log-log scale).On the other hand, the dynamics of the average kinetic energy (shown in a linear scale) does indeed depend on the initial temperature.As it can be seen, the change in kinetic energy induced by the two-body loss is larger in systems with lower initial temperature.
The dependence on the initial temperature of the kinetic energy, which is the first moment of the momentum distribution, is an indication that the two-body loss drives the system into a non-equilibrium state.To confirm this observation, we focus on the evolution of the full momentum distribution of the Jordan-Wigner fermions, which is shown on the right panels for two values of the initial temperature (k B T = J and k B T = J/10).Although the initial distribution is assumed to be thermal, it can be seen that under the two-body loss in both cases it rapidly evolves into a non-thermal distribution.Note that the depletion is most effective (especially at low temperatures) for k near ±π/2.This is because the form factor C 2 kp in Eq. ( 77) is maximum for k = −p = ±π/2, which corresponds to losses of doublons with total zero momentum.The latter are created from "Cooper pairs" ∼ c k c −k of Wigner fermions.In conclusion, even in cases where evolution of the particle filling stays rather close to the results obtained from the phenomenological rate equation (78) down to low temperatures, the system is indeed far from equilibrium as revealed by close examination of other observables like the kinetic energy.

V. LOSSES IN OPTICAL FESHBACH RESONANCE
In this section we describe a multi-component mixture of alkaline-earth atoms with emergent SU(N ) symmetry [53,54] near an optical Feshbach resonance (OFR) 4. Scheme of an optical Feshbach resonance: The system consists of N species of spin-F (N = 2F + 1) alkaline-earth atoms in the ground state via an s-wave potential.Using a laser colliding pairs of atoms in the ground-state are coupled to pairs consisting of one ground-state and one excited-state atom in a molecular state.This state has a onebody coupling parametrized by gα to a bath to which it can be lost or gained.
by means of the following model: Fermions in the ground state 1 S 0 interact via a weak swave potential preserving SU (N ) symmetry.The OFR is described by a SU(N ) symmetry breaking scattering channel which involves an intermediate bosonic molecular state where one of the colliding atoms is in an optically coupled excited state, i.e. 1 S 0 + 3 P 0 .The coupling g J (q) = g J (−q) is the matrix element of the laser-induced transition between two-particles in the ground state and the molecular excited state, i.e. ⟨ 1 S 0 3 P 0 |V las | 1 S 0 1 S 0 ⟩ .In the calculations below, we assume that g J ≪ ∆ m , γ, where ∆ m is the detuning from the excited state and γ is the one-body loss rate of the excited molecular state.In this limit, the loss of the molecular state is due to spontaneous emission.The contribution of stimulated emissions is relatively small.The two-body loss of the system particles is described as the one-body loss of the intermediate molecular state, which is coupled to a bath via H aB whose eigenmodes are the bosonic operators b † q , b q .Since the coupling is linear and the bath is described by a quadratic Hamiltonian, the second-order cumulant expansion is exact and applies even in the limit where the loss rate γ is large.Below we write the Keldysh action after integrating out the bath field following the same steps than in Sec.II.This yields: where In Eq. ( 86), in order to lighten the equations, we have adopted the convention that summation over repeated indices m, n = ± is implied.Next, we integrate out the molecular field, āq,JM , a q,JM using the cumulant expansion.
To this end, we need to obtain the Green's function for the molecular fields, i.e. the inverse of the matrix G −1 0 (t, t ′ ).In order to take into account the distribution function of the molecules correctly, it is necessary [31] to take a step back and work with the discrete version of the path integral.However, the presence of the dissipative terms make the inversion of the matrix cumbersome.One way around this difficulty is to realize that in the parameter regime of interest here, i.e. for large γ(t) and/or large detuning ∆ m , it is possible to neglect the time derivative part of G −1 q (t, t ′ ), which effectively amounts to a Markovian approximation when the molecular field is regarded as a bath itself.In order to understand this, let us neglect the time dependence of the coupling to the bath γ for the time being.Thus, the matrix G −1 q becomes a function of t − t ′ and it can be Fourier transformed, which yields the following matrix in Keldysh space: If |ω| ≪ min{∆ m , γ} we can neglect the frequency dependence of G −1 q (ω).Furthermore, assuming a vanishing number of molecules in the initial state, the correlation functions that determine the second order term in the expansion are given by the inverse of above matrix with ω = 0.In Eq. (87) this amounts to neglecting the term involving the time derivative, i∂ t , which yields Markovian correlations for the molecular field of the form: Using the above expression and Eq. ( 39), we obtain the following effective Keldysh action: where U eff (p, k, q, t) = U 0 + δU (p, k, q, t) is the renormalized interaction: and is the effective two-body loss rate in the limit of strong spontaneous loss on the intermediate molecular states [43,[55][56][57][58][59][60]].In the above expressions ∆ p = ϵ a p + ∆ is the energy of the excited molecular state with total momentum p.Notice that δU is an SU (N )-symmetry breaking interaction and both quantities are perturbatively small in the limit where g J ≪ min{∆ m , γ} of interest here.Hence, perturbation theory to leading order in γ ′ yields the following rate equation for an quantum degenerate gas: where n p (t) is the instantaneous momentum distribution and n c = p n p (t)/Ω is the fermions density.From this result, the phenomenological loss coefficient for a thermal gas can be obtained by replacing the loss-coefficient with its thermal average [56,61,62], Here f M (p, T ) = (2πmk B T ) −3/2 exp(−p 2 /2mk B T ) denotes the Maxwell distribution at temperature T for particles with mass m normalized to unity.Using the average coefficient, the rate equation can be approximated by which is the phenomenological two-body loss rate equation, Eq. ( 1) with γ 1 = 0 and γ 2 = 2γ T for a thermal gas [56,57,59,61,62].

VI. CONCLUSION
In this work, we have discussed the derivation of the Keldysh path integral for open quantum systems using the 2nd order cumulant expansion.Although we have focused on Markovian baths, the method is not limited to the latter and can be extended to describe effects beyond Markovianity.It also does not require the bath to be non-interacting or the bath coupling to be of a particular form.Formally, it requires that the coupling to the bath is weak enough to be accurately treated using second order perturbation theory.However, as we have shown above, when describing two-particle losses, the system-bath coupling can be sometimes conveniently reformulated and a strong loss regime can be also described.
Turning to models relevant for ultra-cold atomic gases, we have studied models of one and two-body losses.Thus, we have shown how two-body losses caused by photo-association of doublons in the one dimensional Bose-Hubbard model can be described within the pathintegral formulation allowing us to obtain a microscopic loss rate equation.The latter has been compared with a previously derived phenomenological loss equation.The microspically derived rate equation shows that the phenomenological equation is mostly accurate at high temperatures and/or lattice fillings close to unity.However, we have shown (see Sect.IV) that even in cases where the phenomenological rate equation appears to be sufficiently accurate, the implicit assumption that the system remains in a thermal equilibrium characterized by a temperature T can be incorrect.This has important implications for the calculation of other physical quantities such as the kinetic energy, for which our theory, which properly handles such deviations from equilibrium, is necessary.
Finally, in Sec.V we have applied the formalism to a model describing an optical Feshbach resonance in multicomponent mixture of ultracold alkaline-earth fermions.We have also shown that the phenomenological rate equation is expected to apply to the high temperature regime.Although we have not fully explored the low temperature regime yet, using the lessons learned with the much simpler one-dimensional Hubbard model, in presence of an OFR we expect that if quantum coherence becomes important, deviations from the phenomenological approach will appear.
Finally, let us mention that, in this work, when dealing with the interactions between the Jordan-Wigner fermions in the lossy one-dimensional Bose-Hubbard model we have used perturbation theory.This is justified because the latter are weak [50] and the studied temperatures are relatively high.However, in future work it will interesting to revisit this problem in order to account for the effect of the interactions in the system, which in one dimension can be done non-perturbatively using bosonization [63,64].Other possible extensions of this work are, as pointed out above, studying effects beyond Markovianity and the effect of strong interactions/correlations in the bath.
oscillation is an artifact of the hard cutoff).Thus, if the dynamics of the system is characterized by frequencies much smaller than D, we can effectively replace the above function by a Dirac delta-function δ(t 1 −t 2 ), which yields, In order to obtain the above expressions, we have used θ(t)δ(t) = θ(0)δ(t) = δ(t) and θ(t)δ(t) = θ(0)δ(t) = 0, as required by the discrete version of the path integral [31].The θ(t) and θ(t) are two regularizations of the step function.
FIG. 3. (Top left panel) Decay of the lattice filling for an initially half-filled system and different initial temperatures.The dependence of the initial temperature is very weak (note the log-log scale).(Bottom left panel) change of kinetic energies with different initial temperatures.Systems with a initial temperature undergo a larger change in the kinetic energy.(Upper right panel) Time evolution of the distribution of Jordan-Wigner fermions for an inverse temperature kBT = J (i.e.βJ = 1).Notice that a large depletion in the distribution near k = ±π/2.(Bottom right panel) Evolution of the distribution of Jordan-Wigner fermions for kBT = J/10 (i.e.βJ = 10).