Theoretical methods to treat a single dissipative bosonic mode coupled globally to an interacting many body system

We present two approaches capable of describing the dynamics of an interacting many body system on a lattice coupled globally to a dissipative bosonic mode. Physical realizations are for example ultracold atom gases in optical lattice coupled to a photonic mode of an optical cavity or electronic gases in solids coupled to THz cavity fields. The first approach, applicable for large dissipation strengths and any system size, is a variant of the many-body adiabatic elimination method for investigating the long time dynamics of the system. The second method extends the time-dependent matrix product techniques to capture the global coupling of the interacting particles to the bosonic mode and its open nature. It gives numerically exact results for small to intermediate system sizes. As a benchmark for our methods we perform the full quantum evolution of a Bose-Hubbard chain coupled to a cavity mode. We show that important deviations from the mean field behavior occur when considering the full atoms cavity coupling [1].


I. INTRODUCTION
The coupling of quantum matter to quantum light has been achieved in numerous experimental platforms. Examples of such realizations are ultracold atomic gases coupled to optical cavities [2,3], electron gases in solids coupled to THz cavities [4][5][6], or superconducting artificial atoms coupled to on-chip cavities [7,8]. These systems open the exciting possibilities to study self-organization phenomena of light and matter [9][10][11]. Novel phenomena arise from the interplay of the long range interactions and dissipative nature induced by the cavity field and the short range interactions between the atomic degrees of freedom.
The framework of most theoretical treatments of coupled atomic cavity systems so far was based on the mean field decoupling of the cavity field and the atoms [9,18,27]. This mean field approach simplifies the numerous technical difficulties introduced by the description of the full atom-photon coupling. Within this approach, the cavity field is assumed in a coherent state and adiabatically eliminated. This results in an effective Hamiltonian for the atoms with a self-consistency equation which is typically solved for the ground state. Deviations of this mean-field treatment have been found by taking the exact coupling between the atomic and photonic states correctly into account for small systems of one or two atoms, or two sites [28][29][30][31][32][33], or in closed systems [34]. This calls for new methods which can also treat larger atomic ensembles globally coupled to bosonic FIG. 1. Sketch of a chain of interacting particles (e.g. atoms or electrons) coupled to a single bosonic quantum mode (cavity fields or phononic modes) . The bosonic mode has a dissipative nature and it is coupled to every site of the chain. The coupling strength can vary from site to site.

fields.
In this work, we develop two methods capable of capturing the exact coupling and the dissipative nature of the combined system. The first method is an extension of the many body adiabatic elimination, valid for large dissipation strengths to the combined system. Within this method any system size can be considered. It is valid for relatively long times for which system dynamics is dominated by virtual processes around the dissipation free subspace. In particular, this method provides insights about the steady state of the system. The second method consists in quasi-exact numerical simulations based on matrix product states (MPS), which can perform efficiently the full quantum time evolution of the arXiv:2004.11807v1 [cond-mat.quant-gas] 21 Apr 2020 coupled system. This method is numerically exact and can deal with small to intermediate system sizes.
Whereas these methods are very generally applicable for quantum many body systems with short range interaction coupled globally to a single dissipative bosonic mode, we benchmark the presented methods for a Bose-Hubbard chain coupled to a cavity mode and transversely pumped with a standing-wave laser beam. We concentrate here on the description of the methods and their performance. The physical effects obtained in this system, which go beyond the established mean field results, are presented in Ref. [1].
In Sec. II we describe the general setup of interacting particles on a chain coupled to a single bosonic field. Further, we describe the model for the interacting atoms coupled to an optical cavity for which the benchmarks are performed. In Sec. III we develop a variant of the many body adiabatic elimination method for the combined atom-cavity system and analyze the obtained steady state. The numerically exact tMPS method for coupled atomic cavity systems is presented in Sec. IV. We discuss in detail its implementation and convergence properties.

II. MODEL
In this work we consider dissipative systems of interacting particles globally coupled to a bosonic field, as sketched in Fig. 1. The particles can for example describe atoms or electrons and the bosonic quantum field can be for example a photonic field of a cavity or a long lived phononic mode. Generically, these systems can be described by a Lindblad equation for the density operator ρ given by [9,18,35,36] where a and a † are the annihilation and creation operators for the bosonic mode. The dissipative term proportional to the dissipation strength Γ takes into account the losses from the bosonic mode. It is described by a Lindblad form of the dissipator where the jump operator is the annihilation operator a of the bosonic mode. In a cavity these can be due to the imperfections of the mirrors and for phononic modes it describes the decay into a bath of phononic modes. Let us note, that a generalization of the developed methods to any set of jump operators acting only on either the bosonic or the particle part of the system is straightforward. The methods that we present in this work can deal with a Hamiltonian of the following form, H = H c + H chain + H ac . Where H c contains only operators of the bosonic field, H chain is an interacting short-range one dimensional Hamiltonian for the many body degrees of freedom. H ac ∝ (a + a † )A couples the bosonic field to a global particle operator, where A is a sum over operators which act on one, or at most two, atomic sites.
FIG. 2. Sketch of the bosonic atoms confined in a onedimensional chain in an optical cavity. The atoms tunnel with the amplitude J and have an on-site interaction of strength U . The coupling of the atoms to the cavity is realized with a retroreflected transverse pump beam. As the lattice spacing is commensurate with half of the wavelength of the cavity mode, the cavity field is coupled to the total imbalance between the odd and even sites of the chain. The strength of the coupling is controlled by the pump amplitude Ω. The cavity is losing photons with the dissipation strength Γ, due to the imperfections of the mirrors.
We will benchmark the developed methods using interacting bosons confined to a chain coupled to a single cavity mode transversely pumped with a standing-wave laser beam, as depicted in Fig. 2. However, the methods are easily adaptable to interacting spins or interacting fermions. In the considered model, the Hamiltonian has the form [9,18,27] The term H c describes the cavity mode in the rotating frame of the pump beam, with a detuning between the cavity mode and the transverse pump beam δ = ω c − ω p . The operators b j and b † j are the bosonic annihilation and creation operators of the atoms on site j and n j = b † j b j . L denotes the number of sites of the chain and the total number of bosonic atoms is N . For the atomic part of the Hamiltonian we have the terms H kin which describes the tunneling processes of the atoms with the amplitude J and the term H int representing the repulsive on-site interaction of strength U > 0. The term H ac gives the coupling between the cavity field and the total imbalance between the odd and even sites of the chain, ∆, with the effective pump amplitude Ω. This coupling is realized due to the assumed commensurability of the cavity mode with twice the periodicity of the lattice spacing within the chain [18].
The methods and results presented in this work are to be put in contrast with the approach of adiabatically eliminating the cavity field by a mean field decoupling of the atoms and the cavity mode [9]. In this crude approximation, after eliminating the cavity field one obtains an effective Hamiltonian for the atoms, which for the system presented in Eqs. (1)-(2), is given by The parameter V c has to be determined selfconsistently as it depends on the expectation value of odd-even imbalance, V c = 2 Ω 2 δ δ 2 +Γ 2 /4 ∆ . In this approach above a certain threshold Ω MF,c √ N the cavity field a takes a finite value and the atoms self-organize into a density modulated pattern either on the odd or even sites of the chain breaking spontaneously the Z 2 symmetry of the effective Hamiltonian. The steady state is a pure state composed of a product state between the atomic and photonic sector ρ MF = |α(∆ eff ), ∆ eff α(∆ eff ), ∆ eff |. The photonic mode is in the coherent state α(∆ eff ) with α(∆) = Ω δ−iΓ/2 ∆ and the corresponding average photon number is n MF = Ω 2 δ 2 +Γ 2 /4 ∆ 2 eff . The atomic state |∆ eff denotes the ground state of the effective Hamiltonian with the self-consistency condition. The effective imbalance ∆ eff is defined as the expectation value of the odd-even imbalance in the ground state of the effective Hamiltonian.

III. MANY BODY ADIABATIC ELIMINATION FORMALISM
In order to understand the long-time behavior of our system in the strongly dissipative regime, we employ the many body adiabatic elimination method [37][38][39][40][41]. In this section we describe the many body adiabatic elimination formalism and how to apply it to the photon mode coupled to the interacting atoms.

A. Derivation of the equation of motion
We assume that we can consider the kinetic energy term, H kin , as a perturbation ( Γ Ω, δ J) compared to the other terms in the Liouvillian L 0 = − i [H c + H int + H ac , ·] + D(·). This approach will give an insight into the effective dynamics of the density matrix in the decoherence free subspace of L 0 , i.e. the space Λ 0 formed by all density matrices ρ 0 which are eigenstates of the superoperator L 0 with vanishing real part of the eigenvalue. The other spaces, Λ m , formed by the right eigenvectors corresponding to eigenvalues with equal non-zero real part are only considered within perturbation theory.
In Fig. 3 we sketch the decoherence free subspace Λ 0 and two different subspaces Λ 1 , Λ 2 , together with the action of the Liovillian L 0 and the perturbation L kin = − i [H kin , ·] which connects the different subspaces. If we consider only contributions from the subspace, Λ 1 , that can be accessed via one hopping event, the effective dynamics for the elements of the decoherence free subspace is given by [38,41] where ρ 0 ∈ Λ 0 is an eigenstate of L 0 with vanishing real part of the eigenvalue, i.e. L 0 [ρ 0 ] = λ 0 ρ 0 with Re(λ 0 ) = 0. The operators P 0 and P 1 are the projectors onto the subspaces Λ 0 and Λ 1 , respectively. In the following we need to determine the elements of the decoherence free subspace Λ 0 and of Λ 1 . Solving the eigenvalue equation belonging to L 0 is already complex for the system we consider. However, we find that a set of right eigenstates of L 0 is given by states of the form ρ = |α(∆); n 1 , . . . , n L α(∆ ); n 1 , . . . , n L | .
At this point we do not assure that these states are physical density matrices. The atomic part is given by Fock states with the odd-even imbalances ∆ = j (−1) j n j and ∆ = j (−1) j n j and its total interaction energies u = U 2 j n j (n j − 1) and u = U 2 j n j (n j − 1). The photons are in a coherent state which depends on the atomic imbalance The corresponding eigenvalues for the right eigenvectors in Eq. (4) are given by For ∆ = ∆ the real part of the eigenvalues is zero. Thus, the states in Eq. (4) with ∆ = ∆ lie in the decoherence free subspace of L 0 . Interestingly, the eigenstates with ∆ = ∆ , but with different interaction energies u = u have purely imaginary eigenvalues. The subspace which can be accessed via a hopping event from the decoherence free subspace is given for the states in which ∆ = ∆ ± 2.
In order to identify the steady states we need to solve ∂ ∂t ρ 0 = 0. If we look at the coefficients of the diagonal terms of the decoherence free subspace, |α(∆); n 1 , . . . , n L α(∆); n 1 , . . . , n L |, in Eq. (7), we ob-serve that their sum is zero, thus the mixed state given by is a steady state of the system. Here we sum over all possible density configurations {n j } and N is the number of these configurations, which is the number of ways one can arrange N identical particles in L sites, In Ref. [26] the authors consider the same model, but they eliminate the cavity field and analyze the obtained effective Liouvillian in the atomic sector. As in their case the effective jump operators are Hermitian, it follows directly that the fully mixed state is a steady state. In contrast, in our analysis we consider the full Liouvillian Eqs. (1) and (2), including the photonic degrees of freedom, and because the jump operator (annihilation operator of the cavity mode a) is not Hermitian we need to perform the complicated many body adiabatic elimination in order to obtain insights into the nature of the steady state.

B. Properties of the steady state
With the many body adiabatic elimination method we obtain a steady state ρ mix which is very different in nature compared with the expected mean field state. The mean field state is the ground state of the effective Hamiltonian and by this a pure state with a coherent state in the photonic sector and a density wave in the atomic sector. In contrast, tracing out the photonic mode in ρ mix leads to a fully mixed atomic sector corresponding to an infinite temperature state, which is very different from a pure ground state. Moreover, the steady state, ρ mix , is a mixture of separable states, thus no entanglement is present between the photons and the atoms, but the strong cavity-atoms coupling is reflected by the fact that in each of the pure states present in the mixture the cavity field is fully determined by the atomic density profile.
This very distinct nature of the steady state is also reflected in the physical observables. Due to the fully mixed atomic sector, the density-density correlations have a flat profile. Therefore, the staggering of the density-density correlations vanishes. ρ mix has a zero expectation value of the cavity field, a = 0, capturing the weak Z 2 symmetry of the system [1], but has a finite expectation value of the photon number. The average photon number for the state ρ mix , Eq. (8), is given by where the sum is taken over the set ∆ ∈ {−N, −N + 2, . . . , N − 2, N } and c ∆ being he number of states with The behavior is consistent with a N −1 scaling of a † a /N . In the inset we have the dependence of the photon number, a † a , on the particle number, N , which seems to saturate at large N . a certain imbalance ∆, given by By plotting the scaled photon number a † a /N at a fixed filling N/L, Fig. 4, we can see that this quantity vanishes as N −1 at large N . This implies that even though the scaled photon density per atom is finite for any finite size system, it goes to zero in the thermodynamic limit, N → ∞. Thus, the many body adiabatic elimination method tells us that in the thermodynamic limit at large dissipation strengths the system is no longer in a superradiant state with a finite a † a /N , but in a state with an average number of zero photons and a fully mixed atomic sector. This is very distinct to the normal state predicted by the mean-field approach. Despite the fact that in both states the average photon number vanishes, the atomic part of the mean-field state is a pure state and not the infinite temperature state as predicted by adiabatic elimination. Therefore, our results also question the nature of the transition predicted by the mean-field approach between the superradiant state and the normal state.

IV. TIME-DEPENDENT MATRIX PRODUCT STATE (tMPS) METHOD FOR COMBINED ATOM-CAVITY SYSTEMS
In this section we describe the numerical exact method based on matrix product states (MPS) we developed to perform the quantum time evolution of the combined system.
A. Details of the tMPS method for the coupled photon-atom system The considered dissipative system of atoms coupled to an optical cavity poses several challenges for its treatment via MPS based methods. The first difficulty arises due to the in principle arbitrarily large dimension of the Hilbert space of the cavity field. The second obstacle is the global coupling of the cavity mode to the interacting atoms. The third challenge is the dissipative nature of the combined system due to the photon losses. In the following we describe how our implementation overcomes all these difficulties. We implement the newly developed algorithm efficiently using the ITensor library [42].
We begin by presenting how the dissipative aspect of the considered models is included in the numerical method. For the simulation of the dissipative many body quantum systems the time-evolution of the density matrix following the Lindblad equation needs to be determined. State of the art are two different routes: the first is the purification approach which relies on the rewriting of the density matrix with a larger dimension [43,44]. The second is the stochastic unravelling of the master equation using quantum trajectories [45,46]. This approach has the advantage of simulating the timeevolution of wavefunctions instead of density matrices at the disadvantage of a stochastic sampling. We have chosen as a first implementation the stochastic unraveling of the master equation. This has in particular two reasons: First, already the representation of the interacting ground state of the bosonic atoms as initial state would have been demanding in the purification approach. Secondly, the additional presence of the large Hilbert space of the photons would have increased the required matrix dimension.
In the stochastic formulation we take good quantum numbers into account for the atomic sector. The efficient combination of the stochastic unraveling of the master equation with the matrix product state methods has been relatively recent and only few groups have efficient implementations taking conserved quantum numbers into account (see e.g. [47][48][49][50][51][52][53]).
In order to apply the stochastic unravelling procedure, the time-evolution of many trajectories of pure states is sampled and finally the results are averaged. The initial states for the trajectories are drawn corresponding to their probability weights in the initial density matrix. Then a stochastic time-evolution is performed for each trajectory which is described in the following: • A random number η is drawn from the interval [0, 1).
• For each trajectory, the time evolution is performed for a time step with a non-unitary time evolution operator, corresponding to the effective Hamiltonian,H = H − i 2 Γa † a. • Since the effective Hamiltonian is not Hermitian, this leads to a decay of the norm of state in time.
FIG. 5. The graphical representation of one time step based on the Trotter-Suzuki decomposition described in the text, Eq. 11. The first (red) site in the graphical representation of the MPS structure corresponds to the cavity mode and the rest to the atomic sites. To be noted that the cavity mode index marked with a red line has a large local dimension. Green boxes represent the application of the two site gates of the atomic terms of the time-evolution after Trotter-Suzuki decomposition followed by an SVD compression step. With orange we depict the large tensor corresponding to the time evolution of the cavity and cavity-atoms coupling terms of the Hamiltonian. Its application is detailed in Fig. 7.
The non-unitary deterministic time evolution is performed until the norm is smaller than a threshold posed by the random number η.
• A quantum jump is performed by applying the jump operator a onto the wavefunction and the state is normalized.
• The described procedure is repeated until the required final time is reached.
One can show [45,46] that taking the Monte Carlo average over all sampled quantum trajectories, the described time evolution reproduces the Lindblad dynamics correctly up to the first order in the chosen time step. In order to achieve convergence in the computed quantities many trajectories are required and the time steps need to be chosen small enough to avoid multiple jumps within one time step. In our case, because the jump operator only acts on the photonic space, we need to sample several hundred trajectories, as discussed in the Sec. IV B.
In order to perform the time evolution within the MPS formalism, we represent the wave function as a matrix product state (MPS) [54], with the first site initially corresponding to the cavity mode and the rest to the atomic lattice using a Fock basis for each site (see Fig. 5). In order to take care of the in principle arbitrarily large Hilbert space of the photonic mode, we introduce a cutoff for the dimension of the local Hilbert space of the photonic site, which is dynamically adapted during the time evolution. This is done by setting a truncation goal of the photonic distribution and the details are given in Sec. IV E. For benchmarking we also present results in which a fixed dimension of the photonic Hilbert space is used.
The global range coupling between the cavity mode and all the atomic sites makes the use of the tMPS implementation for short-range Hamiltonians based on the Trotter-Suzuki decomposition impossible. Thus, in order to take both the global coupling between photons and atoms and the short range interaction of the atoms into account, we develop a variant of the tMPS based on the dynamical deformation of the MPS structure. The dynamical deformation allows one to alter the order of the sites as needed using swap gates [50,54,55]. Previous variants of the MPS time evolution with swap gates dealt with short-range interaction in two dimensional models [55], or spin-boson models [50,56]. Our implementation, in contrast, can efficiently deal with interacting bosonic models globally coupled to the dissipative photonic field. An adaptation to fermionic and spin systems coupled to photonic modes is straightforward.
In the following, we describe our procedure for performing a time step dt with the effective Hamiltonian, H. It is based on the Trotter-Suzuki decomposition of the time evolution propagator in combination with swap gates. The terms are split in order to separate the terms containing the cavity field operators and the remaining terms This decomposition is valid to the order O(dt 3 ) in the time-step. The evolution given by the operator e − idt 2 (H kin +Hint) which only contains the atomic operators is computed as in the standard tMPS algorithm for short-range interactions [57,58] by a further decomposition into two site gates. The two site gates are applied to the MPS followed by a compression step via a singular value decomposition (SVD) in the order sketched in Fig. 5. We mention that the boundary gates at the beginning and the end of the bosonic chain differ from the gates applied in the bulk.
For the operator e − idt (Hac+Hc− i 2 Γa † a) which contains the global coupling to the cavity field, we use the fact that we can decompose H ac such that each term only acts on two sites -even though distant ones- This means that we need to apply two-site operators where the two sites are not neighbors in the initial MPS representation. In order to solve this problem, we adapt the structure of the MPS while applying the timeevolution gates such that we bring the two sites on which the operator acts next to each other. This approach is implemented using swap gates, where the action of the swap gates consists in the swapping of the physical indices of two neighboring MPS matrices, i.e. Here S s,σi is the swap operator and M σ1 ...M s M σi ...M σ L is the weight of the state |σ 1 , ..., s, σ i , ..., σ L in the MPS form with s the index of the cavity mode site and σ i the index for the bosonic atoms. In Fig. 6 we sketch how the swap gate acts on two MPS sites and changes their order. The swap gates are constructed from two Kronecker delta functions, each between indices of the same nature, but different sites, i.e. in Fig. 6(a) we have a Kronecker delta from the cavity index s at the first site to the cavity index s at the second site (red curve) and a Kronecker delta from the atomic index σ at the second site to the atomic index σ at the first site. The next step is the application of the swap gate onto the MPS wavefunction and obtaining a two-site tensor with swapped indices [ Fig. 6(b)]. Finally an SVD decomposition is performed to restore the MPS structure. Thus, using the swap gates we can apply the operator e − idt (Hac+Hc− i 2 Γa † a) onto the wavefunction as a series of two-site gates, as depicted in Fig. 7. No additional error is introduced by the swap gates, except the SVD truncation error.
The implemented time-evolution method has an error of the order O(Ldt 2 ) at a certain final time t, stemming from the Trotter-Suzuki decomposition. However, the stochastic unravelling is only valid up to first order to dt, such that we expect that this limits the choice of the time step. A detailed analysis of the contributing errors and improvements of the method using different timeevolution schemes will be performed in future.
To improve the performance of our tMPS algorithm we take good quantum numbers into account, by noting that our system preserves the total number of bosonic atoms.

B. Numerical convergence
The convergence of the numerical method is controlled by several parameters. The stochastic unravelling of the master equation with quantum trajectories requires an averaging of a sufficiently large number of trajectories. Additionally, the time step dt must be chosen small enough in order to avoid the occurence of multiple jumps in one time step. The next source of error comes from the Trotter-Suzuki decomposition of the time evolution operator. Again this require that the time-step dt is small enough. Finally, we introduce an additional error by representing our wavefunctions as a matrix product state with finite local and bond dimensions. This implies a cut-off, N pho , of the local Hilbert space of the photons and in some situations also for the bosonic atoms. The procedure to dynamically adjust the cut-off, N pho , will be presented in Sec. IV E. Additionally, the introduction of the finite bond dimension using SVD leads to a so-called truncation error. To control the bond dimension of the MPS we impose a truncation error goal , thus, in each compression step, after the application of a time evolution or swap gate onto the MPS, the number of states kept is such that the truncation error is smaller than . Let us note that as in the case of the time-dependent MPS [57,58], the arising errors are not independent and therefore, a careful analysis needs to be performed.
In the following we discuss typical values of the convergence parameters and their influence on the results. The discussion of the truncation error and the von Neumann entanglement of the trajectories is presented in Sec. IV C.
Stochastic error. We estimate the error of having a finite number of quantum trajectories included in the Monte Carlo average by computing the standard deviation of the mean for the measured expectation value of an operator E σ (E(t)) = 1 where R is the total number of samples, |ψ r (t) the time evolved wavefunction of the trajectory labelled by r, and E the statistical average over all quantum trajectories. For the numerical data presented, we show this error with the curves when averages are presented. Typically, we average over at least 500 trajectories, which ensures that for the physical parameters considered the relative error in the expectation value of the photon number is smaller than 3% (see for example Fig. 8 upper panel). For the cases when the photon number is small, a † a 1, either at small coupling Ω or large dissipation strengths Γ, we average over 750 trajectories to obtain the same relative error, as the fluctuations have a greater influence (Fig. 8  lower panel).
Cut-off on the dimension of the local Hilbert spaces. The dimension of the local Hilbert space for the photon can be infinite and thus, a cut-off for its dimension is needed in the numerical implementation. In the following we will denote the cut-off N pho , referring to the maximal number of photons that we can capture and noting that we use N pho + 1 Fock states as we also have the vacuum state. In order to identify more clearly the influence of the cut-off on the results we used a fixed cut-off for the photonic site in this section. However, in Sec. IV E we present a more efficient approach by implementing an adaptive photonic local dimension, since the required cut-off varies considerably in time and with the trajectories. Examples with a fixed cut-off with all other parameters fixed are shown in Fig. 8. For a given set of parameters we observe that above a certain value of the cut-off N pho the average value of the photon number is only slightly varying with increasing the cut-off. In particular for the presented situation its variation for N pho ≥ 35 becomes lower that the error bars of the Monte Carlo averaging. However, choosing a too low cutoff e.g. N pho ≥ 25 (N pho ≥ 4) in Fig. 8, leads to misleading results for both the time-evolution and the reached long time values. Note, that even though the long time value lies around 12 (below 1) photons, taking double this value for the cutoff i.e. 25 (4) is not sufficient. The required cut-off depends very much on the physical parameters. Therefore, one needs to consider each case separately, which can result in very different values for the cut-off.
Since we consider bosonic atoms, the maximal possible local dimension for the atomic sites equals the total atom number plus one for the possibility to have an empty site. We found that often this very large local dimension is needed, which also strongly restricts the total number of atoms which can be efficiently considered. However, in some situations a reduced dimension of the local bosonic site can be taken. For example, we checked that for the parameter sets with Γ/J = 1 and not too large values of the dissipation strength a maximal local dimension of five instead of six is sufficient.
Influence of the time step. The dependence of the results on the value of the time step is more involved. This is due to the fact that the time step controls both the convergence of the stochastic sampling process and the Trotter-Suzuki decomposition. Further, as in the normal time-dependent MPS, the time step interplays with the truncation error in a non-trivial fashion, since a smaller time-step requires more truncations and therefore results in an increased truncation error [58]. Therefore, the values used needs to be adjusted very carefully depending not only on the physical parameters but also on the convergence parameters of the model.
In Fig. 9 we show an example of the variation obtained fixing all parameters beside the time step dtJ. A relatively rapid convergence is seen using time-steps between dtJ/ = 0.01 − 0.05. In particular, the convergence is in agreement with the expected linear behaviour in the time-step dt which suggests a well justified extrapolation method. For small values of Γ, here Γ/J = 1, the error induced by the time-step remains larger than the error of the statistical error. The extrapolated value lies a bit above the shown results at a finite time step. In contrast for the case of large Γ, the statistical error is dominating the results and the extrapolated result lies within the statistical error bars of the smallest time steps.

C. Entanglement of quantum trajectories
One of the most important convergence parameters is the bond dimension used within the SVD compressions. One measure of this is the truncation error , the sum of the neglected eigenvalues of the reduced density matrix in the SVD compression. Another measure of the decay of the eigenvalues is the von Neumann entanglement entropy, S vN . We analyze in the following the behavior of these two quantities for different parameters.
We first look at the dependence on the truncation error of the singular value decomposition performed in the time evolution gates and swap gates shown in Fig. 10. As the truncation error is chosen relatively small, the results only weakly depend on the value of the maximal truncation error . In particular, the deviations are of the order to the statistical error. Thus, we are confident that a truncation error of ∼ 10 −11 − 10 −12 provides an accurate description of the considered states in the matrix product form. One can also monitor the maximal bond dimension needed in the MPS representation in order to achieve the set truncation error goal, as depicted in Fig. 11. The maximal bond dimension increases considerably with lowering the truncation error. However, for the smallest chosen truncation errors the maximal bond dimension saturates. We can observe that the bond dimension needed to describe the system is between 100 and 300 even for a system of size L = 10.
In the following we turn to the von Neumann entropy of the quantum trajectories to monitor the coupling of the photonic and atomic sectors and the correlations within the atomic chain. We note that the entanglement entropy of the quantum trajectories is not a direct measure of the entanglement present in the density matrix resulting from the Monte Carlo averaging process. However, S vN provides valuable information about how well our MPS method captures the entanglement present in the trajectories.
In the following, we consider two different bipartitions of the MPS, as depicted in Fig. 12(a), between the cavity site and the rest of the atomic chain, bond l = 1, and in the middle of the atomic chain, where one half also contains the cavity site, bond l = L/2 + 1. This is motivated by our finding that the maximum of S vN throughout the atomic chain occurs at the bond l = L/2 + 1.
In Figs. 12(b)-(c) we present the time evolution of the entanglement entropy for the Monte Carlo average, the maximum entanglement entropy of the sampled quantum trajectories and for a few single trajectories, for the two considered bipartitions. We observe that for both bipartitions S vN saturates to a finite value in time, both for the average and maximum values. In Fig. 12(d) it can be seen that the values of S vN only change within the Monte Carlo averaging uncertainty for all considered truncation errors. Thus, we can be confident that our method captures the dynamics of our system correctly up to long times.
In the long time limit, the von Neumann entropy takes finite values for both bipartitions and both parameter sets considered. In Fig. 12(b) we see that at low dissipation strength the average entanglement entropy computed between the photon mode and the atoms (l = 1) becomes close to log(2) at long times. This signals a coherent superposition of two states which we attribute to a superposition of states with a different sign of the photon field [1]. The value of the entanglement within the chain is larger which points to the contribution of several states in the superposition. In Fig. 13 we computed the entanglement entropy for different system sizes, for two different parameter sets. We observe that in both cases the entanglement present in the quantum trajectories between the photon mode and the atoms seems stable with the respect to the system size. This further supports the claim that a coherent superposition of two system size independent states contribute, as it is the case for the states with the different sign of the photon field [1]. For the bond l = L/2 + 1 we see that the value at which the entanglement entropy saturates increases with the system size, indicating that the system might be in a gapless phase or that the system size is not yet long enough to cover the correlation length of the gapped state.

D. Finite size effects
As we have seen in the previous subsections, at long times the considered quantities have become almost constant in time. Typically, such a regime is reached long before tJ/ ≈ 50, as shown in the time evolution plots, for example Fig. 8 for the photon number, or Fig. 12 for the von Neumann entropy. Therefore, in this subsection we compare the values at late times for different system sizes, as we interpret these as very good approximations of the steady state values.
In order to evaluate the finite size effects we analyze how the transition from the normal state to the selforganized state takes place for different system sizes. In Fig. 14(a) we scale the photon number and the atomscavity coupling with the number of particles. For a comparison we show both the mean field and the numerically exact tMPS method results. Both show only small deviations with increasing the system size. In particular, in the mean field results the transition to the self-organizes phase starts later and becomes steeper with increasing system size. In the tMPS results, the rise of the photon number also seems to occur for a bit larger scaled pump strength. However, the effect are very small which leads us to the expectation that our main findings will remain valid for large systems.
To further support this, in Fig. 14(b) the scaled photon number is plotted as a function of the dissipation strength for large dissipation strength. The numerical results are compared with the many body adiabatic elimination results, as the state ρ mix can be evaluated for any system size. We observe that the agreement is very good at large dissipation strengths for all values of L considered. The scaled photon number is slightly decreasing for larger systems sizes in both approaches. This is consistent with the expected vanishing of the scaled photon number in the thermodynamic limit behavior for ρ mix , as We have shown in Sec. IV B that the cut-off of the local Hilbert space for the photons is an important convergence parameter of the numerical simulations. In particular, a sufficiently large N pho -typically about triple the average value-is required to capture the dynamics correctly (Fig. 8). Often during the time evolution the photon number varies considerably, as for example for the case in Fig. 8(a). In this case at short times, tJ/ 5, there is a sudden increase of the number of photons in the cavity, larger than the value at late times. This suggests that we need a much larger N pho to accurately describe the photonic state at short times, than we would need at later times. Therefore, we decided to optimize our implementation by adapting the local dimension for the cavity site during the time evolution. This improvement is similar in spirit with the recent developments regarding the time-dependent MPS methods with local basis optimization. Refs. [59,60] apply the local basis optimization idea [61] to the time evolution of the Holstein model of fermions locally coupled to phononic modes. In this approach one rotates the local Hilbert space adaptively into an optimized basis that can be truncated. In our case we find that even without changing the Fock basis for the photons we can dynamically adapt the number of states considered. The investigation whether an optimization of the photonic basis is further improving the algorithm is left for future implementations.
In order to implement such an adaptive cut-off, we monitor the evolution of the photon number distribution by measuring the occupation P m of the photonic Fock states with photon numbers m close to the cut-off value in each time-step. We adapt the cutoff using thresholds for the photonic state occupation. To be more precise, the procedure is as following: At a certain time in the evolution of a single quantum trajectory we have a cut-off N pho (t). Depending whether the photon number should increase or decrease we encounter two different cases: (i) The occupation P N pho of the photonic Fock states with the largest photon number is smaller than a chosen threshold p d . This signals that the cutoff can be decreased. In order to do this, we find the photonic Fock state m * ≤ N pho with the largest photon number whose occupation is above the threshold, i.e. P m * ≥ p d . We change the cut off of the local dimension of the cavity site of the MPS such that the maximal photon number is N pho (t + dt) = m * + 1 the for the time step t + dt.
(ii) The occupation P N pho of the photonic Fock states with the maximum photon number is larger than a second chosen threshold p i ≥ p d , i.e. P N pho ≥ p i . This signals that the photon number should increase. We increase the local dimension of the cavity site of the MPS at the next time step N pho (t + dt) = N pho (t) + 2.
The numerical parameters that control the convergence of the method are now the two threshold values, p d and p i .
In Fig. 15 and Fig. 16 we check our procedure of adapting the local dimension for the cavity site by comparing with the results for a fixed converged cut-off. In Fig. 15(a) and Fig. 16(a) we represent the Monte Carlo average of the photon number for two different sets of sampled trajectories with a fixed cut-off and the Monte Carlo average of the photon number with an adapted cutoff for different p i and p d . We can observe that, except for the case with p i = 5×10 −2 and p d = 10 −3 , the results for an adapted cut-off agree, within the Monte Carlo averaging error, with the ones for a fixed cut-off. The evo-lution in time of the cut-off can be seen in Fig. 15(b) and Fig. 16(b). We note that at short times, tJ/ < 4, we always keep the photonic cut-off fixed and relatively large in order to capture the sudden increase in the number of photons in the cavity. But at later times we can observe that in all considered cases the cut-off is approximately 25% smaller than the initial value. This implies a significant speed up of the tMPS method, as the large local dimension of the cavity site being one of the bottlenecks of the method. As a rough estimate, for the parameters used in Fig. 15 the runtime was with 50% smaller compared with the case with a fixed cut-off and for the parameters used in Fig. 16 with 25% smaller. In Fig. 15(c) and Fig. 16(c) we plot the occupation of the photon number states at the final time, tJ/ = 49.75, to check the agreement of the entire photon number distribution. A very good agreement is found except for the case with p i = 5 × 10 −2 and p d = 10 −3 . We also plot in Fig. 15(c) and Fig. 16(c) the photon number distribution at tJ/ = 1.75, close to the peak in the photon number, to show that at short times many photon number states are occupied. We note that we also verified the accuracy of this method at the level of single quantum trajectories, not only by analyzing the Monte Carlo average. The improvement brought by this new development is dependent on the physical parameters of the model, but, roughly, has a more important impact when the average photon number is larger, where the difference between the maximum photon number at short times and the steady state value is larger. For example for the parameters in Fig. 15 we manage to lower the local dimension of the cavity site with more than 10 states compared to the fixed cut-off previously used, but for the parameters in Fig. 16 we have lowered the local dimension with 5 states, due to the lower photon number in the cavity.

V. CONCLUSIONS
To summarize, we described in detail two methods capable of tackling both short and global range interactions of an interacting many body system coupled to a single dissipative bosonic mode. We benchmark the methods with the example of a Bose-Hubbard chain coupled to an optical cavity.
The first method is based on the many body adiabatic elimination formalism and applicable for strong dissipation strength. We show how to derive the steady state of the system in the limit of large dissipation strength within the many body adiabatic elimination formalism. The resulting state is a highly mixed state and the reduced density matrix in the atomic sector corresponds to an infinite temperature state. Our results can be evaluated for any system size and we analyzed its physical properties and the thermodynamic limit. In particular, we showed that the obtained steady state has a very different nature compared to the expected mean field state.
As a second algorithm we developed a quasi-exact tMPS method in order to determine the full quantum evolution towards the steady state of the interacting bosonic chain coupled to the cavity. This implementation deals with all the challenges posed by the atomcavity system: We employ the stochastic unravelling of the master equation to simulate the Lindblad equation, Eqs. (1)- (2). The global coupling of the cavity to the atoms is tackled via the dynamical deformation of the MPS structure with swap gates. The efficient simulation of the very large photonic Hilbert space is ensured by its dynamically adapted cut-off. We analyze carefully the convergence of the method for different parameter sets. In particular, we monitored the time dependence of the von Neumann entropy of the quantum trajectory in order to ensure that we properly capture the entanglement between the cavity and the atoms and within the atomic chain.
Both methods open the possibility to treat many body systems coupled to a dissipative bosonic mode beyond the often applied mean field methods. The presented algorithm is easily adapted to fermionic or spin many body systems. Therefore, the methods will have a wide range of application and we expect that many new physical findings will rely on these methods.