Dissipation-induced correlations in 1D bosonic systems

The quantum dynamics of interacting bosons in a one-dimensional system is investigated numerically. We consider dissipative and conservative two-particle interactions, and integrate the master equation describing the system dynamics via a time-evolving block-decimation (TEBD) algorithm. Our numerical simulations directly apply to stationary-light polaritons in systems where atoms and photons are confined to the hollow core of a photonic crystal fibre. We show that a two-particle loss term can drive an initially uncorrelated state into a regime where correlations effectively inhibit the dissipation of particles. The correlations induced by two-particle losses are compared with those generated by an elastic repulsion. For the considered time range, we find a similar behaviour in local density-density correlations but differences in non-local correlations.


Introduction
The preparation of quantum systems in entangled many-body states is an indispensable resource for quantum information science [1] and allows one to access the rich physics of strongly correlated many-body systems [2]. Despite its reputation as the main adversary of quantum information schemes, dissipation was recognised recently as a versatile tool for the generation of specific many-body states [3,4,5].
The preparation of correlated quantum states via dissipation requires that the target state itself is metastable. An intriguing experiment [6] with cold molecules showed that a dissipative two-particle interaction between bosons in 1D gives rise to correlations that inhibit the loss of particles. The formal analysis of this effect is based on a generalised Lieb-Liniger model [7,8,9] and suggests that the two-particle loss term effectively results in a repulsion such that inelastic collisions between particles are avoided. In the limit of strong interactions [8], the two-particle losses give rise to a Tonks-Girardeau gas [10] where bosons behave with respect to many observables as if they were fermions.
Recently, it was shown that 1D systems of interacting bosons can be realised via dark-state polaritons that arise in light-matter interactions under conditions of electromagnetically induced transparency [11,12,13,14,15]. A possible experimental realisation of the polariton schemes is the setup described in [16,17], where photons and atoms are simultaneously confined to the hollow core of a photonic-crystal fibre, see Fig. 1. Another approach comprises the experimental setup in [18], where atoms are trapped near the surface of an optical nanofibre using evanescent waves. The nature of the polariton-polariton interaction in these approaches can be tuned by external parameters and can be elastic [13] or dissipative [14,15], and the interaction strength is maximal for a purely dissipative interaction [14]. Polariton schemes [13,14,15,19,20,21] are extremely promising candidates for the realisation of strongly correlated photon states.
Here we consider bosons in 1D and investigate the quantum dynamics after a sudden switch-on of elastic and inelastic interactions via time-dependent density matrix renormalisation group techniques. The master equation of the system and the uncorrelated initial state are chosen such that our simulations directly apply to polariton systems [13,14,15]. All details of our model are discussed in Sec. 2, and numerical results are presented in Sec. 3. In particular, we are interested in the preparation of a correlated, non-decaying state via dissipative two-particle interactions that we discuss in Sec. 3.1. Non-local correlations are investigated in Sec. 3.2, where we show that dissipation and repulsion give rise to different results. Finally, a summary is provided in Sec. 4.

Model
We consider bosons of mass m in a one-dimensional system of length L that experience a two-particle contact interaction that can be either conservative, dissipative or a mixture of both. First we describe the full master equation governing the time evolution of the system in Sec. 2.1. A physical realisation of this model via interacting dark-state polaritons [14,15]

Master equation
The dynamics of bosons in our 1D system is described by the master equation where H eff is a non-hermitian Hamiltonian, Here ψ(z) is the field operator obeying bosonic commutation relations and the first term in Eq. (2) describes the kinetic energy. The parameterg is the complex coupling constant and its real part gives rise to a hermitian contribution to H eff in Eq. (2) that accounts for elastic two-particle collisions. On the other hand, the imaginary part ofg together with the term in Eq. (1) constitutes a two-particle loss term that can be written in Lindblad form [22] as The last term in Eq. (1) is defined as and causes diffusion of the field ψ. Throughout this paper we choose such that the losses associated with the diffusion term are small on a timescale set by the kinetic energy term in Eq. (2). The master equation (1) without the last term L D ̺ defined in Eq. (5) can be identified with the dissipative Lieb-Liniger model [8]. For a fixed number of particles N ph , essential features of the system are described by the complex Lieb-Liniger parameter The spectrum of the effective Hamiltonian in Eq. (2) was analysed in [8] via the Bethe-Ansatz. It was found that the ground state of the so-called gaseous states (all momenta converge to finite values for |G| → ∞) approaches the Tonks-Girardeau gas in the limit |G| → ∞. This state is metastable since in the Tonks-Girardeau gas limit two particles never occupy the same position in space, and hence the two-particle contact interaction does not contribute to the dissipation of particles. This feature can be formally derived from the master equation (1). The two-particle loss term Im(g/2)D[ψ 2 ] in Eq. (1) causes a time dependence of the density of particles according to ∂ t n(z) = 2 Im(g)g (2) (z, z) n(z) 2 , wheren(z) = ψ † (z)ψ(z) is the particle density operator and is the second order correlation function. Note that the kinetic energy and the diffusion term L D ̺ in Eq. (1) will, in general, induce an additional time dependence of n(z) . The kinetic energy term conserves the total number of particles but gives rise to polariton fluxes from and to position z that would redistribute the polariton density. In regions where the system and the particle density are homogeneous, these fluxes can safely be neglected. The impact of the diffusion term L D ̺ is negligible for the considered time scales since we chose a small value for the parameter D in Eq. (6). In the strongly correlated regime and for a homogeneous system, the gaseous ground state of H eff obeys [8] g (2) (z, z) and hence g (2) (z, z) ≪ 1 for |G| sufficiently large. The combination of Eqs. (8) and (10) implies that a purely dissipative interaction term supports metastable states.
γ32 γ31 Figure 1. (a) Considered setup of N a atoms confined to an interaction volume of length L and transverse area A. Ω ± are the Rabi frequencies of the classical control fields, andÊ ± are the quantum probe fields. (b) Atomic level scheme. γ ij is the full decay rate on the |i ↔ |j transition, δ and ∆ label the detuning of the probe fields with states |3 and |4 , respectively, and ε is the two-photon detuning.

Physical realisation
An example for a physical system that is described by the master equation (1) is shown in Fig. 1(a), where photons and atoms are simultaneously confined to the hollow core of a photonic-crystal fibre. Since the light-guiding core of the optical fibre is of the same order of magnitude as the optical wavelength, the fibre represents a onedimensional waveguide for the optical fields. The level scheme of each atom inside the fibre is shown in Fig. 1(b). The dynamics of the system can be described in terms of dark-state polaritons [23] representing bosonic quasi-particles that arise in light-matter interactions under conditions of electromagnetically induced transparency (EIT) [24]. The two counter-propagating control fields Ω ± [see Fig. 1(a)] of equal intensity give rise to stationary light [25,26] where the polaritons experience a quadratic dispersion relation like bosons in free space. Moreover, the coupling of the probe fields to the transition |2 ↔ |4 induces a two-particle contact interaction between the polaritons that can be conservative [11,13], dissipative or a mixture of both [14,15]. It can be shown [14,15] that the dark-state polariton dynamics can be described by the master equation (1), where the effective mass m and the complex coupling constantg depend on the detunings δ, ∆ [see Fig. 1(b)] and the intensity of the control fields. In particular, we point out that a proper adjustment of the detuning δ allows one to fulfil Eq. (6), and the choice of the detuning ∆ determines the relative contribution of elastic and inelastic processes to the two-particle interaction.
In an experiment, the evolution of the polaritons under Eq. (1) has to be preceded by a loading process that could be realised as follows [13]. Initially the probe field modes are empty and only the control field Ω + is switched on. Then a probe pulse copropagating with Ω + enters the fibre under conditions of electromagnetically induced transparency. Here we assume that the transition |4 ↔ |2 is far off-resonant during the loading process such that its influence is negligible. The dynamics of the probe field inside the medium can be described concisely in terms of dark-state polaritons [23]. As compared to the pulse propagating in vacuum, the polariton pulse inside the medium is spatially compressed and travels with a reduced group velocity v g ∝ Ω 2 + . Once the pulse is entirely inside the fibre, the group velocity of the polariton pulse is reduced to zero by an adiabatic switch-off of the control field Ω + . The latter process maps the polaritons into a stationary spin coherence, and all probe field modes are empty. Eventually both control fields are switched on adiabatically and with equal intensity, Ω ± = Ω 0 . The time evolution of the corresponding dark-state polaritons is then determined by the master equation (1) [14,15]. Note that the adiabatic switch-on can be much faster than all typical timescales in the master equation (1). At this stage, the detuning ∆ of the probe field with the |4 ↔ |2 transition can be adjusted by a proper choice of the frequency of the counterpropagating control fields. This frequency can be different from the one used for the single control field during the loading process.
The correlations that built up under the evolution of the master equation (1) can be measured if the intensity of the control field Ω − is reduced. The latter process converts the stationary polaritons into a propagating pulse with controllable group velocity. This procedure maps spatial correlations of the trapped pulse into temporal correlations of the output light that can be detected via standard quantum optical techniques.

Discretised model
The dynamics of the system is studied via the time-evolving block decimation algorithm [27,28,29] for density operators. The discretised version [30] of Eq. (1) is given by a generalised Bose-Hubbard model augmented by loss terms, Here H 1 and H 2 correspond to the kinetic energy term and the conservative contribution to the two-particle interaction, respectively, and J and U are defined as The parameter ∆z is the lattice constant of the discretisation grid. The terms L 1 ̺ and L 2 ̺ in Eq. (1) are the discretised versions of the diffusion term in Eq. (5) and the dissipative contribution to the two-particle interaction, respectively, and The discretised version Eq. (11) is a good approximation of the continuous master equation (1) if the smallest wavelength λ min involved is much larger than the grid spacing ∆z. Initially λ min is determined by the momentum components of the initial state, but the evolution under Eq. (11) will change the momentum distribution of the sample. In particular, the suddenly switched-on two-particle repulsion term leads (among other processes) to a population of higher momentum states. It was found [31] that the lattice model is a good approximation provided that U/J ≪ 1. If the latter inequality is violated, lattice artefacts occur. We find that these effects do not occur for a purely dissipative two-particle interaction (U = 0) even if Γ 2 /J ≈ 1. For all numerical simulations we consider a grid with N s = 500 sites, and the physical state space at each site l is spanned by the four Fock states {|0 l , |1 l , |2 l , |3 l } allowing for a maximum occupation of three particles at each site. We verified the convergence of our algorithm by varying the bond dimension χ of the matrix product state. In our simulations, the truncation errors ǫ trunc are bounded by ǫ trunc ≤ 10 −4 for χ = 100.
The choice of the initial state for the numerical integration is motivated by the loading process of an optical pulse into the fibre as described in Sec. 2.2. We consider an initial state with a mean number of N 0 = 5 particles and a spatial density distribution as shown in Fig. 2. Here we neglect distortions and correlations due to the optical nonlinearity that may build up during the loading process. The initial shape of the polariton pulse is then determined by the slowly varying envelope of the probe pulse entering the system [23]. Furthermore, we assume that the initial polariton pulse is derived from laser light and thus described by a product of coherent states. In the discretised model we thus model the initial state by a coherent state at each site,  where l |c l | 2 = 5 and |c l | 2 = ∆z n(z l ) with z l = l∆z. We point out that |ψ 0 in Eq. (17) is an uncorrelated state with g (2) (z, z ′ ) = 1 for all values of z, z ′ in the sample.

Numerical results
Next we describe the results of the numerical integration of the master equation (11). Starting from the initial state in Eq. (17), the interaction is suddenly switched on at t = 0. We consider various scenarios where the interaction is purely dissipative or predominantly conservative. In Sec. 3.1, we discuss whether two-particle losses alone are able to drive the system into a correlated regime where dissipation of particles is suppressed. As discussed in Sec. 2.1, the inhibition of two-particle losses is linked to local correlations in the sample. The second part Sec. 3.2 is concerned with nonlocal correlations. In particular, we compare non-local correlations created by a purely dissipative and a predominantly conservative interaction.

Inhibition of two-particle losses
The two-particle decay term in Eq. (1) in general leads to a loss of particles from the system. The rate at which the particle density diminishes is determined by Eq. (8) and depends crucially on the local correlations g (2) (z, z) of the system. Figure 3 shows the temporal evolution of g (2) (z 0 , z 0 ) at the centre z 0 = L/2 of the sample for different values of the Lieb-Liniger parameter |G| according to the numerical integration of Eq. (11). The values of |G| correspond to the mean particle number N 0 = 5 of the initial state, see Eq. (17). All solid lined lines in Fig. 3 belong to a purely dissipative two-particle interaction, Re(g) = 0. On the contrary, the dashed lines in Fig. 3 (a) and (b) correspond to a predominantly conservative interaction |Re(g)/Im(g)| = 10.
In the case of a purely conservative interaction (Im(g) = 0), it was shown [31] that local quantities reach a stationary state on a timescale T loc = /(gρ), where ρ = N ph /L is the density of particles. Note that the Lieb-Liniger gas does not thermalize globally [32].
Here we generalise this definition to the case of complexg and define Due to the enhanced numerical complexity of TEBD for density operators and numerical limitations, our simulations cannot cover the full time range T loc . In the following we scale time in the presentation of numerical results with The comparison between the dashed and solid lines in Figs. 3 (a) and (b) shows that dissipative and repulsive two-particle interactions are equally effective in the repulsion of particles on a timescale proportional to the inverse interaction strength |g|. Note that for the polariton system discussed in Sec. 2.2, |g| is at least one magnitude larger for a purely dissipative interaction [14,15] as compared to the conservative case. It follows that small values of g (2) (z 0 , z 0 ) are reached much faster in absolute time for a dissipative two-particle interaction. We emphasise that the minimal values for g (2) (z 0 , z 0 ) in Fig. 3 will decrease further in time until their equilibrium values are reached. An estimate of these equilibrium values can be obtained via the work presented in [31]. Here the authors performed TEBD simulations for the unitary time evolution of bosons in 1D with two-particle repulsion. In particular, this system could be described by a pure quantum state which reduces the numerical complexity as compared to our simulations of the full density operator. Therefore, the simulations in [31] could be carried out up to times t = T loc . Given that the initial dynamics of g (2) (z 0 , z 0 ) for two-particle losses and for conservative interactions is very similar, one may suppose that the equilibrium values of g (2) (z 0 , z 0 ) will also be almost the same. For this case we estimate that g (2) (z 0 , z 0 ) for our system should lie within the range [0.1, 0.2] for |G| = 10 and |G| = 20 at t = T loc .
Next we determine the temporal evolution of the particle density at the centre of the cloud via Eq. (8) and the numerical results for the temporal evolution of g (2) (z 0 , z 0 ). The result is represented by the solid lines in Fig. 4 and in good agreement with the direct evaluation of the particle density n(z) via the numerical integration of the master equation (11) (dotted lines). The small deviations are due to the diffusion term in Eq. (5) that results in particle losses and that were omitted in Eq. (8). Note that all curves in Fig. 4 correspond to a purely dissipative two-particle interaction, Re(g) = 0. In order to visualise the slowdown of two-particle losses due to the decrease of g (2) (z 0 , z 0 ) in time, the dashed line in Fig. 4 shows the temporal evolution of n(z) according to Eq. (8) with g (2) (z 0 , z 0 ) = 1 for all times. The slowdown of two-particle losses is most pronounced for |G| = 100 where a quasi-stationary regime is reached with g (2) (z 0 , z 0 ) close to zero. We find that the decrease of the mean number of particles N / N 0 follows the same curves as the particle density at the centre of the cloud. It follows that more than 90% of the initially present particles are left at t/τ c = 4 for |G| = 100. Figures 4 (b) and (c) correspond to |G| = 10 and |G| = 20, respectively, and demonstrate a clear slowdown of two-particle losses. On the basis of the estimates for the equilibrium values of g (2) (z 0 , z 0 ) obtained via the results presented in [31], we suppose that 50% (60%) of the initial particles are still present at t = T loc for |G| = 10 (|G| = 20). On the contrary, g (2) (z 0 , z 0 ) will not drop to values close to zero for |G| = 1 and hence we expect that a significant number of particles is lost at t = T loc .

Non-local correlations
The spatial dependence of g (2) (z 0 , z) as a function of z (z 0 = L/2 is the centre of the sample) is depicted in Fig. 5. Several snapshots in time are shown indicating that nonlocal correlations exhibit oscillations that propagate outwards. The oscillations for the longest evolution times in Fig. 5 exhibit several minima that lie below unity. It follows z/L z/L g (2) (z 0 , z) g (2) (z 0 , z) z/L z/L g (2) (z 0 , z) g (2) (z 0 , z) that a dissipative two-particle interaction alone to some extend gives rise to spatial order of the particles. Note that non-local correlations reach a quasi-stationary regime on a timescale that is much larger than T loc . In particular, we point out that the correlations shown in Fig. 6 are not Friedel oscillations [2] that characterise the Fermionic nature of particles in the Tonks-Girardeau regime. In Sec. 3.1 we concluded that repulsion and dissipation are equally effective in the creation of small values of g (2) (z, z). On the other hand, dissipative and conservative two-particle interactions give rise to different non-local correlations as shown in Fig. 6. The solid and dashed lines correspond to a purely dissipative and a mostly conservative two-particle interaction, respectively. It can be seen that the non-local correlations propagate at the same velocity in space in the dissipative and conservative case, but the oscillations are more pronounced for a repulsive two-particle interaction. In addition, two-particle losses give rise to a broader dip in g (2) (z 0 , z) as a function of z as compared to the conservative case. Note that the oscillation amplitudes in g (2) (z 0 , z) are also damped by the diffusion term in Eq. (5).

Discussion and summary
In this paper we investigated the quantum dynamics of bosons in 1D after a sudden switch-on of two-particle losses. The master equation and the uncorrelated initial state are chosen such that our simulations directly apply to stationary-light polariton systems [13,14,15]. Apart from a (small) diffusion term, the considered master equation coincides with the dissipative Lieb-Liniger model [8]. In particular, we investigated the dissipation-induced preparation of correlations that effectively inhibit two-particle losses. We found that a metastable regime where losses are substantially suppressed can be prepared for large values of the Lieb-Liniger parameter |G| ≫ 1. Formally, it was shown [7,8] that the two-particle loss term effectively results in a repulsion between particles such that they never occupy the same position in space. Therefore, dissipation via the two-particle contact interaction is inhibited. A physical explanation for this counterintuitive result can be given for a discrete lattice model [6,7]. There it can be regarded as a manifestation of the quantum Zeno effect, where the two-particle losses play the role of a continuous measurement that always projects the system onto states with less than two particles in each lattice site before higher particle numbers can build up. Furthermore, we compared the correlations induced by dissipative and repulsive interactions, respectively. We found that dissipation and repulsion are equally efficient for the generation of local correlations g (2) (z, z) ≤ 1 on a timescale proportional to the inverse interaction strength. On the other hand, non-local correlations induced by dissipation show different features than those created by two-particle repulsion. As compared to an elastic repulsion between the particles, we find that two-particle losses give rise to a broader dip in g (2) (z 0 , z) as a function of z (z 0 is the centre of the sample). Both repulsion and dissipation give rise to spatial oscillations in g (2) (z 0 , z) that spread in time and indicate to some extend a spatial ordering. We find that these oscillations are more pronounced for a repulsion between the particles as compared to the dissipative case.