Dissipation induced Tonks-Girardeau gas in an optical lattice

We present a theoretical investigation of a lattice Tonks-gas that is created by inelastic, instead of elastic interactions. An analytical calculation shows that in the limit of strong two-body losses, the dynamics of the system is effectively that of a Tonks gas. We also derive an analytic expression for the effective loss rate. We find good agreement between these analytical results and results from a rigorous numerical calculation. The Tonks character of the gas is visible both in a reduced effective loss rate and in the momentum distribution of the gas.


I. INTRODUCTION
A Tonks-Girardeau gas is a one dimensional (1D) gas of bosons where a repulsive interaction dominates all other energy scales. It was predicted long ago [1,2] that repulsion can mimic the Pauli exclusion principle and that, for one dimension, the wavefunction for bosons and fermions would be the same, except for a simple transformation. Not only that, but the excitation spectrum of these onedimensional hard core bosons can also be mapped to that of a fermionic system. A few years ago, experiments observed a Tonks-Girardeau gas using ultracold gases in optical lattices with strong repulsive interactions [3,4,5].
In a recent experiment [6] we showed that strong dissipation in the form of two-body losses can also simulate a Pauli exclusion principle, fermionizing a system and transforming it into a dissipative but long-lived strongly correlated gas. This equivalence was demonstrated for molecules of 87 Rb loaded in an optical lattice. For deep enough lattices, these particles exhibit strong two-body decay rates which, in appropriate units, exceed the average kinetic energy of the particles. The particles then avoid coming close together and behave like impenetrable bosons. This happens both in the continuum case of 1D tubes and when the tubes are modulated by a deep optical lattice -two situations which resemble the experiments with elastically interacting bosonic atoms in Refs. [4,5] and Ref. [3], respectively.
At least in the lattice case, the equivalence between strong dissipation and a Pauli exclusion principle can be understood in terms of the Zeno effect. Following the discussion in Ref. [6], Fig. 1 depicts a stable configuration with two molecules on neighboring sites. This configuration is connected via hopping to states with double occupancy. These states decay at a rate Γ which is much larger than the hopping amplitude J/h of the particles. Treating this as a typical three-level system from Quantum Optics, one concludes that particles stay on their original sites with only a minor loss rate of O(J 2 /h 2 Γ).
In this work we present a rigorous theoretical analysis of this system. The outline is as follows. We begin in Section II by introducing a master equation that models two-body losses for a single species of bosonic particles -for instance the molecules in Ref. [6]-. We will partic-a) 1: Two molecules sit on neighboring lattice sites. The molecules can tunnel with a hopping amplitude J/h which is much weaker than the decay rate of the resulting doubly occupied site into the vacuum. The effective model is that of impenetrable particles decaying with a much weaker O(J 2 /h 2 Γ) decay rate.
ularize the model to the case in which particles are confined by an optical lattice and explain how dissipation becomes the dominant term. In Section III we will show that in the limit of strong losses the master equation can be replaced by an effective model in which the rapid twobody decay has been eliminated. The dominant terms of the effective model are identical to the Hamiltonian of an elastic Tonks-Girardeau gas. The residual effect of losses is a slow perturbation that can be obtained analytically. In Section IV we compare both results to exact numerical simulations of the full master equation. We use Matrix Product Density Operators [7] (MPDO) to verify that both the density and momentum distributions of the particles closely resemble those of a Tonks gas. A second signature of the hard-core boson dynamics is a slowdown of the inelastic losses, for which we find good agreement between numerical and analytical estimates. In the last part of our work we offer further details about the methods and derivation used to obtain the effective Tonks gas models and the slowdown of losses. Finally, we summarize this paper in Section VI.

A. Master equation
We model the dynamics of the particles using a Markovian master equation where H is the Hamiltonian describing the unitary part of the evolution and D is a dissipator associated with the losses due to inelastic collisions. In the absence of a lattice potential [8] where H s (x) = −h 2 ∇ 2 /2m + V trap (x) is the singleparticle Hamiltonian and Ψ(x) is the bosonic field operator. The strength of the interparticle interactions is g 3D = 4πh 2 a/m, where m is the mass of a particle and a is the scattering length. Re(a) describes elastic collisions, whereas Im(a) ≤ 0 describes inelastic collisions that lead to losses.
To assert the consistency of this model, let us estimate the decay rate of the particles. When taking expectation values over the density operator, n(x) = Ψ † (x)Ψ(x), the previous master equation becomes For a condensate of molecules this gets the expected form wheren(x) is the density of bosons at a given point and the decay rate is proportional to the square of the density and to K 3D = −2 Im(g 3D )/h.

B. Optical lattice
As explained in the introduction, in this paper we want to model the experiment with molecules in a deep optical lattice [6]. For low temperatures and tight confinements, we can expect that the particles will accommodate to the motional ground state of each lattice site. Note that if this is true for individual atoms, it is even more so for molecules, because having twice the mass and twice the polarizability, they are better trapped by the same optical potential. Under these conditions, it was shown in Ref. [9] that it is convenient to expand the bosonic field operator Ψ † using Fock operators, a † k , that create particles on the k-th lattice site: The w(x) are Wannier wavefunctions associated to the states localized on each site. In particular, for the arrays of one-dimensional lattices used in Ref. [6], we can separate the wavefunction of these bosonic modes into a product of three wavepackets, w(x) = w (x)w ⊥ (y)w ⊥ (z), a less confined longitudinal one, w , and two transverse ones, w ⊥ , which are very tight due to the perfect decoupling between adjacent tubes. The advantage of the Wannier functions is that we can now perform a tight-binding approximation: anywhere outside the single-particle terms in the Hamiltonian, the overlap of different Wannier functions is neglected. This procedure transforms the unitary part of the master equation into a Bose-Hubbard model [9], and discretizes the dissipative terms as well Dρ =h Γ 4 Regarding the notation, the sum k, l extends over nearest neighbors along the same tube, |k − l| = 1. The tunneling amplitude between neighboring lattice sites is denoted by J. Finally, the on-site interaction matrix element contains both a real and imaginary part which contribute to the unitary and the dissipative parts of the master equation, respectively. The imaginary part of the interaction constant governs the decay of the number of particles per site, n k = a † k a k d dt In the cases that we will study this rate will be larger than the speed at which particles tunnel to neighboring sites, Γ J/h. To facilitate the calculations, we will group the terms in the master equation according to their strength. We introduce a superoperator V that contains the tunneling part (H J ∝ J) and a superoperator L int that describes the elastic (H el ∝ U r ) and inelastic inter- The term "superoperator" refers to the fact that D, V, and L int are linear operators acting on density matrices, not on pure states. Finally, in the following it will be important to realize that V is of order J, whereas L int is of order |U | |J|, and dominates the evolution. The following sections exploit this difference of scales to create new and simpler effective models that describe the dynamics of the molecules.

A. Second order effective theory
Our goal is to develop an effective master equation that is equivalent to (10) in the limit of strong dissipation,hΓ J, in which L int dominates. We will sketch the main ideas and refer the reader to Sect. V for more details. The process begins by identifying the eigenvalues, λ i , and eigenspaces of the dominant term, L int = i λ i P i . We then decompose the density matrix into a sum of contributions from these eigenspaces Finally, we will determine some approximate evolution equations for the different ρ i in the presence of a nonzero hopping term, V.
We only need to study three eigenspaces, corresponding to the dissipation-less states, λ 0 = 0, and to the states most immediately connected to them by the hopping. The most relevant term of our density matrix, ρ 0 ∼ O(1), is given by states with zero or one particle per site. These states do not decay in the absence of hopping and have the property that they belong to a space with zero or one particle per site. We thus have with a projector expressed in the Fock basis of occupation numbers where L is the number of sites in the lattice. The next states that we need to consider have a pair of particles on some site. As we will see later, we have two sets of states depending whether the double occupation is on the left or on the right side of the matrix coherences where we have introduced In the absence of hopping these states decay with eigenvalues λ 1a = −Γ/2 − iU r /h = −iU/h and λ 1b = −Γ/2 + iU r /h = iU * /h, respectively By neglecting higher order contributions to the density matrix, it is possible to integrate formally the projected master equation and obtain an effective model for the dominant term, ρ 0 (t). After some manipulations one arrives to the following model [Sect. V C] We will now write down explicitly and discuss the meaning and implications of these different terms.

B. Tonks-Girardeau Gas
The most important term in our effective model (16) is given by the Liouville operator L 1 = P 0 VP 0 . As shown in Sect. V D this superoperator is equivalent to a Hamiltonian of hard-core bosons, also known as a Tonks-Girardeau gas. Therefore, to lowest order with bosonic hard-core operators c † k and c k that have implicitely the restriction of one particle per site This is the main result of our paper. Namely, that a strong dissipation such as the two-body losses from our system can lead to coherent evolution. As it was mentioned in the introduction, the same result can be obtained in an alternative way. By establishing an analogy between losses and a continuous measurement, it is intuitively clear that the strong dissipation causes a Zeno effect which suppresses the process of two particles coming together and being lost. In this regime, the dynamics of the molecules must be given by a Hubbard model where doubly occupied states have been projected out.
In practice, the hard-core bosons model implies a very simple dynamics that can be tested numerically, as we do in Section IV. However, verifying the same thing with real particles in an optical lattice represents a challenging experiment. Nevertheless, one can easily check two phenomena: first, that in the regime of strong decay, hΓ J, no significant losses take place in a time scale of order 1/Γ; and second, that the effective loss rate can be estimated using the second part of our effective model (16c). This is the goal of the following sections.

C. Second order corrections
With a lengthy calculation we can rewrite the second order terms in Eq. (16) as an effective Liouvillian L 2 with both a Hamiltonian and dissipation It is convenient to introduce an operator that destroys a pair of particles in neighboring sites With these pairs the Hamiltonian part can be written with an effective strength (1/λ 1a =hi/U ) Note that this Hamiltonian contains both effective nearest neighbor interactions and a three-site hopping of the form c † k+1 n k c k−1 . Both terms are typical of the Bose-Hubbard model in the limit of strong repulsive interaction U r |J| [10]. Moreover, these Hamiltonian corrections disappear when the lossy particles do not interact elastically on-site, U r = 0.
The dissipative term is equally simple, and has a loss coefficient In the limit of weak elastic interaction between molecules we may write which shows that the decay rate is indeed O(J 2 /h 2 Γ) as anticipated.

D. Effective losses
Let us write the evolution of the total number of particles,N , under the effective master equation (16) This expression in general cannot be simplified any further, at least not without some assumption about the state with which we compute the expectation value.
The experiments described in Ref. [6] start from a state with exactly one molecule at each site and evolve only until approximately half of the particles are lost. Here further approximations are possible: first we treat the system as homogeneous and second we assume that the populations of different sites are uncorrelated [16]. We thus obtain wheren is the density and z = 2 is the coordination number of our lattice. This approximation leads to a rate equation characteristic of two-body processes A similar equation was derived previously in Ref. [6] using a more restricted theory than our effective master equation (16).

IV. NUMERICAL RESULTS
In order to study the quality of the approximations used in the effective model, we have performed numerical simulations of the full master equation (10) using MPDO [7]. This is a method that approximates the density matrix ρ(t) using a matrix product state structure. As described in Ref. [7], this variational ansatz is well suited to simulating evolution of the state under a Liouvillian like V + L int , which can be decomposed into a sum of local or nearest-neighbor terms. In our simulations we have worked with up to 48 sites and open boundary conditions, setting U r = 0 so that the different effects cannot be attributed to a repulsive interaction. We have experimented with different cutoffs, from two to four particles per site, verifying that they give similar results.
A. Comparison with the Lossless Tonks Gas Fig. 2 shows numerical results forhΓ/J = 4000. The system was initially prepared such that the N = 24 central sites of the lattice contain one particle per lattice site, while all other lattice sites are empty, similar to the state experimentally prepared in Ref. [11]. We then let the system evolve for a time t = 2h/J and measure both the density and the momentum distributions, which are plotted in Figs. 2(a) and 2(b), respectively. The solid lines were obtained with the full master equation (10) and differ only slightly from the dashed lines obtained with the lossless Tonks model of Eq. (17). The observed difference in the distributions is largely due to the reduced particle number.
In order to quantify how small the remaining difference of the position distributions is, we define where N (t) is the total number of particles at a given time, n x is the number of particles for the site x and n Tonks x is the same for the lossless Tonks gas. A similar measure can be defined in momentum space. Fig. 3(a) shows these quantities. For J/Γ → 0 the distributions converge to those of a lossless Tonks gas, as expected from the effective model.
Another important quantity in this comparison is the fraction of the state in the hard-core bosons subspace of states with at most single occupation This quantity is shown in Fig. 3(b). Comparison with part (a) shows that 1 − F tends to be much smaller than . This is because loss causes the system to evolve into an incoherent mixture of systems with different total particle number. Each of these systems remains to a pretty good approximation in the Tonks gas subspace, but the position and momentum distributions look different so that the mere scaling with total particle number yields larger deviations in .

B. Effective Loss Rate
Our analytic model suggests that the loss rate of particle number in a homogeneous system can be approximated by a two-body decay Eq. (28), which corresponds to an evolution of the density given bȳ This result applies to the case in which different sites are not very much correlated and the losses Γ dominate over the hopping, J.
In order to test this prediction, we have performed MPDO simulation of Eq. (10) using a uniformly filled lattice with L = 40 sites and N = 40 particles. We studied the evolution for a wide range of values of J/hΓ, with some examples shown in Fig. 4(a). After an initial transient that vanishes on a timescale ≈ 1/Γ, a two-body decay Eq. (28) fits the data fairly well. However, while for J/hΓ 1 the effect of the transient is negligible, for larger J/hΓ there are better ways to fit the resulting curves. One possibility, used in Ref. [6] to fit the numerical data, is to include a free parameter t 0 describing an offset on the time axis However, we found that a more accurate model is a modified decay equation with an exponentially modulated decay coefficient Here the exponential term with λ ∼ Γ represents an heuristic approximation to the transients that we have neglected in developing the effective model [See Sect. V C]. The solution of this differential equation still has two fit parameters Fig. 4(b) shows the loss rate coefficient κ from these two fits and Fig. 4(a) shows the quality of the fits for small values of the hopping. For larger values, though, the fitting becomes numerically unstable. In this regime a fit (32) behaves better. What we do not show in the previous plots is that the long term behavior of the system no longer follows the simple two-body decay laws from Eqs. (32) and (34). The reason is that at low densities there are enough correlations that we can no longer use the simple models from Section III C. In this regime ofhe small densities, a coarse grained description becomes approximately equivalent to the Lieb-Liniger model [12] but with inelastic interactions. As we have shown elsewhere [6], the decay at long times is then expected to follow the law dn/dt ∝n 4 .

V. DETAILED CALCULATIONS
The goal of this work is to find an effective model for the particles in the lattice, which works in the limit of fast dissipation, J hΓ. Our main tool to understand this limit is a generalization of Kato perturbation theory [13], also known as adiabatic elimination, to the superoperators V and L int . Section V A explains how our calculations relate to this broader scope. Readers less interested in this aspect may skip to Section V B where the actual derivation begins.

A. General ideas
Both V and L int are linear operators that act on an appropriate space of matrices. Even though within this space the superoperators are not Hermitian, the dominant term L int has infinitely many eigenstates forming a discrete spectrum of well separated points beginning at λ 0 = 0 and spanning through the left half of the complex plane [ Fig. 5b]. Each of these points represents a family of matrices that, in the absence of hopping term, decay at a rate given by the real part of the eigenvalue, Re(λ i ) ≤ 0. In the limit of strong dissipation, the hopping V will couple these eigenspaces very weakly, with an amplitude J/h much smaller than the typical eigenvalue separation, Γ.
Precisely in this limit of weak hopping we will be able to write analytic expansions of the eigenstates, eigenvalues and evolution equations for the perturbed superoperator (V + L int ). The ideas is to use Kato's resolvent method [13] with an expansion of the Liouvillian which uses a complete set of pseudo-projector operators In ordinary perturbation theory of Hermitian operators we find eigenspaces with eigenvalues which are well separated along the real axis. A weak Hermitian coupling, smaller than the energy separation between spaces ε |λi − λj|, only causes small shifts in these energy levels. (b) Our Liouville operator Lint is not Hermitian but has well separated eigenvalues in the space of density matrices (We have set Ur = 0 for simplicity). The real part of these eigenvalues now represents the decay rate of such matrices. The hopping of particles J acts as a weak perturbation J hΓ, that slightly changes both the real and the imaginary part of the original eigenvalues. In particular, the states with one particle per site, ρ0, now acquire contributions of O(J/Γ) 1 and O(J/Γ) 2,3,... by coupling directly, ρ1, or indirectly, ρ2,3,..., to faster decaying subspaces. Both in (a) and (b), for second order approximations of the unperturbed eigenspaces, λ0, we may neglect all spaces with an indirect coupling.
These operators, built from the right and left eigenvectors of L int , are non-Hermitian but this fact poses no difficulties in generalizing the usual perturbation theory.
In particular, the space of density matrices that had zero decay rate for J = 0 will transform into the perturbed eigenstates and eigenvalues The decay rate will no longer be zero but remain small and these states will acquire a small admixture of the matrices belonging to the unperturbed eigenspaces that had larger decay rates Notice the weight of the different contributions ρ n depends on how many applications of V we have to perform to connect ρ 0 to these other eigenspaces. Furthermore, we understand that it is the coupling of ρ 0 to the decaying states through V what makes the initially stable space lose particles. A similar reasoning can be applied to the rapidly decaying eigenspaces,ρ 1,2,... (J), which will be modified by the coupling V. However in this case the decay rates will remain of order Γ with small corrections from the hopping of particles. We therefore conclude that if we prepare our particles in any initial state, after a short transient t ∼ 1/Γ, most of the state will be captured by the eigenspace with the lowest decay rate (38). Following the structure given in Eq. (39) we will decompose our evolved state in a series of contributions from the unperturbed eigenspaces where, according to Eq. (39), the high order contributions have a vanishingly small size We are therefore allowed to perform a self-consistent approximation which consists on, first, writing the evolution equations for each subspace ρ n (t), then impose that all contributions ρ n≥2 0 and finally integrate the remaining equations until we reach an effective model for the lowest order term.

B. Local Projections
As sketched before, our study of the losses in the lattice begins by finding out the eigenspaces of the unperturbed superoperator, L int . This task is greatly simplified by the fact that L int = k L loc,k is a sum of commuting local superoperators L loc . We can thus focus on diagonalizing one of these superoperators on a single lattice site.
We will now introduce some notation. Since we are interested in density matrices as elements of a Hilbert space on which the superoperators act, we will introduce a basis of vectors of this space. For a single site we will define the Fock states and introduce the Hermitian conjugate (n , m | given by (n , m |n, m) = δ nn δ mm . In this basis L loc becomes a bidiagonal, nonsymmetric operator where ξ x = x(x − 1). The kernel of this operator is given by states with 0 or one particle per site, so that we can write the space of non-decaying density matrices as in Eq. (12). In addition to the right eigenvectors, (L loc − λ n )|v n ) = 0, we will also search the left eigenvectors, (w n |(L loc − λ n ) = 0, which have the property (w n |v m ) = δ nm . With these families of operators we will construct a set of pseudo-projector operators of the form We face one problem, though, which is that L loc acts on an infinite-dimensional space, with occupation numbers that can be arbitrarily large. We argue that for our purposes it suffices to truncate the Hilbert space to occupations n, m ≤ 2. The reason is that our initial states will all belong to the space ρ 0 given above and, since we will neglect contributions from third and subsequent order couplings [i. e. ρ 2,3,... = 0 in Eq. (40)], we will obtain at most a double occupation per site. The self-consistency of our approximation will be evident at the end.
Using the previous truncation, L loc becomes a 9 × 9 block-diagonal and banded matrix with projections onto eigenspaces and corresponding eigenvalues (46d) P loc 1a and P loc 1b are connected to P loc 0 by a single application of the hopping term, V, and P loc 2 by two applications of V. Note also that P loc 0 contains not only a projector onto single or zero occupancy, but also a term that destroys two particles at a site, emptying it. This will be essential later on.

C. Adiabatic elimination
As explained in Section III, we will consider that our states can be described by a density matrix with contributions coming from the eigenspace with the lowest decay rate, ρ 0 , and the two eigenspaces connected to it, ρ 1a,1b . These contributions are respectively obtained by applying the following pseudo-projectors onto ρ(t) P 0 = P loc 0 ⊗ · · · ⊗ P loc 0 , (47a) Connecting to the previous notation of projectors onto states with zero and one particles per site, it will be useful to realize that the zeroth order projector can be written as follows Since we will neglect higher order couplings, we can use Eqs. (10a) and (36) to write evolution equations for the density matrices in the form We have abbreviated V ij = P i VP j and introduced the notation that the index c runs through {1a, 1b}. The terms ρ c can be integrated out of the model using the fact that our states are initially prepared in the slow decaying manifold ρ c (0) = 0. Formal integration of Eq. (49b) yields which after integration by parts becomes We neglect the remaining integral, because it is of higher order in J/hΓ than the previous term, which is evident from the fact that dρ 0 /dt ∝ J in Eq. (49a). Insertion of ρ c (t) into Eq. (49a) yields The first line of this equation represents our effective model (16) and will be discussed in the following section. The second line is a transient that decays at a rate ∝ Γ. Therefore we obtained that the system converges to the slow-decaying eigenspace in a time t ∼ 1/Γ, much shorter than the typical time scaleh/J at which the effective model operates.

D. Hard-core bosons
Let us analyze the lowest order contribution to our effective model (52), given by L 1 = V 00 . Using the expression in Eq. (48), we obtain In other words, this Liouville operator is equivalent to a Hamiltonian in which we have projected out all states with double occupation. This is the hard-core bosons or Tonks-Girardeau gas model presented in Sect. III B.

E. Second order losses
We are now going to consider the second order Liouvillian L 2 from Eq. (16) We can expand this expression Using the property ρ 0 = Q 0 ρ 0 Q 0 , one realizes that the only relevant terms are L 2 ρ 0 = P 0 Aρ 0 /h 2 with We now consider the final projection with P 0 . Following Eq. (48), this pseudoprojector contains two operations: the first one keeps terms proportional to Q 0 H J Q 1 H J Q 0 , while the second one acts on the terms that create a doubly occupied site on each side of the density matrix, that is Q 1 H J Q 0 ρ 0 Q 0 H J Q 1 . Introducing T = Q 1 H J Q 0 /(−J), It is now time to rewrite everything in terms of hard core boson operators. We notice the following equivalence This arises from the fact that Q 1 projects onto a state with a single pair. Therefore, the tunneling term only contributes with processes that take two neighboring particles (C k ) and create a pair in one of the sites. Notice that C k already enforces the projection Q 0 . Using this notation we can simplify our expressions even further thus arriving to the final model which is studied in Sect. III C.

VI. CONCLUSION
We have shown analytically and confirmed numerically that strong, inelastic interactions can induce a Tonks gas dynamics for a cloud of molecules trapped in an optical lattice. The particles act like hard-core bosons, with dissipation playing the role of a strong repulsion. This effective model is completed with a reduced loss rate, γ eff ∝ J 2 /Γ which is much slower than both the tunneling amplitude, J/h, and the original loss rate, Γ.
Even with the small losses, the state of the system can at all times be described as an incoherent mixture of strongly correlated Tonks gases with different total particle number. In this respect, being based on the idea of using dissipation to create strong correlations, our paper connects to recent works which suggest using dissipation to engineer states and phase transitions [14,15].
We acknowledge financial support of the German Excellence Initiative via the program Nanosystems Initiative Munich and of the Deutsche Forschungsgemeinschaft via SFB 631. J. J. G.-R. acknowledges financial support from the Ramon y Cajal Program of the M. E. C. and the projects FIS2006-04885 and CAM-UCM/910758.