Reservoir engineering with localized dissipation: dynamics and pre-thermalization

Reservoir engineering lattice states using only localized engineered dissipation is extremely attractive from a resource point of view, but can suffer from long relaxation times. Here, we study the relaxation dynamics of bosonic lattice systems locally coupled to a single squeezed reservoir. Such systems can relax into a highly non-trivial pure states with long-range entanglement. In the limit of large system size, analytic expressions for the dissipation spectrum can be found by making an analogy to scattering from a localized impurity. This allows us to study the cross-over from perturbative relaxation to a slow, quantum-Zeno regime. We also find the possibility of regimes of accelerated relaxation due to a surprising impedance matching phenomena. We also study intermediate time behaviors, identifying a long-lived"prethermalized"state associated that exists within a light cone like area. This intermediate state can be quasi-stationary, and can very different entanglement properties from the ultimate dissipative steady state.


I. INTRODUCTION
In quantum information applications, states with entanglement and other non-classical properties serve as basic resources. A powerful approach to generating such states is to couple the system of interest to a tailored dissipative environment, so that the resulting system steady state is the desired non-trivial quantum state; this approach is known as reservoir engineering [2,3]. There is by now a growing body on work on utilizing such techniques, ranging from the stabilization of systems with a few degrees of freedom (see e.g. Refs. [4][5][6][7][8][9], to the preparation of many body states via controlled system-wide dissipation (see e.g. Refs. [10][11][12][13][14][15][16][17]. Work has also showed that in some cases many-body states can be stabilized by controlling dissipation in a limited, localized spatial region only [18][19][20][21]. Such an approach was recently implemented in a superconducting circuit experiment to stabilize a bosonic Mott insulator [22]. Our recent work in Ref. [1] analyzed a particularly striking example of this local approach to reservoir engineering: it demonstrated that an entire class of free boson lattice systems with a generalized chiral symmetry can be stabilized in this manner by making use of the lattice symmetry and a single, localized, squeezed dissipative reservoir. This allows the stabilization of non-classical, often highly non-locally entangled Gaussian pure states [1]. While the ability to stabilize non-trivial states using a localized coupling to engineered dissipation is extremely attractive in terms of the needed experimental resources, an obvious potential drawback is that the timescale for relaxing to the steady state can be extremely long (making one more susceptible to unwanted dissipation and decoherence). The relevant relaxation time in such schemes typically scales with system size. Intuitively, this makes sense: correlations are generated locally at the dissipative site, and take time to propagate throughout the lattice. As a result, a long intermediate time regime exists before the final steady state is reached. Understanding this dynamical phenomena is of interest for many reasons. It would allow one to optimize the relaxation process and speed when local dissipation is used in reservoir engineering. In addition, the system's path from its initial configuration to the dissipative steady state may exhibit different qualities and new physics that is present in neither.
In this article we study the relaxation dynamics in this kind of reservoir engineering setup. We focus on the simple but paradigmatic case considered in Ref. [1]: a bosonic lattice system described by a quadratic hopping Hamiltonian, coupled to a Markovian reservoir on just a single site. As we show, the relaxation here can be fully characterized by the eigenvalues of the system's dynamical matrix. Further, the corresponding relaxation rates have a strong parallel to the physics of scattering off a localized potential and the Friedel sum rule. Analyzing the behavior of these modes as local dissipation strength is increased, we see a transition from perturbative relaxation of the system's original eigenstates to a quantum Zeno-like decoupling of the dissipative site. While these limiting cases could easily be anticipated, we also find that under some circumstances, there can be a surprising resonant enhancement of mode relaxation rates for intermediate coupling strengths; this can be interpreted as an impedance matching phenomenon. We also consider the overall evolution of the system from its initial state, describing the intermediate-time behavior in large systems. We find the existence of a quasi-steady state, whose form and duration is dictated by ballistic propagation physics. This quasi-steady state describes sites where correlations have had time to propagate from the dissipative site, but have not had time to reach the system boundary and then return. Surprisingly, we find that this quasi-steady state can have a very different form and pattern of entanglement than the final steady state of the system.
The remainder of the article is organized as follows. In Section II we outline the basic model of the class of dissipative systems we consider. In Section III, we find the eigenmodes of our system's dynamical matrix and corresponding "dissipation spectrum", and analyze its behavior. We provide several numerical examples and one analytical solution of a sample system. Finally, in Section IV, we turn to the intermediate-time behavior of the system, where the state differs both from the initial state and the final stabilized state. We show the emergence of a quasi-steady state within the light cone corresponding to ballistic propagation of correlations from the dissipative site.

II. MODEL
We consider the same class of model studied in our previous work [1], though to start, do not impose any sort of symmetries. We consider a generic particle conserving quadratic bosonic Hamiltonian, whereâ n (â † n ) is the annihilation (creation) operator for a boson on site n, and the Hamiltonian matrix H consists onsite potentials H n,n = V n and hopping elements H m,n = J m,n . The site labels n describe an arbitrary d-dimensional lattice with N sites, and we do not assume any symmetry or translational invariance to begin with. This Hamiltonian describes a range of bosonic systems, including coupled arrays of superconducting cavities or mechanical oscillators. We linearly couple a single "drain" site, marked by n 0 to a Markovian, Gaussian reservoir at zero temperatrure. Using standard input-output theory [23], the Heisenberg-Langevin operator equations of motion are theṅ where Γ is a dissipation rate parameterizing the strength of the coupling to the reservoir andζ(t) is a Gaussian white noise operator describing its vacuum fluctuations. We will eventually takeζ(t) to describe squeezed vacuum fluctuations, but for now keep things general (as the system's characteristic relaxation rates are independent of the nature of the noise). The Hamiltonian of Eq. (1) can also be written in diagonal form using its energy eigenmodesb i and corresponding eigenstate wavefunctions ψ i [n]:Ĥ Without loss of generality, we label these energy eigenmodes so that ε i+1 ≥ ε i . We will also take the spectrum to be non-degenerate with ψ i [n 0 ] = 0 for all modes; in other words, we focus on the portion of the spectrum that is coupled to the drain. Any other modes are unaffected by the dissipation dynamics (see discussion in [1]).
In the eigenmode basis, including the coupling to the reservoir, the operator equations of motion take the forṁ where the dynamical matrix A is given by Here are the phase and magnitude of the coupling between the energy eigenmode i and the reservoir.

III. DYNAMICAL MODES AND CHARACTERISTIC DISSIPATION RATES
For a generic Markovian system, it is common to characterize relaxation time scales by considering the the eigenvalues of the Liouvillian which governs the evolution of the system's reduced density matrix [24]. Here, the linearity of our system makes life much easier. We can fully characterize the system's dynamics and relaxations by simply diagonalizing the dynamical matrix A defined in Eq. (5). Its eigenvalues are the characteristic mode frequencies of the linear dynamics of Eq. (4). The non-Hermitian nature of A means these mode frequencies will be complex, with non-positive imaginary part corresponding to loss. We therefore calculate the left-eigenvectors V ij and eigenvalues λ i of A, As shown in Appendix A, the left-eigenvectors are given by and its eigenvalues λ i are the solutions of the self-consistency equation These eigenvectors define a set of operators whose time evolution can be immediately integrated, Note that these operators are not a set of independent, canonical annihilation operators, but instead satisfy Despite this, they will be useful in the analysis which follows.

A. Dissipation Spectrum
Our goal is to use Eq. (9) to calculate the dynamical matrix eigenvalues λ i in the limit of a large system. Note that the evolution described in Eq. (2) is analogous to scattering off a localized impurity, except now, the impurity potential is imaginary (i.e. effectively non-Hermitian). Thus, similar to the standard treatment of potential scattering [25], we look for dynamical eigenvalues λ i that are just a small shift of the original energy eigenvalues ε i , The real and imaginary parts of this shift corresponds to an energy shift δν i , and an inverse lifetime γ i /2, The sum in Eq. (9) can then be rewritten In the case of a large system, we generally expect the shift and the relaxation rate to be quite small, |δλ i | ∼ 1/N . The sum in Eq. (14) is then dominated by the resonant terms, |ε i+n − ε i | ∼ |δλ i |. Further, we generically expect the unperturbed density of states and the coupling to the drain,Γ i+n , to remain roughly constant for these resonant contributions. We thus approximate the spectrum and the local coupling for modes near i as where ∆ i = (ε i+1 − ε i−1 )/2 is the energy spacing near ε i . We can then approximate Eq. (14) as this resonant portion, This is equivalent to assuming a constant density of states, and eigenmode wavefunctions that have a constant amplitude at the drain site. The sum in Eq. (16) can be immediately evaluated to find This is the usual expression for an s-wave scattering phase shift from a localized impurity potential [26][27], with the potential taken to be imaginary, V → iΓ/2. The phase shift determines the energy shift of the modes, and so the resulting imaginary portion gives the relaxation rate, Eq. 18 is a key result of this paper. It demonstrates that the relaxation rate associated with mode i is controlled both by the local level spacing ∆ i of the unperturbed spectrum near E = ε i , as well as the coupling rate Γ i between mode i and the localized dissipative reservoir. ∆ i will generically scale as 1/N ; while the same is true forΓ i if the mode i is extended. In this case the relaxation rate γ i will also scale as 1/N . For a localized eigenmode i, we instead expectΓ i ∼ exp[−r i /ξ] for the mode's distance from the drain r i and its extent ξ.
We observe that the relaxation rates of the system are determined by the parameterΓ i /∆ i . This can be understood as the ratio of the mode's dwell time (∆ −1 i ) to its dissipative lifetime (Γ −1 i ). We see two limits for this rate, In the weak coupling limit,Γ i ∆ i , the original modes are perturbatively coupled to the dissipative bath, yielding γ i Γ i , a result that would be expected from a standard Fermi's Golden Rule calculation. In this limit the relaxation time of any localized mode grows exponentially with its distance from the coupled site. In contrast, for extended modes the dissipation rate will scale as 1/N . The strong dissipation limit,Γ i ∆ i , can be understood in terms of quantum Zeno physics [28]: the coupled site is measured by its bath faster than it can interact with these modes, and their dissipation is suppressed.
It is also interesting to consider the behaviour of γ i in the regime where there is an approximate matching of timescales,Γ i ∼ ∆ i . In this case, there is a form of impedance matching, as the propagation rate of waves arriving at the coupled site matches that of waves radiated into the bath. To understand the relaxation rates at this regime, we must take into account the non-resonant portions of the sum in Eq. (16). We define the remainder This represents the correction to approximating the system as having constant wavefunctions and level spacing.
Recall that we expect δλ is extremely small in the limit of large system size: δλ ∼ ∆ i ∼ 1/N . If the approximations of Eq. (15) hold, the resonant portions of the sums in the square brackets will cancel out, implying that the quantity inside the brackets is non-singular as δλ → 0. The first term above then will be proportional to δλ. It follows that in this limit, R i will be dominated by the last term on the second line of Eq. (20). Hence, to leading order in 1/N : It follows that the only dependence on the strength of the dissipation Γ in Eq. (21) is through an overall factor. We thus define the rescaled remainder, The dimensionless parameter r i is independent of the strength of the dissipation, and can be calculated for a given system and drain site. It characterizes how strong the corrections to the approximation in Eq. (15) are for any particular mode i. Note that while the resonant termS i , as defined in Eq. (16), is determined purely by the effective dissipation and local level spacing for each energy eigenmode, the remainder r i depends on the global form of the entire spectrum, as well as fluctuations from mode to mode of the wavefunction at the drain-site. This means that the particular values of r i depend on the the specifics of the Hamiltonian involved. However, as j |ψ i [n 0 ]| 2 = 1, we can interpret Eq. (22) as a weighted average of the terms in the summation. This allows it to capture the corrections from the non-uniformity of level spacing and drain-site wavefunctions over the entire spectrum. Recalling that ∆ i is the energy spacing near i, we have 1 ≤ |ε j − ε i |/∆ i , we then have by definition |r i | ≤ 1/π. Furthermore, as the denominator of most terms in the sum is of the order of magnitude of the spectrum, we might expect |r i | = O(1/N α ) for some positive α.
Combining Eqs. (14), (16) and (20) to (22), we find the full expression for the relaxation rates of the eigenmodes, Note that if one knows the original system's eigenstate energies and wavefunctions at site n 0 , then Eq. (23) allows a calculation of the dissipation spectrum for an arbitrary dissipation strength Γ. It thus allows calculation of the dissipation spectrum in a simple manner, without having to rediagonalize the dynamical matrix for each different In the second row, we plot the ratio γi/∆i for three different values of the dissipative coupling Γ (see legend). The energy spacing, ∆i = (εi+1 − εi−1)/2, which is proportional to the group velocity, is plotted in Fig. 2 for these lattices. The circles are calculated by diagonalizing the dynamical matrix (see Section II). The crosses are calculated from Eqs. (22) and (23) using properties of the non-dissipative (Γ = 0) system only. We see remarkable agreement. We observe that while increasing the coupling from Γ = 0.5J to Γ = 2J increases the relaxation rate of all modes, in these cases a further increase to Γ = 4J creates Zeno-like suppression for most modes as the εi = 0 gains a macroscopic relaxation rate.
choice of Γ. Further, for some systems it is possible to analytically calculate the spectrum and r i ; we provide such an example in the next section.
To test its validity, we have directly compared the dissipation spectrum calculated from the approximate expression in Eq. (23) against an explicit numerical diaonalization of the dynamical matrix for a range of 1D and 2D models. Representative results are shown in Fig. 1. One finds a very good agreement with the approximate dissipation rates. We find that γ i retains the qualitative behavior discussed above, including the perturbative regime at small Γ and Zeno behavior at large Γ. Equation (23), however, gives us a full picture of the relaxation spectrum's behavior, as well as the transition between the two regimes. Using this expression, we find: We see that for modes with |r i | |ψ i [n 0 ]| 2 , there is a logarithmic enhancement of in the relaxation rate atΓ i 2 = ∆i π . It is also interesting to consider the ratio γ i /Γ i as a function of Γ; this ratio measures how different a mode's dissipation rate is from the simple, Fermi's Golden Rule estimate. For |r i | ≥ |ψ i [n 0 ]| 2 / √ 3, the ratio is a monotonically decreasing function of Γ. In this case, the mode essentially transitions directly from the perturbative regime to the Zeno regime. For |r i | < |ψ i [n 0 ]| 2 / √ 3, however, γ i /Γ i initially increases, signifying the kind of impedance matching discussed above. Examples of this behavior are shown in Fig. 2 for several systems. We see the evolution of the spectrum take on quite different forms for different systems, with resonant enhancement appearing for none of the modes, for several modes at the same coupling strength, or for different modes at different values of the coupling. We now apply the results of the previous section to a specific model of a one dimensional ring with nearest neighbour hopping (Fig. 3), pierced by a non-zero flux. In this case, we can use Eq. (23) to analytically calculate the dissipation spectrum to leading order in system size.
The system Hamiltonian in this case is given by: We further assume periodic boundary conditions,â N +1 ≡â 1 , and restrict the ring flux ϕ to the interval [0, π]. The system is diagonal in a plane-wave (i.e. momentum) basis,Ĥ where k = 2π N × 0 . . . N − 1. To order the modes so that ε i ≤ ε i+1 , we relabel Labelled this way, every two consecutive eigenmodes have opposite wave-numbers, and in the absence of a flux (i.e. ϕ = 0), are degenerate. These modes have energieŝ To simplify the calculation we choose ϕ = π 2 , where We note that due to translational invariance, b j ,â † Translational invariance guarantees that all modes couple equally to the drain site.
We can now calculate the remainder of Eq. (22), taking N 1, We find, in this case, r i = O 1 N . Note that this is a feature of the ring system, and not generically true, as we have seen for different systems in Fig. 2.
We next calculate the dissipation spectrum to first order in 1/N . For Γ < 4J, we find a critical momentum, where the group velocity matches the dissipation rate, For the relaxation rate, we see an impedance-matching like phenomenon, as discussed above, at this point, In the overdamped case, Γ > 4J, there is no impedance matching. Instead, we see the relaxation rates dropping off away from the middle of the spectrum, In this regime, we also find that Eq. (9) is satisfied by a central rate with a macroscopic relaxation rate, At large Γ, this mode becomes localized to the drain site, and effectively detaches from the rest of the ring. These results are plotted in Fig. 4, along with a numerical calculation for finite systems of several sizes.

IV. INTERMEDIATE TIME BEHAVIOR
Having developed a full understanding of the dissipation spectrum associated with our local reservoir engineering setup, we now turn examining the more global features of the system's evolution from an initially prepared state to the final, dissipation-induced steady state.
As described in Section III, the eigenvalues λ i of the dynamical matrix that characterize the system's time evolution can be usefully expressed as a dissipation-free energy plus a complex shift, λ i = ε i + δλ i ; the imaginary part of this shift encodes the relaxation rate associated with a particular mode. Our analysis revealed that these shifts, including the relaxation rate, are generally small and inversely proportional to system size λ i ∼ 1/N . These basic features imply that the system's relaxation to the steady state can be broken into two parts.
• First, for a relatively long period initial period (whose duration τ ∝ N ), we can to a good approximation ignore the dissipation-induced contributions to the dynamical matrix eigenvalues, and simply replace them by the corresponding energy eigenvalue: where J ∝ i ∆ i is the energy scale for the system's dynamics, i.e. the hopping rate in the systems we consider. During this initial period, the system evolution is well described by the non-dissipative dynamics generated by the system's coherent Hamiltonian. At a heuristic level, particles can be injected into the system at the drain site, but will then propagate ballistically (i.e. according toĤ).
• For longer times, the dissipative contribution to mode eigenvalues is non-neglible, and we have exponential decay associated with relaxation to the steady state: At a heuristic level, this corresponds to a timescale long enough that particles injected from the drain site have had enough time to traverse the system and reflect off its boundaries. This process continues for a long time until the steady state is reached.
We now attempt to see more explicitly how the above picture manifests itself in the state of the system. As the form of the dissipation-induced lattice state will depend on the form of the bath noise, we specify now to a Gaussian bath with a squeezed form. We take the input operator of Eq. (2) to have correlators To understand the intermediate dynamics, we return to Eq. (4) for the evolution of an energy eigenmode, rewriting it asḃ During the initial evolution period, the simplest approximation would be to neglect the damping effect of the bath (last term) and only keep the driving term. This would then correspond to a picture where the bath simply drives the system with correlated pairs of particles, but does not modify its dynamics or response properties.
To get a slightly more accurate approximation of this early-period evolution, we can instead exactly solve the above equation forb i (t) in terms ofâ n0 : The last two terms describe the driving and damping of the drain site by the bath as it is modified by the response of the rest of the lattice (e.g. fluctuations may enter, bounce around the lattice several times, and then finally emerge at n 0 ). This response is non-Markovian at short time scales, when these dynamics are sensitive to the finite bandwidth of the lattice; and at long times, when the discrete, non-uniform density of states in the lattice is significant. However, during the intermediate regime of Eq. (36), a Markovian approximation is sufficient to qualitatively capture this effect. We thus take where J ∝ J is some factor dependent on the effective density of states at the drain.
Combining Eqs. (39) to (41), we arrive at an effective equation of motion, This equation shows simple linear dynamics with a source term. Each mode is coupled to the source via its effective dissipation, with an overall suppression as the drain siteâ n0 detaches from the rest of the lattice at large Γ. This additional factor captures some of the effect of the change in the system's dynamical eigenmodes, given in Eqs. (8) to (10). These are significantly different from the system's original eigenmodes even at short and intermediate times.
The effective equation can be immediately solved. Taking the initial state to be the vacuum for simplicity, we find We observe that the two-mode correlations depend on their energy difference, and evolves in multiple stages: • Initially, while |ε i − ε j |t 1, the correlation grows linearly, b † i (t)b j (t) ∝ t.
• After a time approximately equal to the corresponding rate, |ε i − ε j |t 1, the correlation saturates with • For the duration of the intermediate time regime, until t ∼ N/Γ, these correlations are then independent of time up to a rotating term.
• Finally, at long times, further equilibration occurs, associated with the dissipative contribution to the dynamical mode eigenvalues (which are neglected here).
The anomalous correlations behave similarly with regards to the sum of the energies. In the intermediate regime, these behave similarly, coming into the light-cone at the same time and growing in the same way. In the long time regime, the dissipative dynamics drive this chiral system into a mirror-correlation steady state, and so the same-site correlations vanish while mirror correlations grow.
In Fig. 5 we plot the intermediate-time behavior of the energy eigenmode correlations of the simple system described in Section III B, the one dimensional ring with flux threaded through it. We observe the behavior outlined above.
We can also find the intermediate time correlations in real space. Within the approximation of Eq. (41), they are If the modes of the system are extended and have a plane-wave character, a ballistic behavior pattern emerges. If we take the wavefunctions to behave as for some propagation speed c, then â † m (t)â n (t) is simply the Fourier transform of the term in parentheses, evaluated at τ m , τ n . However, as these correlators are cut off at width t, while narrower features are largely defined by the time-independent post-saturation correlations. Thus, in the presence of a characteristic propagation speed c, we can expect a behavior of the form â † m (t)â n (t) = 0 max |m − n 0 |, |n − n 0 | ct Constant |m − n 0 |, |n − n 0 | ct.
We show this behavior for a one dimensional ring in Fig. 6, observing the behaviors outlined above. We observe a clear light-cone, with correlations cut off at a distance of 2Jt from the drain. Within this light cone, the correlations asymptotically approach a fixed pattern. Notably, this pattern is quite different from the steady state, exhibiting same-site anomalous correlations â nân ∼ â nâ−n which vanish in the long term.

V. CONCLUSIONS
We have analyzed here the dynamical properties of a bosonic lattice locally coupled to a single Markovian and Gaussian reservoir, with a particular focus on a squeezed bath. In calculating the the spectrum of its relaxation rates, we found that it is largely defined by the ratioΓ i /∆ i , crossing over from a perturbative coupling between the system's modes and the bath to a Zeno-like suppression when the drain is strongly coupled. We have also shown that the regime of intermediate-valued bath couplings can exhibit rich behavior, including resonant amplification of some modes' relaxation rates. Using Eqs. (22) and (23) this behavior can be analytically calculated for a given system, allowing these parameters to be chosen to e.g. optimize the over relaxation time.
We have also explored the intermediate time behavior of such systems in the case of a large lattice, proving the existence of a distinct intermediate time regime. There, we found that a pre-thermalized correlation pattern emerges, which may be quite different from both the initial state and steady state.
The local squeezing bath we suggest here could be experimentally realized as a modification of existing experiments [29,30]. As such, the dynamics we have calculated could be directly observed, and may prove useful, e.g. in generating desired entanglement in microwave cavities. Further modifications, such coupling to a non-Markovian bath or the addition of particle interactions, may reveal even richer physics in such lattices.
These have satisfying Eq. (7). Then, from Eq. (4) we finḋ and so we may immediately integrate to findb i (t).