Kinetics of Many-Body Reservoir Engineering

Recent advances illustrate the power of reservoir engineering in applications to many-body systems, such as quantum simulators based on superconducting circuits. We present a framework based on kinetic equations and noise spectra that can be used to understand both the transient and long-time behavior of many particles coupled to an engineered reservoir in a number-conserving way. For the example of a bosonic array, we show that the non-equilibrium steady state can be expressed, in a wide parameter regime, in terms of a modified Bose-Einstein distribution with an energy-dependent temperature.

First experimental realizations of many-body reservoir engineering are starting to appear: it has been used to stabilize a circuit QED based Mott insulator in a 1D chain of eight transmon qubits against photon losses [19] as well as to dissipatively prepare quantum states in a three-transmon array [20].In that way, reservoir engineering contributes to the implementation of quantum simulators, especially in cases where the naturally available dissipation would not drive the system to the right many-body ground state.
In this context, one very important scenario concerns the case where the coupling to the reservoir conserves the total number of particles [20].Only then is it possible to cool without ending up in the vacuum state.Here, we are interested in providing a general framework to quantitatively describe the kinetics of many particles being scattered among states thanks to the (weak) interaction with an arbitrary non-equilibrium reservoir.The nonequilibrium nature is significant: even the steady-state distribution will depend on details of the interaction and of the reservoir noise spectrum, in contrast to the case of a thermal heat bath encountered, e.g., in certain approaches on particle-conserving photon equilibration [21][22][23].
We will introduce a perturbative approach to derive the steady-state distribution in momentum space.For the bosonic case treated here, this turns out to be a "deformed" Bose-Einstein distribution with an energydependent effective temperature.Illustrating the general theory in the specific case of a linear array of (harmonic) bosonic modes, we furthermore observe features such as negative-temperature states and prominent accumulation of particles at certain momenta during the time-evolution.The physics we encounter is partially reminiscent of cavity optomechanics [24], except that we are now dealing with a many-particle system instead of a mechanical resonator.Furthermore, the coupling to that system has been engineered to preserve the number of excitations, somewhat analogous to the unconventional quadratic coupling encountered in some optomechanical systems.Model-For concreteness we focus on bosons in a 1D array of modes with nearest-neighbor hopping.This physical situation can be realized, e.g., using linear elements of the superconducting circuits toolbox (capacitively coupled cavities).Each mode is further coupled to a lossy driven cavity via a density-density interaction [see Fig. 1 (a)].The driven, lossy cavities will act as an engineered non-equilibrium reservoir that mediates transitions between the eigenstates of the array while simultaneously conserving the number of excitations; this ensures a well-defined chemical potential.We note that this type of coupling between system and reservoir can be achieved, e.g., via the anharmonic transmon terms when considering a single cavity coupled dispersively to an array of transmons.That was realized in Ref. [20], using the reservoir to cool a Bose-Hubbard chain in the atomic limit, i.e., a three-site chain with at most two excitations.An alternative implementation would be an optomechanical array consisting of a chain of mechanical modes, each of which is coupled quadratically [25] to a driven cavity.
The Hamiltonian describing the system is given by where Ĥarray = ω 0 ) with ω 0 the frequency of the modes, J the hopping strength, and bj the annihilation operator of mode j.For the lossy cavities, the Hamiltonian is Ĥcav = ω c j â † j âj , with ω c the frequency, and âj the annihilation operator of a photon.Most importantly, we assume a densitydensity coupling between the array and the reservoir, of the form Ĥint = χ j â † j âj b † j bj , with χ the coupling constant.As we show later, such a coupling allows one to cool the array while conserving the initial number of excitations present in the system.Alternatively, one can consider a situation where only the first site of the array is coupled to a driven reservoir.This setting leads qualitatively to the same cooling dynamics.However, the situation studied here presents the advantage of being homogeneous, having a well-defined thermodynamic limit, and the cooling power scaling with the length of the array.Finally, the classical drive on the reservoir cavity is given by Ĥdrive = Ω 0 j exp(−iω d t)â † j + H.c., with ω d the driving frequency and Ω 0 the amplitude of the driving field.
The cooling dynamics of the array is best understood by diagonalizing Ĥarray and switching to a frame rotating at the drive frequency ω d .To diagonalize the Hamiltonian of the array, we introduce the mode operators βk defined via the relation bj = k ϕ k (j) βk where ϕ k (j) is the mode function determined via the boundary conditions of the problem (periodic or open).Then the interaction Hamiltonian is given by from which one sees that a transition from mode k to k is possible without exchanging particles with the reservoir.We note that for a finite-size system the wave vector k has a discrete set of values, which are different for open and periodic boundary conditions.
To understand how such an interaction can give rise to cooling dynamics, it is useful to displace the reservoir field, i.e. âj → a 0 + dj , where a 0 is the classical amplitude of the field and the operator dj describes the quantum noise affecting reservoir j [26,27].Here, we neglect the fast transient dynamics and consider the reservoir to be in the steady state.This leads to a 0 = iΩ 0 /(i∆ − κ/2) with κ the energy decay rate.In the displaced frame the Hamiltonian of the system is where ω k = ω 0 + 2J cos(k) and ∆ = ω d − ω c is the detuning between the reservoir cavity and the drive.Using Eq. ( 3) one can derive Golden rule rates describing transitions from mode k to k , which generalizes optomechanical cooling ideas [26,27] to the many-body case, We introduce the power spectrum of the non-equilibrium reservoir noise, ]. Now we are in a position to describe the cooling dynamics of the array, by formulating a set of coupled kinetic equations that describe the evolution of the average occupation of the eigenmodes of the array: with n(ω k , t) = β † k (t) βk (t) .In the presence of weak interactions between particles, one would have to supplement Eq. ( 5) by addition of two-particle scattering terms [28].
As it can be seen from Eq. ( 4), the sign of the detuning ∆ determines whether the noise S(ω) [see Eq. ( 4)] peaks at either negative (∆ > 0) or positive (∆ < 0) frequency.As a consequence, by choosing the detuning ∆ to be negative and assuming ω k > ω k , we have Γ k→k > Γ k →k in Eq. ( 5).This corresponds to particles in high-energy modes being preferentially scattered into low-energy modes by having the reservoir absorb the excess of energy, leading to cooling.The reverse situation, with positive detuning ∆, means the noise pumps energy incoherently into the array, scattering particles to higher energies.Another property of Eq. ( 4) that influences the dynamics described by Eq. ( 5) is the Bose enhancement factor; the rate at which a boson is scattered to a state with occupation n is enhanced by a factor n + 1. Results -In the long-time limit, the system approaches a steady-state distribution.In general, this does not correspond to any thermal equilibrium distribution, as it can be seen in our numerical results [see Fig. 2 (a)] obtained on the basis of Eq. ( 5).An important aim in any non-equilibrium scenario is to characterize the resulting steady-state distributions.
First insights can be obtained by recalling that even for a noise source that is not in thermal equilibrium one can always define an effective temperature associated to a single transition frequency by using the Stokes relation In the limit where ∆ 2 + (κ/2) 2 ω 2 we can expand S(±ω) in powers of 1/[∆ 2 + (κ/2) 2 ].We also expand tanh[β eff ω/2] around β eff ω → 0.Then, we find that the effective inverse temperature for a single transition frequency is given by where we defined the equilibrium temperature Eq. ( 6) shows that all transitions have the same temperature β 0 .In this limit, the reservoir effectively acts as a thermal reservoir and we expect the steady state to be described by Bose-Einstein statistics with an inverse temperature β 0 and a chemical potential µ: We will now demonstrate that, in the vicinity of this limiting case, the steady-state distribution can be understood as a certain deformed version of the Bose-Einstein distribution; in other words we look for a description in terms of small corrections around the equilibrium statistics.To do so, we turn to a perturbative solution of Eq. (5).
The perturbation theory is best set up by considering the continuum limit of Eq. ( 5) for the steady state.We have with ω k = ω k + ω [see Fig. 1 (b)].Here we introduced the density of states (for the 1D array, The idea of the perturbative treatment is to look for a series representation of the distribution n(ω k ): n(ω k ) = i n i (ω k ), with n i ∝ (β 0 ω k ) i (see supplemental information for more details).We find where n BE (ω k ) is the Bose-Einstein statistics and c is a constant of integration.We can fix both µ and c in Eq. ( 8) by enforcing particle number conservation,   10)] comparing the steady-state solution with the perturbative solution.In the limit β0J 1, the perturbative solution describes accurately the steady state of the system.For these simulations, we considered the experimentally relevant situation of open-boundary conditions and the parameters L = 100, N = 50, χ = J/1000, and Ω0 = J/10, unless specified otherwise.
k n(ω k ) = N ; we first obtain µ by imposing normalization for n BE and then extract c from the same condition applied to the next order.
When we compare Eq. ( 8) to a generalized Bose-Einstein statistics with an energy-dependent inverse temperature (setting β 0 → β(ω k )), we conclude, to leading order: In Fig. 2 (a), we display the inverse temperature as a function of energy.One can see that the curvature depends on the ratio between the detuning ∆ and the decay rate κ.We note that Eq. ( 9) in addition to reduce to β 0 for a small expansion parameter [see Eq. ( 6)] also reduces to β 0 at the special point ∆ = − √ 3κ/2.To assess the accuracy of our perturbative solution, we quantify its distance from the true steady-state solution of Eq. ( 8).We employ the Kullback-Leibler divergence [29] (relative entropy) which describes the amount of information lost when approximating one distribution by another.Here, we expressed our perturbative solution in terms of the normalized density p(ω k ) = n(ω k )/N [with N the total number of particles and n(ω k ) given by Eq. ( 8)] and compare it to the steady-state distribution p Figure 2 (c) shows the Kullback-Leibler divergence as a function of detuning and decay rate.The numerical results were obtained by time-evolving the kinetic equations until the system settles into the steady state.
We observe the pertubative solution in Eq. ( 8) approximating the non-equilibrium steady state well in a wide range of parameters.The approximation becomes better, as expected, in the limit β 0 J 1 (for β 0 J → 0, the effective temperature becomes energy-independent, and Eq. ( 8) reduces back to the Bose-Einstein distribution).We stress that Eq. ( 8) is much closer to the real solution than the Bose-Einstein distribution itself; indeed, in a large parameter regime we have While we focussed on the specific case of a onedimensional chain coupled to a non-equilibrium reservoir with a particular noise spectrum, our perturbative treatment is applicable to any spectrum and for noninteracting bosons in any dimension, in arbitrary potentials, with suitably modified mode functions and density of states.In all of these cases, it can serve to extract the steady-state non-equilibrium distribution brought about by reservoir-induced particle scattering.
Equation ( 8) represents our main analytical result.We now discuss some distinct properties of the nonequilibrium distribution.
In thermal equilibrium, detailed balance ensures that the Bose-Einstein distribution is independent of the density of states.This is no longer true here, out of equilibrium.As we show in Fig. 2 (b), signatures of the density of states (with its characteristic divergence at the band edge in 1D) can be observed when the noise spectrum itself [see Eq. ( 4)] is sufficiently narrow i.e, for κ/J 1.In that case, prominent features (almost non-analytic) appear in the distribution at energies separated from the band edge by the detuning.These features are even more pronounced during the temporal dynamics (see below).
When deriving the kinetic equation [Eq.( 4)], we pointed out the presence of the bosonic enhancement factors.The effects of such an enhancement can be observed when the particle density increases, particularly in the most strongly occupied regions of the band [see Fig. 3 (a)].For red detuning, we obtain many-particle ground state cooling.To assess its efficiency, we plot in Fig. 3 (b) the difference between the occupation of the single-particle ground state and the highest excited state in the band, δn = n high −n GS .The bosonic enhancement serves to increase the efficiency of ground state cooling.
In stark contrast to optomechanics, where the spectrum of the mechanical resonator is unbounded, we are dealing with a physical system that possesses a bound (many-particle) spectrum.This has an important consequence: even for blue detuning (∆ > 0), and without any recourse to non-linear cavity dynamics, the system is still stable, as can be seen in Fig. 3 (c).In this "heating" regime, particles accumulate at the upper band edge.Such states can be described by a negative temperature [see Fig. 3 (d)].Similar negative-temperature states were previously obtained using localized spin systems [30][31][32] and ultracold atoms trapped in optical lattices [33].Creating stable negative states with optical lattices is challenging as it requires relatively complex state preparation.By contrast, here, one only needs to choose a positive detuning.
Finally, we investigate the time-evolution.We are interested in describing the mechanism that leads to the accumulation and depletion of particles in the distribution, as shown in Fig. 4. As we noted earlier, these features are present in the steady-state distribution [see Fig. 2 (b)], so we expect them to be generic also during the evolution.To consider the simplest possible situation, we assume the initial distribution to be uniform.Such a uniform distribution in k-space can be realized, e.g., by incoherently loading the whole array with bosons at constant density (or quenching from a Mott insulator).At early times, the incipient deviations from the uniform distribution can be obtained perturbatively:  where we have replaced n(ω k , t) by its initial uniform value n 0 .The integral in Eq. ( 11) is maximal (or minimal) when the divergence of the density of states aligns with the maximum of the spectrum S(ω) [or S(−ω)].
Physically, this translates into having a large incoming rate of particles for states with ω k = 2J + ∆ and a large outgoing rate for states with ω k = −2J − ∆ (assuming ∆ < 0).As a consequence, particles accumulate (deplete) at ±∆ from the band edges [see Figs.In experimental implementations, two factors may come into play to distort the long-term dynamics from the idealized scenario -we mention them briefly, though a detailed study would be beyond the scope of the present work.First, there is the loss of array bosons at a constant rate (e.g.intrinsic cavity decay or mechanical dissipation).This will simply lead to an overall depletion of the total particle number N (t) ∼ e −Γt (i.e. the steady state will slowly evolve, due to the modification of bosonic enhancement factors).Second, if there are interactions between the bosons (e.g.attractive interactions for a chain of transmon qubits instead of linear cavities), these will lead to additional scattering that is compatible with ground-state cooling but will drive the distribution closer to thermal.Conclusion -We have investigated the scattering kinetics of bosonic particles subject to what is probably the most widely used engineered reservoir, namely a driven lossy cavity.We have shown that the resulting nonequilibrium steady-state distribution can be perturbatively analyzed in a wide parameter regime and interpreted as a deformed Bose-Einstein distribution with an energy-dependent temperature.In this work, we have focussed on the illustrative case of a one-dimensional array, with each site coupled to a driven cavity.However, our treatment remains valid and we expect similar observations in general for non-equilibrium reservoirs coupled in a particle-conserving manner to arbitrary non-interacting many-body systems, i.e. in higher-dimensions, for arbitrary lattices and band structures, and also for fermionic systems, with straightforward modifications.

Figure 1 .
Figure 1.(Color Online) (a) Schematic representation of an array of bosonic modes, where each site experiences densitydensity coupling to a non-equilibrium reservoir (a driven, lossy cavity).(b) A particle at energy ω1 is scattered to a state of energy ω2 < ω1 by emitting a photon into the reservoir.The rate at which particles scatter is proportional to the noise spectrum S(ω) of the reservoir.

Figure 3 .
Figure 3. (Color Online) Bosonic distribution in momentum space (long-time limit).(a) Impact of the Bose enhancement factor (∆ = −J, κ = 3J).As the particle density increases, the low-energy levels become more occupied during the cooling dynamics.(b) Efficiency of ground state cooling.We plot the difference between the average occupation of the ground state and the highest excited state, δn = nGS−n high (κ = 3J).(c) Transition towards negative temperature.For blue detuning (∆ > 0), the system exhibits stable negative-temperatures steady states.(d) Schematic representation of states with positive and negative temperatures.