Reservoir Engineering for Classical Nonlinear Fields

Reservoir engineering has become a prominent tool to control quantum systems. Recently, there have been first experiments applying it to many-body systems, especially with a view to engineer particle-conserving dissipation for quantum simulations using bosons. In this work, we explore the dissipative dynamics of these systems in the classical limit. We derive a general equation of motion capturing the effective nonlinear dissipation introduced by the bath and apply it to the special case of a Bose-Hubbard model, where it leads to an unconventional type of dissipative nonlinear Schr\"odinger equation. Building on that, we study the dynamics of one and two solitons in such a dissipative classical field theory.

In that context, an important and natural question to ask is whether the dissipation induced by the reservoir can be made to preserve particle number in the system.In many of the most relevant platforms for quantum simulation [20,21], this is nontrivial, since excitations of qubits or bosonic particles like photons and phonons are naturally destroyed by dissipation.At the same time, when particle number is conserved interesting complex ground states can be reached, e.g., forming a Bose-Einstein condensate or some highly entangled state.Theoretical proposals have investigated that issue in some detail [4,17,[22][23][24].
First successful steps have been taken, particularly in an experiment that demonstrated the cooling of a qubit chain to a Mott-insulating state [15].As a result, interesting questions have been raised regarding the kind of many-body states that can be reached using such an approach [25,26].Nonequilibrium Bose condensation of photons is another important direction making use of particle-conserving dissipation [27][28][29].
When we specialize to a concrete system, we will study the classical limit of a Bose-Hubbard chain where each site is coupled via its particle number to a driven and dis-Figure 1.(a) The physical scenario of a bosonic system (in the classical limit) coupled to a driven reservoir, in a particleconserving way.(b) Particular illustrative realization considered in the text, a Bose-Hubbard chain (blue bj) with sites coupled to driven cavities (orange bosonic modes aj) with coupling strength χ and to each other with hopping rate J.
See main text for details.(c) Spectrum of a dissipative cavity, where the incoming drive can be up-scattered in frequency through interaction with the bosonic many-body system, extracting energy.The detuning ∆ and the dissipation rate κ determine the properties of this driven reservoir.
sipative cavity (see Fig. 1).This can potentially be implemented using superconducting qubits tunnel-coupled to each other, simulating the Bose-Hubbard model [30,31], provided they are additionally coupled by cross-Kerr interactions to driven microwave resonators forming the reservoir.There are already works towards particleconserving interactions with engineered reservoirs, e.g., an experiment showing particle conserving autonomous cooling of a Bose-Hubbard chain [14] and a recent experiment preparing a Mott insulator in a Bose-Hubbard chain [15].Further works in that direction include an experiment showing Bose-Einstein condensation in an optical cavity [27] and theory work for parametric coupling to generate light with a chemical potential [23].The particular system we study in the present article has been considered previously by some of us in the quantum regime, deriving kinetic equations describing the scattering of particles by the coupling to the reservoir, see Ref. [24].
Since many of such systems used for quantum simulations aided by reservoir engineering are bosonic in nature, it is a meaningful and interesting question to ask how they would behave in the classical limit.Suppose we have particle-number-conserving dissipation, brought about by reservoir engineering, and acting on a bosonic many-body system (see Fig. 1): What are the effective classical equations of motion describing this scenario?Analyzing those may give us new physical insights that can then even be relevant for interpreting the original quantum dynamics.Moreover, in the classical limit of interacting bosonic systems, we can find collective excitations like solitons.How do those excitations behave in the presence of particle-number conserving dissipation?Can we identify how their motion described in terms of collective variables is modified in the presence of such unconventional dissipation?
In the present article, we set out to answer these questions.We will present three main results.First, in Sec.II we derive a general set of classical effective equations of motion for any bosonic many-body system coupled to a driven reservoir in a particle-conserving manner.These equations are still completely general, in that the Hamiltonian of the bosonic system can be arbitrary, but they already show the shape of the effective dissipation and allow us to analyze the dependence of the decay rate on parameters like the detuning between external drive and reservoir modes.The description found there can be viewed as a generalization of some formulas known from laser-cooling, e.g., in cavity optomechanics, to the case of arbitrary many-body systems, treated in the classical limit.As a second result, we specialize to a lattice of coupled bosonic modes in Sec.III, and we derive a nonlinear Schrödinger equation in the presence of particle-conserving dissipation (both on the lattice and in the continuum limit [Sec.III A]).In that equation the dissipation term is nonlinear, in contrast to the conventional dissipative nonlinear Schrödinger equation, where linear dissipation is considered.Finally, we apply these equations to describe the motion of solitons.There, as a third main result (Sec.III B), we show how the effective equations of motion for the collective coordinates of a soliton are modified by the new dissipation terms.We compare our analytical derivations with numerical simulations.Lastly, we briefly discuss the interaction of two solitons in Sec.III C and conclude with Sec.IV.

II. SYSTEM-RESERVOIR COUPLING WITH PARTICLE CONSERVATION: EFFECTIVE DISSIPATIVE EQUATION
We start out with a general situation, where a bosonic quantum many-body system is coupled to a reservoir in a particle-conserving manner.This will allow us to derive a formula of wide applicability, describing the effective dynamics of such a system, without yet specializing to a particular case.We consider an arbitrary bosonic lattice model.Each system site is coupled independently to a driven dissipative reservoir mode âl .Importantly, to ensure particle-number conservation, the coupling is of the density-density type.Thus, the full Hamiltonian is of the form Here, ĤS describes the system containing the bl modes, which we do not need to specify further at this point.The reservoir modes (â l ) and their interaction with the environment and external drives is contained in ĤB .Extensions like correlated reservoirs (couplings of site l to arbitrary â † k âm ) or disorder in the reservoir modes would be fairly straightforward to implement in the formalism we are going to discuss and do not change the overall structure.Couplings between system and reservoir, with strength χ, of the type postulated in Eq. ( 1) arise naturally as cross-Kerr interactions between optical or microwave modes.
In this article, our goal is to stay in the classical limit of large amplitudes, replacing b → b.We will now derive an effective set of classical equations of motion upon elimination of the reservoir modes.Let us denote the classical Hamiltonian of the system itself as H, which is a function of the classical mode amplitudes b l and their conjugates b * l .The equation of motion in absence of the reservoir is obtained from the normal ordered Hamiltonian H in general as i ḃl = ℏ −1 ∂H/∂b * l (see Appendix A for more details), with the "Wirtinger derivative" [32] construction, i.e., assuming b l and b * l to be independent variables for the purposes of the derivative.The constant ℏ appears in this classical equation only because we chose to define the classical field amplitudes b l directly as the limit of the quantum operators (and this retains ℏ as a conversion factor between energy and frequency).Using that convenient notation, we can now rewrite the dynamical equation in the form where we denote the temporal derivative using ḃl = db l /dt.Likewise, there is the equation of motion of the driven reservoir mode.Following the standard derivation for a driven, dissipative mode [33], we find with ω l = −∆ − iκ/2 + χ|b l | 2 the effective frequency of mode a l , η the strength of the external drive, and κ the energy decay rate.For convenience we have chosen to work in a frame rotating at the drive frequency ω L of the reservoir modes and introduced the detuning ∆ = ω L −ω c between the drive and the frequency ω c of a reservoir mode .For simplicity we assume ∆, η, and κ to be the same for every lattice site l.
Due to the dynamics of the system amplitudes b l , the reservoir mode frequency ω l will become time-dependent, thus finding the solutions of Eq. ( 3) is a non-trivial problem.We can, however, look for perturbative solutions of the form a l = δa l − iη/ω l , where we chose the second term in analogy to the steady state of the driven dissipative cavity (for χ = 0) but taking into account the time-dependent frequency shift of the cavity due to the interaction with the system.To determine δa l , we first substitute our ansatz for a l in the right-hand side of Eq. (3) , which leads to i ȧl = ω l δa l .Calculating explicitly the time derivative of our ansatz leads to i ȧl = −η ωl /ω 2 l + i (δa l ).Equating the two obtained expressions for i ȧl allows us to find the approximate expression δa l ≈ −η ωl /ω 3 l , where we have neglected the term (δa l ) which is proportional to ω2 l /ω 4 l , ωl /ω 3 l .This is justified within the weak-coupling limit, where In this form the weakcoupling limit can be ensured via a large decay rate or detuning.
Inserting this approximate solution back into the equation of motion for b l and expanding up to second order in χ completes the desired elimination of the reservoir.We note that the expansion in χ introduces no additional approximations but simply ensures consistency with the previously employed weak-coupling limit.The result is the classical effective equation of motion for an arbitrary bosonic system coupled to a reservoir in a particleconserving manner valid in the weak-coupling limit.This is the first main result of the present work.We used Eq. ( 2) to get The second term on the right hand side of Eq. ( 4) is an effective change of the nonlinearity induced by the coupling to the reservoir with strength ( Setting ∆ < 0 (red detuned drive) in Eq. ( 5) leads to an effective attractive interaction.We have also used an appropriate rotating frame by omitting a trivial static shift of the harmonic frequency given by ℏχη 2 /(∆ 2 + κ 2 /4)b l for brevity.The main effect, which will be the focus of our discussion, is the unusual particle-conserving dissipation represented by the last term of Eq. (4).To verify that the latter corresponds to a dissipative term, we consider the change in the system's energy given by with As long as γ > 0, Eq. ( 7) predicts that the energy E S will decrease, confirming the interpretation of the parameter γ as a (dimensionless) decay rate.Dissipation arises in the red-detuned regime, i.e., for ∆ < 0, which yields a positive γ.While ∆ allows one to control the sign of γ, the strength of γ depends on the interplay between the drive strength η, the decay rate κ of the reservoir mode, and the coupling χ to the system.We show in Fig. 2 the dependence on the effective parameters.
Even though we concentrate on classical dynamics, it is insightful to understand the origin of the decay mechanism in the original quantum mechanical picture, where the incoming red-detuned photon is up-scattered to the higher resonance frequency of the reservoir mode by absorbing an energy excitation of the bosonic system, thereby leading to dissipation.At the same time, the total particle number, N = l |b l | 2 is conserved if H S conserves N , which can be readily verified by taking the time derivative of N and employing Eq. ( 4).
Formulas like Eq. ( 7) are known from the theory of laser cooling [34], and its recent implementations , e.g., in cavity optomechanics [35].Indeed there is a direct connection of the particle-conserving coupling studied here to quadratic optomechanics, where one couples the light mode to the square of the mechanical displacement [36][37][38], albeit now with the mechanical oscillator replaced by a many-body system.

III. MICROSCOPIC MODEL: THE ANHARMONIC CHAIN
To derive the dynamical consequences of Eq. ( 4), we will in the following invoke a specific illustrative example, namely a chain of anharmonic bosonic modes, corresponding to the classical limit of the Bose-Hubbard chain.This will already produce very rich behaviour which we are able to analyze both numerically and semianalytically.
The classical system Hamiltonian for this chain is where J is the hopping amplitude between the sites of the chain and α denotes the anharmonicity (negative for attractive interactions).Such a chain can describe various systems: Bosonic cold atoms in a 1D lattice with a Hubbard interaction, coupled superconducting transmon qubits (with α < 0 in that case), as well as coupled optical cavities with a Kerr interaction inside the cavities.As we will explore in some detail below, its continuum limit can give rise to solitonic solutions.Without loss of generality we assume J to be positive in Eq. ( 8) [39] A. Nonlinear Schrödinger Equation with Particle-Conserving Dissipation We now focus on the 1D anharmonic chain described by Eq. ( 8).To obtain an equation describing the effective dynamics of the system, we use Eq. ( 4) with H S = H chain .This yields an effective nonlinear Schrödinger equation on a lattice, which includes particle-conserving dissipation (in contrast to the nonlinear Schrödinger equation with linear loss [40]): where γ is given by Eq. ( 7), g = α + δg is the renormalized coupling strength of the nonlinearity due to the reservoir [see Eq. ( 5)], and n = 1, . . ., L labels the lattice sites.In the following, we consider periodic boundary conditions, i.e., b j = b j+L (j ∈ Z).In our analytic derivation we further consider open boundary conditions where b 0 = b L+1 = 0.An alternative derivation of Eq. ( 9) is presented in Appendix B where we start from the general exact solution.
Within the effective dynamics, the amount of relevant parameters is reduced from formerly six to three, with J setting the overall frequency scale and the dimensionless parameters (g/J, γ) determining the qualitative behaviour (see also Fig. 3 for an illustration on how different microscopic parameters can lead to the same effective dynamics).
For large arrays (L ≫ 1) and assuming that the lattice spacing δx = 1 is smaller than typical wavelengths, i.e., the solution varies slowly on the scale of the lattice spacing, it is possible to consider the continuum limit of Eq. ( 9).To this end, we introduce the continuous function Ψ(x, t) with x ∈ [0, L], which obeys Ψ(x = n, t) = b n (t) for n = 1, . . ., L. To obtain the continuous limit of Eq. ( 9), we replace the finite difference [Ψ(x + δx) − 2Ψ(x) + Ψ(x − δx)]/δx 2 by the second order derivative of Ψ(x) with respect to the position x ∂ 2 Ψ/∂x 2 .This yields the particle-conserving dissipative nonlinear Schrödinger equation (PCDNSE): with the parameters γ [see Eq. ( 7)], g, and J defined above [see Eqs. ( 7) and ( 9 5) and ( 7) as functions of the detuning: The correction of the effective onsite interaction δg (black solid) as well as the dissipation parameter γ (dashed gray).The dissipation parameter is positive in the red detuned domain (∆ < 0) and changes sign in the blue detuned domain (∆ > 0).Both the correction of the onsite interaction, as well as the dissipation parameter scale with the second order of the bath interaction but have different limiting behaviour at large detuning ∆. 10) constitutes the second main result of our work.It will form the basis for our analysis in the remainder of this work, deriving the physical consequences of particle-conserving dissipation for this paradigmatic 1D model.
Since the nonlinear Schrödinger equation without dissipation is known to harbour solitons, we investigate their dynamics in the presence of particle-conserving dissipation.For this purpose, the sign of the dissipation parameter γ, which is positive (negative) for red (blue) detuning, is of the utmost importance [see Eq. ( 7)].We show in the following that in the red-detuned regime, wave packets in the form of soliton solutions can be stabilized and deaccelerated.In contrast to the red-detuned regime, they are accelerated in the blue-detuned regime and any perturbations away from stable solutions can be amplified, leading to wave packets that can break up into several parts.In the following we mainly focus on the dynamics in the red-detuned regime.
To verify the validity of Eq. ( 10) we simulate both the original equations describing the classical dynamics of a discrete lattice of bosonic modes coupled to driven cavities [see Eqs. ( 2), (3) and for the Hamiltonian (8) as well as Appendix B] and the continuum particleconserving dissipative nonlinear Schrödinger equation given by Eq. (10).We compare the resulting dynamics for several different parameter choices.
As an initial condition we choose a stable soliton, that will be introduced in more detail in Sec.III B. We set up the simulation such that the extent of the soliton is some fraction of the total chain length L, and we keep this fraction constant while we vary L in order to discuss the deviations from the continuum limit (which is attained perfectly when L → ∞).
In Fig. 3 we compare the dynamics of the discrete model before eliminating the reservoir modes [see Eqs. ( 2), ( 3) and ( 8) as well as Eqs.(B2) and (B3) for further details] with the PCDNSE.We plot the solutions for different array lengths L [see Fig. 3(a)] and different reservoir parameters [see Fig. 3(b)].Our results show that for long arrays, L ≫ 1, and within the weakcoupling limit both dynamics are in good agreement.Due to the intensity-intensity coupling to the reservoir the site occupation needs to be accounted for in the effective coupling strength such that the weak-coupling limit corresponds to having χb max 2 ≪ |i∆ + κ/2| and χJb max 2 ≪ |i∆ + κ/2| 2 , where we have defined b max = max l |b l |.We note that b max can be related to the occupation of the soliton divided by the width of the soliton.As it can be observed in Fig. 3(b) if the weak coupling conditions are not fulfilled, the discrete and continuous solutions do not agree anymore because integrating out the cavities is no longer justified.
As anticipated, the continuum approximation becomes worse for solitons that (initially) extend only over a few sites, e.g., about 10 sites for L = 400.However, even in this case the continuum equation already agrees reasonably well with the site occupations of its discrete counterpart [see Fig. 3(a)].We stress that within the continuum and weak-coupling limits different parameters that produce the same effective nonlinear coupling strength g and decay rate γ lead to equivalent dynamics.For this reason, we focus in the following on a better understanding of the dynamics of the continuum equation [see Eq. ( 10)].

B. Variational Ansatz
Employing the PCDNSE allows us to use (semi-)analytic methods.In particular we use a variational approach to describe the dynamics of collective coordinates [40][41][42][43], here applied to obtain the dynamics of a single soliton.We consider the ansatz with the collective coordinates given by the amplitude ψ(t), the center position x 0 (t), the velocity v(t), the width w(t), as well as a parameter quantifying the (quadratic) phase dispersion d(t) and the global phase φ(t).A straightforward calculation (see Appendix D) yields the equations of motion for the collective coordi-   (10)] and the discrete equations of motion including the reservoir dynamics, see Eqs. ( 2) and (3).In panel (a) we compare simulations for different array lengths to the PCDNSE and in (b) different cavity parameters leading to the same effective evolution.Only when the continuum limit or weak-coupling limit are not satisfied (small L or ∆ = −2 J, respectively) the discrete dynamics are in disagreement with the dynamics of the PCDNSE.Unless otherwise specified (see legend) we use L = 800, κ = J, ∆ = −0.1J,η = J and χ, α such that we get the effective parameters g = −0.1 J, γ = 0.05.The distributions are taken after an evolution for Jt = 50(L/400) 2 of an initial stable soliton [see inset of (b)] with height ψ(0) = 1 and velocity v(0) ≈ 0.48.The PCDNSE uses a length of L = 400.For the evolution according to the Langevin EOM, we assume that the cavities are initially in their unperturbed steady state η/(κ/2 − i∆).
nates (EOMs): The effects due to reservoir engineering emerge from all the terms proportional to γ.When setting γ = 0, we can verify that these equations coincide with those found in the literature for soliton dynamics [40].
As a first step in our analysis, we determine the nontrivial stable soliton using Eqs.( 12)- (17).To find the properties of the stable soliton, we impose both the width and the amplitude to remain stationary: ẇ(t) = ψ(t) = 0.This condition implies that the dispersion vanishes, i.e., d(0) = ḋ(t) = 4J π 2 w 4 (t) + 2gψ 2 (t) π 2 w 2 (t) = 0.This immediately leads to a relation between width and amplitude of the stable soliton given by and imposes the quadratic phase dispersion to be zero, d(t) = 0.The initial conditions of the remaining parameters v(t = 0), x 0 (t = 0), φ(t = 0) can be chosen independently.
It can be shown that both the width and amplitude of the stable soliton are time-independent.To this end, we use that the particle number, is conserved.Substituting Eq. ( 18) into Eq.( 19) and solving for w leads to the width of the stable soliton: The time-independent amplitude can be obtained by substituting Eq. ( 20) into Eq.( 19) and solving for ψ, which yields We can also derive the shape of the stable soliton using the ansatz postulated in Eq. (11); this would consist in minimizing the energy of the unperturbed nonlinear Schrödinger equation for a fixed particle number In the attractive regime, i.e., g < 0, the last term of Eq. ( 22) related to the onsite interaction is negative and, thus, competes with the remaining positive terms allowing for a stable soliton.Furthermore, Eq. ( 22) reveals the dependence of the energy of the stable soliton on the collective variables, which demonstrates that it is possible to change the energy of the stable soliton without breaking it apart by changing its shape (amplitude ψ and width w) or reducing the dispersion d and or velocity v.
As we show below, this can be achieved by leveraging the coupling to the engineered reservoir.By substituting w = N/2ψ 2 [Eq.(19)] into the energy Eq. ( 22) we find that the kinetic energy term of the group velocity v is independent of the soliton shape (determined by ψ).Minimizing Eq. ( 22) with respect to ψ and d for attractive onsite interaction (g < 0) at fixed N and v leads to d = 0, and Eqs. ( 20) and ( 21).The particular shape of the soliton [Eqs.(20) and ( 21)] arises due to the competition between part of the hopping term and the onsite interaction, in particular the third and last terms of Eq. ( 22), respectively.The energy of the stable soliton is given by The first term of Eq. ( 23) is the kinetic energy which depends quadratically on the velocity v.The second term corresponds to a potential energy that depends on the particle number N as well as the system parameters g and J.The energy of the stable soliton [see Eq. ( 23)] is linked to the evolution of the global phase via dE SS /dN = v ẋ0 − φ [see Eqs. ( 24) and ( 26)].

Dissipative dynamics of the stable soliton
With these considerations in place we begin by studying the dynamics of a stable soliton.Solitons in Bose-Hubbard-type models have been discussed in the quantum domain [44,45] without the particle-conserving dissipation we consider here.Because the PCDNSE conserves the particle number but dissipates the energy, the dynamics need to change the velocity to influence the energy as it can be seen from Eq. ( 23).This can be simply understood using the EOMs [see Eqs. ( 12)-( 17)] for a stable soliton with the (dimensionless) velocity damping rate given by Note that the velocity damping rate of the stable soliton depends on the ratio g/J because this fraction determines the particular shape of the stable soliton, i.e., the amplitude ψ and width w [see Eq. ( 13)].
Since we need an attractive interaction (g < 0) to get a stable soliton, we conclude that the velocity damping rate Γ is positive when the original effective damping parameter γ is also positive, i.e., when ∆ < 0 (red-detuned regime).Therefore, the main insight these equations provide is that the velocity v of the soliton is damped at a rate of ΓJ for ∆ < 0. This highlights that with the particle-conserving coupling to the engineered reservoir it is indeed possible to de-accelerate and eventually even freeze the motion of a stable soliton while conserving its shape.On the other hand when ∆ > 0, one can observe instead exponential acceleration.
We show this behavior in Fig. 4 where we compare the relative (de-)acceleration of the soliton v/vJ as function of the damping rate Γ as predicted by the variational ansatz to a numeric simulation of the PCDNSE.We simulate numerically the evolution up to Jt = 4 and then extract the soliton coordinates using a least squares fit to the ansatz given in Eq. (11).We then estimate v/J using the finite difference [v(Jt) − v(0)]/Jt.We stress that while in the fully displayed red-detuned domain the numerical results agree nicely with the predicted linear dependence on the damping rate Γ, in the blue-detuned regime there is a sudden increase of the deviation from the predicted linear dependence (leftmost purple square in Fig. 4).We attribute this deviation to the amplification of a small departure from the stable soliton, which is possible in the blue-detuned domain [see Eq. ( 16)].
In the simulation numerical errors can lead to a small deviation from the stable soliton which can result in a non-vanishing dispersion d ̸ = 0 that then gets amplified for γ < 0 [see Eq. ( 16)].We verified this by reducing the absolute and relative tolerances of the simulation and thereby the potential deviation from the stable soliton, resulting in final numerical values for the anti-damping rate in agreement with the variational calculation (see purple star in Fig. 4).Due to the amplification, even small deviations from the stable soliton can eventually lead to a break-up of the soliton.As (external) perturbation can also occur in the real world, the blue-detuned regime needs to be treated with great care if the goal is to prepare stable solitons.

Shape stabilization
In contrast to the shape breakdown in the bluedetuned regime, which we discussed above, the reddetuned domain supports the stabilization of solitons when they start off with a small deviation from the stable shape.We now analyze this behaviour in some more detail.Since the EOMs conserve the particle number, the initial deviation from a stable soliton can be quantified by the initial phase dispersion d(0) and initial relative difference in amplitude δ = ψ(0)/ψ SS − 1.For simplicity we focus on the second type of deviations here, i.e., deviations in the initial amplitude.
To study the phenomenon of soliton stabilization, we need to consider the long-time limit of the evolution.While this can be done efficiently using the variational approach, we must make sure that we are not in a regime where small amplitude deviations from the stable soliton characterized by δ result in the wave packet breaking up into multiple components.To identify the regime of validity, we employ a hybrid approach: We first calculate a short-time evolution using the computationally more expensive PCDNSE to ensure the stability of the soliton.If the soliton retains the correct shape, we then use the  .Velocity damping rate as a function of gγψ 4 (0).We see good agreement between the variational approach and numerical results of the particle conserving dissipative nonlinear Schrödinger equation (PCDNSE).The deviation at the smallest value is due to the soliton breaking up in the simulation, underlined by the datapoint given by the star where the relative and absolut tolerances where decreased to 10 −13 and 10 −12 respectively from 10 −8 for both.The inset shows that the deviation |δ(Jt = 4)| = |ψ 2 (Jt = 4)/ψ 2 (0) − 1| (pentagons, diamonds) as well as d(Jt = 4) (triangles) from the stable soliton stays small for the evolution.The PCDNSE simulation is for a stable soliton with ψ(0) = 1, x(0) = L/8, v(0) ≈ 0.49, and L = 600 and for a duration Jt = 4.
collective coordinates to estimate the long time evolution.
We illustrate and justify this approach in Fig. 5 where we compare the PCDNSE and variational approach.Based on the short-time evolution of ψ(t) in Fig. 5(a), one can identify that the soliton breaks up into multiple wave packets for δ = 0.3 within the PCDNSE.Thus, we conclude that the variational approach is not suited in this case.
For the remaining values of δ = −0.1,0.01 we additionally compare the dynamics on a short time scale (after a previous long evolution) in panel (b) and verify that the PCDNSE solitons still agree with the spatial form of the general soliton ansatz (11) in panel (c).Both underline the validity of the hybrid approach.
In Fig. 5(d) we plot the maximal deviation of the soliton amplitude from the steady state value within a timespan of duration 5×10 4 /J to confirm that for red-detuned dynamics the solitons converge to a state where the velocity v and dispersion d tend to zero, while the amplitude ψ(t) and width w(t) tend to the values associated to a stable soliton, i.e., ψ SS and w SS [see Eqs.(20) and (21)].
Even for the smaller |δ| values, the PCDNSE simulation and the collective coordinates approach also deviate due to higher order dispersion terms that are not accounted for in our collective coordinate ansatz.Because the higher order dispersion terms do not influence the (modulus square) shape they can be incorporated in the variational approach to improve the collective coordinates to higher order in future work.The simulations in Fig. 5 furthermore suggest that the ansatz we employed shows the the continued oscillation around the steady state amplitude after a prolonged evolution time of Jt = 2 × 10 6 and (c) displays the modulo square distribution after this evolution.To quantify the long time dynamics we compare the envelope of the deviation of the maximal amplitude ψ(t) deviation from the steady state amplitude ψSS on a long timescale in these cases (d).To this end, we calculate the envelope as the maximal deviation within a time-window of length 5 × 10 4 .The relaxation process is sketched in the inset of (d).We use g = −0.1J,N = 1, γ = 0.1, L = 10wSS, and only deviate from the stable soliton via δ ̸ = 0 (i.e.d = 0).
still gives an estimate for the timescale of the relaxation as well as the qualitative dynamics of the relaxation.The relevance of the hybrid approach becomes evident in the blue-detuned regime even in the case where the initial soliton is stable; in this regime a slight perturbation from the stable soliton solution can be amplified and lead to a break up of the wavefunction as we discussed in the previous section.
Additionally, the blue-detuned regime reveals the main drawback of the variational approach employed here, i.e., the complete breakdown of the solution as soon as the space spanned by the ansatz is insufficient to describe the dynamics.However, we want to stress again the great simplification that is possible as long as the solutions stays within this solution space, i.e. replacing the non-linear partial differential equation (10) with a set of ordinary differential equations ( 12)-( 17).This enabled our analytic results like the derivation of the velocity damping rate for a stable soliton [see Eqs. ( 24)-( 26) and Fig. 4].
We emphasize that the breakdown of the collectivecoordinates approach is a problem inherent to taking the highly restricted ansatz given by Eq. (D1), and is not due to the particular application to the model of this article; this is particularly visible as the approach works better for bigger γ (for fixed g/J).
Combining the insights gained above for the reddetuned domain, regarding the velocity damping and shape stabilization, as well as the separate finding that large deviations lead to a break up of the wavepacket, we suspect that the eventual steady states for the break up scenario can be described by several spatially separated stable solitons.

C. Two Solitons
After investigating single solitons in great detail in the previous section we want to use these insights now to investigate the collision of two solitons.To this end, we investigate the effect of the interaction by comparing the evolution of a single traveling soliton using the solution obtained using the collective coordinates to the evolution of two solitons propagating in opposite directions, evolved (numerically) using the PCDNSE [using Eq. ( 22)].In particular, we calculate the total energy numerically and compare it to the energy of a single soliton, calculated analytically according to the collective coordinates ansatz [see Eqs.(23) and (25)].
We display this in Fig. 6 where we see that initially, i.e., before the interference begins (see inset), the prediction agrees with the expectation of two separate solitons, i.e., the energy of the two solitons is twice the energy of the single soliton.During the evolution, when the solitons collide, the (unperturbed) energy of the two solitons can even increase temporarily.We ascribe this behaviour to interaction terms that go beyond the unperturbed single soliton [see Eq. (22)].However, it is apparent that during the soliton collision energy is being dissipated beyond what would be expected for separated solitons.We attribute this to the increase of the variation of the dissipation term ∝ Im Ψ * ∂Ψ ∂x 2 Ψ as the solitons interfere with each other.
As a consequence, in a system of many solitons, energy is predominantly dissipated whenever two solitons collide.Beyond that, there is the usual velocity-damping for each individual soliton, already discussed above.

IV. CONCLUSION
In this article we have analyzed the behaviour of a bosonic many-body system subject to particle-conserving dissipation, treated in the classical limit.Particleconserving dissipation is important in modern bosonic quantum simulators for platforms such as photonic systems that can be employed to simulate the interesting physical scenario of particle-conserving equilibration dynamics when subject to an engineered reservoir of this manner.The results for the classical limit obtained here will provide a guide to understanding such experiments in the limit of large amplitudes.We have derived a general equation of motion describing the effective dissipation.We have then applied it to the special case of a Bose-Hubbard chain, where we have obtained a particular new version of a dissipative nonlinear Schrödinger equation.
We used numerics and a collective coordinate approach to study the behaviour of solitons in such a scenario and showed how reservoir engineering can relax solitons to a stable shape as well as accelerate and de-accelarate moving solitons.Lastly we briefly investigated the interaction of two solitons, where we found that extra energy is dissipated during the collision of solitons.
Further works could build upon these insights to study the breakup of wave packets into soliton trains and the interactions of multiple solitons as well as their long-term behaviour, all in the presence of particle-conserving dissipation.The behaviour of a 2D bosonic field coupled to this kind of dissipation remains completely unexplored territory as well.
Following the standard derivation for driven dissipative cavities (see [33] for a review) and considering the regime where the quantum fluctuations are negligible compared to the classical amplitudes we find the equations of motion of the classical amplitudes A n and B n : with the detuning ∆ = ω L − ω c between the cavity and the drive.The terms proportional to the drive rate η and the decay rate κ arise due to integrating out the environment following the standard procedure for a driven dissipative cavity (with an interaction between the environmental modes ê(ω) with frequency ω and the cavity of the form ê(ω)â † + h.c..We note that in our notation all constants apart from ∆, α are positive.First we note that only the hopping term leads to a change in the site occupation For a stable soliton we can relate B max (t) = ψ(t) = N/2w SS , i.e. the square root of a density given by the number of excitations in the soliton devided by the width of the soliton.
Taking the ansatz such that we can integrate both sides of the equation to find Where we approximate the integral t 0 e g(τ ) dτ with g To get the correct order of these terms we consider  where we used that we can add a real number inside the imaginary part to solve the integral, and that Ψ needs to either vanish at the boundaries (for open boundary conditions) or Ψ(0) = Ψ(L) for periodic boundary condi-tions.Therefore we demonstrated that the particle number N = dx |Ψ| 2 (x, t) is constant.

Appendix D: Variational Ansatz
In this appendix we provide the necessary details to understand the derivation of the equations of motion for the collective coordinates.Without any coupling to the environment (γ = χ = 0) the particle conserving dissipative nonlinear Schrödinger Eq. (10) becomes the standard unperturbed nonlinear Schrödinger equation which supports soliton solutions.With this motivation we study the evolution of a generalized soltion ansatz 2 d(t) + iφ(t) sech x − x 0 (t) w(t) (D1) within a variational approximation where the timedependent parameters are called collective coordinates [41,42].This approach was also employed to study soliton solutions in different generalized non-linear Schrödinger equations in recent works [40,43].Within this approach we first calculate the conservative Lagrangian for the unperturbed nonlinear (the non-conservative term of the particle conserving dissipative nonlinear Schrödinger Eq.).Evaluating the integrals of the right hand side combined with some algebra leads to the equations of motion of the main text ( 12)-( 17).

2 𝜂 2 /𝜅 4 ]Figure 2 .
Figure 2. Effective parameters as defined in Eqs.(5) and(7) as functions of the detuning: The correction of the effective onsite interaction δg (black solid) as well as the dissipation parameter γ (dashed gray).The dissipation parameter is positive in the red detuned domain (∆ < 0) and changes sign in the blue detuned domain (∆ > 0).Both the correction of the onsite interaction, as well as the dissipation parameter scale with the second order of the bath interaction but have different limiting behaviour at large detuning ∆.

Figure 3 .
Figure 3.Comparison of the PCDNSE [Eq.(10)] and the discrete equations of motion including the reservoir dynamics, see Eqs. (2) and (3).In panel (a) we compare simulations for different array lengths to the PCDNSE and in (b) different cavity parameters leading to the same effective evolution.Only when the continuum limit or weak-coupling limit are not satisfied (small L or ∆ = −2 J, respectively) the discrete dynamics are in disagreement with the dynamics of the PCDNSE.Unless otherwise specified (see legend) we use L = 800, κ = J, ∆ = −0.1J,η = J and χ, α such that we get the effective parameters g = −0.1 J, γ = 0.05.The distributions are taken after an evolution for Jt = 50(L/400) 2 of an initial stable soliton [see inset of (b)] with height ψ(0) = 1 and velocity v(0) ≈ 0.48.The PCDNSE uses a length of L = 400.For the evolution according to the Langevin EOM, we assume that the cavities are initially in their unperturbed steady state η/(κ/2 − i∆).

Figure 4
Figure 4. Velocity damping rate as a function of gγψ 4 (0).We see good agreement between the variational approach and numerical results of the particle conserving dissipative nonlinear Schrödinger equation (PCDNSE).The deviation at the smallest value is due to the soliton breaking up in the simulation, underlined by the datapoint given by the star where the relative and absolut tolerances where decreased to 10 −13 and 10 −12 respectively from 10 −8 for both.The inset shows that the deviation |δ(Jt = 4)| = |ψ 2 (Jt = 4)/ψ 2 (0) − 1| (pentagons, diamonds) as well as d(Jt = 4) (triangles) from the stable soliton stays small for the evolution.The PCDNSE simulation is for a stable soliton with ψ(0) = 1, x(0) = L/8, v(0) ≈ 0.49, and L = 600 and for a duration Jt = 4.

Figure 5 .
Figure 5. Shape stabilization of ill prepared solitons.The maximal amplitude of the PCDNSE (dashed lines) is compared to the prediction of the collective coordinates (solid lines) in (a),(b) for the different deviations δ = ψ(0)/ψSS − 1 [see legend in (c)].While for δ = −0.1,0.01 both solutions oscillate around the steady state amplitude, for δ = 0.3 the soliton shape breaks down and this oscillation is visible in neither solution.Therefore, we investigate the dynamics for long times for the smaller δ = −0.1,0.01 (b)-(d).Panel (b)shows the the continued oscillation around the steady state amplitude after a prolonged evolution time of Jt = 2 × 10 6 and (c) displays the modulo square distribution after this evolution.To quantify the long time dynamics we compare the envelope of the deviation of the maximal amplitude ψ(t) deviation from the steady state amplitude ψSS on a long timescale in these cases (d).To this end, we calculate the envelope as the maximal deviation within a time-window of length 5 × 10 4 .The relaxation process is sketched in the inset of (d).We use g = −0.1J,N = 1, γ = 0.1, L = 10wSS, and only deviate from the stable soliton via δ ̸ = 0 (i.e.d = 0).

Figure 6 .
Figure 6.Influence of the dissipative dynamics due to the interaction of two stable solitons.The figure depicts the energy of the two solitons compared to the single soliton solution as a function of time during the interaction for different γ (see legend).The inset shows |Ψ(x, t)| 2 for γ = 10 −2 where the values range from 0 (dark) to ≈ 3.62 (light).We use g = −0.1Jand the initial solitons have ψ(0) = 1 (N ≈ 9), and |v(0)| ≈ 0.5 with opposite signs and are placed 5w(0) apart on a space of size 20w(0).The figure shows that the interaction of the two solitons can enhance the energy dissipation.