Superadiabatic thermalization of a quantum oscillator by engineered dephasing

Fast nonadiabatic control protocols known as shortcuts to adiabaticity have found a plethora of applications, but their use has been severely limited to speeding up the dynamics of isolated quantum systems. We introduce shortcuts for open quantum processes that make possible the fast control of Gaussian states in nonunitary processes. Speciﬁcally, we provide the time modulation of the trap frequency and dephasing strength that allow preparing an arbitrary thermal state in a ﬁnite time. Experimental implementation can be done via stochastic parametric driving or continuous measurements, readily accessible in a variety of platforms.

The fast control of quantum systems with high fidelity is broadly acknowledged as a necessity to advance quantum science and technology. In this context, techniques known as shortcuts to adiabaticity (STA) have provided an alternative to adiabatic driving with a wide variety of applications [1]. STA tailor excitations in nonadiabatic processes to prepare a given state in a finite time, without the requirement of slow driving. The experimental demonstration of STA was pioneered in a trapped thermal cloud [2], soon followed by implementations in Bose-Einstein condensates [3], cold atoms in optical lattices [4], and low-dimensional quantum fluids [5]. More recently, STA have also been applied to Fermi gases, both in the noninteracting and unitary regimes [6,7]. Beyond the realm of cold atoms, STA have been demonstrated in quantum optical systems [8], trapped ions [9], nitrogen-vacancy centers [10], and superconducting qubits [11,12]. Their application is not restricted to quantum systems and classical counterparts exist [13,14] of relevance, e.g., to colloidal systems [15].
A variety of related control techniques fall under the umbrella of STA. Prominent examples include counterdiabatic or transitionless quantum driving [16][17][18], the fast-forward technique [19,20], reverse-engineered dynamics using Lewis-Riesenfeld invariants [21], as well as the use of dynamical scaling laws [14,[22][23][24], Lax pairs [25], variational methods [26], and Floquet engineering [27]. The use of STA in the quantum domain is severely limited to isolated systems, in which sources of noise and decoherence are considered an unwanted perturbation [28,29]. Applications to finite-time thermodynamics have thus been limited to the speedup of strokes in which the working substance is in isolation and decoupled from any external reservoir [7].
Controlling heating and cooling processes would pave the way to the realization of superadiabatic heat engines and refrigerators based, e.g., in an Otto or Carnot quantum cycle [30][31][32][33][34][35]. Hence, the possibility to speed up the dynamics of open quantum systems is highly desirable in view of applications to cooling, and more generally, in finite-time thermodynamics. In this context, the nonadiabatic control of composite and open quantum systems using STA remains an exciting open problem on which few results are available [34][35][36][37][38][39].
In this work, we introduce STA with open dynamics and apply them to the superadiabatic cooling and heating of a thermal harmonic oscillator. We show that the required control protocols are local and involve only the driving of the trap frequency and the dephasing strength. They can be achieved using stochastic parametric driving, thus harnessing noise as a resource.

I. MODEL
We shall consider a single particle in a driven harmonic trap, with Hamiltonian and a density matrix evolving according to a master equation of the form the derivation of which will be provided below. The case with constant dephasing strength γ admits the Lindblad form with position operatorx as the single Hermitian Lindblad operator. This naturally arises as the high-temperature limit of quantum Brownian motion. The dynamics with an arbitrary time dependence γ t is generally non-Markovian. We shall show how a time-inhomogeneous Markovian dynamics [40], corresponding to γ t > 0, can be engineered by tailoring noise as a resource. The dynamics along the process is assumed to remain Gaussian, with a density matrix in coordinate space of the form where A t , B t , C t are time-dependent coefficients to be determined from the master equation, with N t = √ 2(A t + C t )/π being the normalization factor. This form includes coherences during the dynamics and represents a family of dynamical processes that, as shown below, allows for a fast and controlled thermalization.
As a relevant example, we consider the driving in a finite time t f of an initial thermal state, parameterized by the trap frequency and inverse temperature (ω 0 , β 0 ), to a different thermal state with (ω f , β f ). For the Gaussian variational ansatz (3) to describe the exact dynamics of the master equation (2), the following consistency equations are to be satisfied (see Appendix B): The parameter 2hB t /m = −(Ȧ t +Ċ t )/[2(A t + C t )] ≡ t is homogenous to a frequency and directly follows from these equations. The boundary conditions are given from the initial and final states that we choose to be thermal. As detailed in Appendix A, it follows that A 0 = mω 0 2h coth(hω 0 β 0 ), B 0 = 0, and C 0 = −mω 0 /[2h sinh(hω 0 β 0 )]. Similarly, for the final state to be thermal, the coefficients at time t f should reduce to the values Because the initial and final states are taken as equilibrium states, they are stationary. This imposes the additional boundary conditionsȦ 0/ f =Ḃ 0/ f = C 0/ f = 0. We also require thatÄ 0/ f =B 0/ f =C 0/ f = 0 at initial and final time-the latter conditions are auxiliary, but guarantee a smooth variation of ω t and γ t .
A protocol speeding up the evolution from the thermal state characterized by (ω 0 , β 0 ) to (ω f , β f ) is obtained by explicitly specifying both the time dependence of γ t and ω t , as directly given by the consistency equation (4), according to Engineering a shortcut to thermalization between Gaussian states thus requires the ability to control both the frequency and dephasing. The control of the harmonic frequency ω t is performed with routine in a variety of setups and has been used to implement STA in isolated quantum systems, e.g., with trapped ultracold atomic systems [2,3,[5][6][7]. The requirement of a time-dependent dephasing γ t makes the dynamics open. It can be experimentally implemented from the microscopic picture provided below.

II. ENGINEERING OF TIME-DEPENDENT DEPHASING RATES
To modulate the dephasing strength γ t > 0 in the laboratory, we propose two different strategies: (i) harnessing noise as a resource [41,42] or (ii) via continuous measurements, which have been implemented in, e.g., trapped ions [43] and solid-state qubits [44], respectively.
(i) The master equation (2) can be obtained from implementing the stochastic Hamiltonian characterized by the Wiener process W t = W 0 + t 0 ξ t dt defined in terms of the real Gaussian process ξ t . While such a stochastic process is not differentiable, all integral quantities can be defined from the Wiener increment dW t = ξ t dt. The noise-averaged expressions follow from the moments ξ t and ξ t ξ t , that we choose to be zero and δ(t − t ), respectively, to describe a real Gaussian white-noise process [45].
The evolution of a quantum state dictated by the stochastic Hamiltonian (7) is described by a master equation that we derive below. For a small increment of time dt, the wave function can be written as with dW t defined in the Itô sense, i.e., fulfilling (dW t ) 2 = dt and dW t dt = 0 [46][47][48]. A Taylor expansion of the exponential then gives the only nonzero terms being first order in dt or (dW t ) 2 . Further, in the Itô calculus, the Leibnitz chain rule This gives the evolution of the density matrix ρ st = |ψ t ψ t | as which preserves the norm at the level of each individual realization. We then take the average over the realizations of the noise and denote the ensemble ρ t = ρ st . Using the fact that the average of any function F t of the stochastic process vanishes, F t dW t = 0 [46], we find that the evolution for the ensemble density matrix ρ t as dictated by the master equation (2).
(ii) Alternatively, the same evolution can be induced via continuous quantum measurements [44,[49][50][51][52] whenever the strength of the measurement is time varying. Consider a quantum system subject to a continuous quantum measurement of the observableÂ. Its evolution is known to be described by the stochastic nonlinear master equation [50,51] where dW t denotes a random Gaussian real random variable of zero mean and variance dt. The characteristic measurement time with which observableÂ is monitored, denoted τ m , can 033178-2 be controlled by changing the measurement strength. The deterministic part of the evolution L(ρ st t ) includes a nonunitary term of the standard Lindblad form, while the so-called innovation term reads The latter is nonlinear in the state ρ st t and represents the measurement backaction on the system resulting from the acquisition of information during the measurement process. A specific trajectory is associated to a given realization of the Wiener process dW t , and characterized by fluctuations of the measurement outcomes given by dr A = Â (t ) + √ τ m dW t . When the observer does not have access to the measurement outcomes, the system is consistently described by the state ρ t = ρ st t , which results from averaging over an ensemble of trajectories, and that satisfies dρ t = L(ρ t )dt. So monitoring the position operator (Â =x) with a time-dependent measurement strength such that 1 8τ m = γ t (obtained, e.g., by applying feedback on the system [52]), effectively generates the master equation (2).
To sum up, the engineering of a prescribed modulation in time of the dephasing strength γ t can be achieved via stochastic parametric driving or continuous measurements, provided that γ t > 0. Interestingly, both techniques allow modulating γ t independently from the frequency ω t , which contrasts with the time-dependent Markovian quantum master equation derived by driving the coupling of a system to a thermal bath [53]. Our scheme can be readily implemented in a single trapped ion [54], in which the creation of an open dynamics with artificial environment [55,56] or via the addition of noise [43] have been experimentally demonstrated.

III. CHARACTERIZATION OF THE DYNAMICS
The evolving density matrix can be diagonalized at all time according to ρ t = n p n,t |ψ n,t ψ n,t |, the eigenvalues and eigenfunction being (see Appendix C and Ref. [57]) where H n denotes the Hermite polynomial defined from . The effective inverse length k t and dimensionless constant u t that characterize the control trap are detailed in Appendix C and below. Interestingly, the evolving density matrix ρ t can be interpreted as a thermal state σ t rotated through a unitary transformation U x,t ≡ e −iB tx 2 by noting that The density matrix σ t , with coordinate representation x|σ t |x = N t e −A t (x 2 +x 2 )−2C t xx , corresponds to the instantaneous thermal state of a harmonic oscillator with effective frequencyω t and inverse temperatureβ t provided assuming oscillators of equal mass. The effective inverse length is then explicitly given by k t = √ mω t /h. By construction, the two states share the same eigenvalues and σ t = n p n,t |n t n t |, with the probability now written in terms of a thermal probability at all times, p n,t = e −β thωt n /Z t , the partition function being Z t = 1/(1 − u t ), and u t = e −ε t . However, the eigenvectors are different and |n t =Û x,t |ψ n,t correspond to the well-known Fock states of the reference, time-dependent harmonic oscillatorH t -whose parameters are distinguished with a tilde.
At all times of evolution, we have A t = (k 2 t /2) cothε t and . So the control frequency and dephasing strength can be recast in the form where the control parameter depends on the scaling factor κ t ≡ k 0 /k t = √ ω 0 /ω t and temperatures as These are our main results. The combined modulation of the trap frequency and the dephasing strength is sufficient to engineer finite-time shortcuts to thermalization. Equation (16) gives the correction of the control of the trap frequency ω t with respect to a reference oneω t that needs to be experimentally implemented for the preparation of the thermal state in a finite, prescribed time. A comparison of these results with the ones reported for isolated systems with γ t = 0 [24] shows that the control parameter in Eq. (18) not only depends on the scaling factor but also accounts for the change of temperature through an additional, nontrivial term.
Our scheme can be implemented by choosing an interpolating ansatz between the boundary conditions imposed on the state, i.e., ω 0/ f and β 0/ f that define A 0/ f and C 0/ f . For illustration, we choose (19a) in the form of a fifth-order polynomial, f (τ ) = 10τ 3 − 15τ 4 + 6τ 5 , to ensure a smooth dynamics. The modulation of the control parameters ω 2 t and γ t readily follow from Eqs. (14)- (18). Figure 1 illustrates the control frequency and dephasing strength corresponding to phase-space compression and expansion protocols, discussed below. Short control processes require trap inversion (a negative squared frequency), which can be achieved experimentally via, e.g., a painted potential [58] or a digital micromirror device [59]. They also rely on a dephasing strength of larger amplitude, which can be experimentally more challenging to achieve. We 033178-3 propose to use the maximum of the dephasing strength, denoted γ max as a measure of the "cost" to implement the process by a technique such as stochastic parametric driving. We show in Appendix D that the maximum dephasing strength scales inversely with the process time, as illustrated in the insets of Fig. 1(b). We further use the relative entropy, defined as S(ρ t σ t ) = Tr(ρ t ln ρ t ) − Tr(ρ t ln σ t ), as a measure of the distance of the engineered state ρ t to the effective thermal state σ t along the dynamics. It can be written as (20) where the overlap of the eigenfunctions m t |n t is given explicitly in Appendix E following [60][61][62]. Figure 1 illustrates the relative entropy between the engineered and the thermal state. The insets show the maximum relative entropy for different process times, evidencing that the state is going further away from a thermal distribution for shorter protocols. The shape of the dephasing strength and the relative entropy is independent on the duration of the process, which only influences their maxima, γ max and s max , respectively, as illustrated in Fig. 1(b).

IV. SUPERADIABATIC PROTOCOLS
Processes satisfying β f ω f = β 0 ω 0 conserve the mean phonon number and are often referred to as phase-space (density) preserving. The inverse temperature and frequency can be related to two physical lengths, namely, the particle characteristic length, given by the de Broglie wavelength, l dB =h √ β/(2m), and the trap characteristic length, λ HO = √h /(mω). Their ratio l dB /λ HO = √ βhω/2 is conserved for phase-space-preserving transformation. STA in closed systems are limited to phase-space-preserving cooling techniques, such as adiabatic cooling. These processes preserve the von Neumann entropy S t = −Tr(ρ t ln ρ t ). By contrast, cooling and heating processes altering the phase-space density, and the number of populated states lead to an entropy change [63] and require an open dynamics. STA for open processes thus allow reaching arbitrary thermal states (ω f , β f ) from an initial thermal state, as schematically represented in Fig. 2, along with the variation of entropy. The sign of the dephasing strength determines the variation of relative energy sign(ε t ) and entropy change.
In particular, a positive dephasing strength yields a monotonic increase of entropy. Indeed, the rate of change of the von Neumann entropy reads Protocols restricted to γ t 0 allow only STA for thermalization to high-temperature states (heating), with S = S f − S 0 > 0. Whenever values of γ t can be engineered, this restriction is lifted. The maximum dephasing strength, illustrated in Fig. 3, is specific to each scenario. It follows a different behavior when changing the trap frequency at constant temperature or vice versa. This is not surprising since the two parameters correspond to different physical phenomena, as discussed above. The plateau observed when decreasing the temperature at a fixed trap frequency [ Fig. 3(c)] might be set by the trap size, which constrains the size of the particle. Interestingly, a given final phase-space density can be reached from different dephasing strengths, even when starting from a same initial state. In other words, processes yielding to β f ω f from β 0 ω 0 can have different implementation costs according to γ max , as illustrated in Fig. 3(d).
In conclusion, we have introduced shortcuts to adiabaticity with an open dynamics engineered to control the thermalization of a quantum oscillator. The resulting protocols are expected to be broadly applicable as their implementation requires only a time modulation of the harmonic frequency and the dephasing strength, accessible, e.g., via stochastic parametric driving or continuous quantum measurements. Our results can be directly applied to non-Markovian dynamics whenever the amplitude and sign of the dephasing strength can be engineered. Extension to obtain generalized Gibbs states and for interacting systems is under investigation.

ACKNOWLEDGMENTS
It is a pleasure to acknowledge discussions with Ángel Rivas and Bijay K. Agarwalla.

APPENDIX A: THERMAL STATE OF A HARMONIC OSCILLATOR
It is well known that the thermal state of a harmonic oscillator is Gaussian in the coordinate representation. For the sake of completeness, we briefly sketch the derivation below. For a time-independent harmonic oscillator, the HamiltonianĤ = where k −1 = √h /(mω) denotes an effective length characteristic of the harmonic oscillator. The coordinate representation of the thermal density operator then reads We use Mehler's formula [57], to rewrite the sum with the Hermite polynomials as a Gaussian, yielding 033178-5 where we have defined u = e −βhω . Finally, with the explicit form of the partition function Z = Tr(e −βĤ ) = 1/(1 − u), this coordinate representation also takes the form The same derivation holds for a time-dependent Hamiltonian and provides the initial and final coefficients A and C given in the main text. We verify that the normalization factor is k

APPENDIX B: CONSISTENCY EQUATIONS FROM THE EVOLUTION OF THE GAUSSIAN ANSATZ
The master equation (2) in the coordinate representation reads For the Gaussian ansatz given in Eq. (3) of the main text, the real and imaginary parts of the evolution equation respectively givė from which the consistency equations (4) directly follow.

APPENDIX C: INSTANTANEOUS DIAGONALIZATION OF THE DENSITY MATRIX
We look for the eigenvalues p n,t and eigenfunctions |ψ n,t that diagonalize the density matrix ρ t at any time. For the sake of simplicity, we omit the time dependence in the notation below. By definition, the eigenvalues fulfill p n > 0 and n p n = 1, so we choose to write them as p n = u n (1 − u), where u can be seen as an exponential e −ε . We verify below that the functions By identification, we obtain the reverse transformation corresponding to the physical setting being for u > 0 and 0 < −A/C < 1, and The mean phonon number easily follows as n = ∞ n=0 np n = u 1−u , and the von Neumann entropy S(ρ) = −Tr(ρ log ρ) reads