Engineered Swift Equilibration for Brownian objects: from underdamped to overdamped dynamics

We propose a general framework to study transformations that drive an underdamped Brownian particle in contact with a thermal bath from an equilibrium state to a new one in an arbitrarily short time. To this end, we make use of a time and space-dependent potential, that plays a dual role: confine the particle, and manipulate the system. In the special case of an isothermal compression or decompression of a harmonically trapped particle, we derive explicit protocols that perform this quick transformation, following an inverse engineering method. We focus on the properties of these protocols, which crucially depend on two key dimensionless numbers that characterize the relative values of the three timescales of the problem, associated with friction, oscillations in the confinement and duration of the protocol. In particular, we show that our protocols encompass the known overdamped version of this problem and extend it to any friction for decompression and to a large range of frictions for compression.


I. INTRODUCTION
Shortcuts To Adiabaticity (STA) emerged in quantum mechanics as fast protocols for state-to-state transformations that would otherwise require the slow and therefore time-consuming modification of a control parameter of the system to reach the desired final state following a quasi-adiabatic trajectory [1]. Many strategies have been proposed to set up non-adiabatic routes to reach the same final state through the use of dynamical invariants [2], counter adiabatic driving [3][4][5], reverse engineering methods [6][7][8], fast-forward techniques [9,10], Lie algebraic approaches [11,12], and optimal control [13][14][15][16] to name but a few. Slow processes (adiabatic in quantum mechanics jargon) and thus STA are quite common to prepare the state of the system in a wide variety of domains including atomic and molecular physics [17,18], quantum transport [19][20][21], solid state [22], many-body physics [23][24][25], classical mechanics [26] and statistical physics [27][28][29]. STA also have applications in the design of optimal devices, as recently proposed in optics [30] and in internal state manipulation for interferometry [31]. STA therefore enjoy a large domain of applications, and the number of experiments demonstrating their efficiency is soaring.
Recently, these techniques have given birth to new protocols in statistical physics. Thermodynamic transformations that connect two different equilibrium states are not in most cases quasi-static and thus necessarily visit out-ofequilibrium states. Operating such transformations in a finite and short amount of time, potentially much shorter than the relaxation time of the system, is crucial for many applications, in particular in micro and nano devices or engines [32][33][34][35][36][37], triggering a number of works considering how STA could boost engines, among which [38][39][40][41]. As in quantum physics, performing this kind of quick transformations requires to devise an appropriate driving of the intermediate out-of-equilibrium dynamics. A recent example has been provided with protocols to compress or decompress an isolated 3D harmonically trapped cloud of atoms in an arbitrarily short amount of time [27]. More importantly, such an approach has been generalized to systems in contact with a thermostat which are ubiquitous in thermodynamics. The so-called Engineered Swift Equilibration (ESE) protocols have been introduced to study the isothermal compression of a colloidal particle [28]. The idea was then adapted to encompass the shift of the cantilever of an Atomic Force Microscope [29], and has been generalized in [42] to underdamped processes using a non-conservative driving force. Here, our interest also goes to the underdamped dynamics of a Brownian particle, under the proviso that the driving force, used to manipulate the system but also to confine the particle, is conservative.
In section II, we address the general case of ESE protocols in the underdamped regime, for a non-isothermal transformation and generic potentials. Contrarily to the works [40,42], we resort to conservative drivings, through potentials that only depend on space and not on velocity. This general framework leads to lengthy equations, that are significantly simplified in the case of transport-free harmonic potentials. We next restrict to this class of transformations in section III and show how to obtain a fully explicit isothermal protocols, choosing the shape of one characteristic quantity of the particle density function and deducing from it the appropriate evolution of the control parameter. Finally, in section IV we exhibit and analyze thoroughly the "phase diagram" of such protocols. We proceed to show that it largely depends on whether the transformation is a compression or a decompression, and work out the various properties of these protocols, such as existence, cross-over to the overdamped regime, position-velocity decoupling or also transient negativity of the stiffness. We supplement this by a discussion on the shape of the temporal evolution of the stiffness, through the comparison between the relevant timescales of the problem. We also analyze the robustness of the phase diagrams with respect to the shape of the protocol, and comment on the change occurring in the protocol when its duration is decreased. We conclude in section V.

II. GENERAL FORMALISM
The ESE protocol brought to the fore in [28] addressed the case of an overdamped confined Brownian object. While the overdamped limit is suited for colloids in a solvent like water, it is desirable to study the generalization of the idea to underdamped situations, when inertial effects no longer are negligible, such as for an Atomic Force Microscope tip where friction is on purpose reduced as much as possible [43], or for the study of a levitated nanoparticle in air where friction can be tuned through gas pressure [44]. Generically, when viscous friction is not high compared to the other characteristic frequencies of the problem, one should include the velocity degrees of freedom in the description in addition to positional ones; the overdamped approximation, on the other hand, assumes that the former are equilibrated at all times. To extend the ESE method proposed in [28] to the underdamped description of an object immersed in a thermal bath trapped in a confining potential, we introduce the probability density function K(x, v, t) of the position x and velocity v of the particle. It obeys the Kramers equation where U (x, t) is the confining potential, mγ the damping coefficient in the fluid with m the mass of the particle, k B the Boltzmann constant and T the temperature of the bath. At thermal equilibrium, this probability density function is simply given by the Boltzmann law To connect an initial equilibrium state characterized by the potential U i (x) and the temperature T i to a final equilibrium state (U f (x), T f ), we assume that the probability density function keeps a Gaussian form in v during the transformation This ansatz, inspired by previous works on Boltzmann equation where it results from the use of Boltzmann Htheorem [27], remains operational here. In Eq. (3), 1/B plays the role of a kinetic temperature, that should at equilibrium coincide with that of the bath (T ), but is otherwise distinct. It is worth emphasizing here that the bath temperature can be time-dependent. In the colloidal realm for example, this is achieved by an appropriate random shaking of the confinement potential [45,46], which creates an effective temperature for the Brownian object while the true bath temperature remains constant.
In the spirit of ESE techniques, we do not impose the control function/parameter U (x, t) and T (t) beforehand to study the response of the system through the functions A, B and D. On the contrary, we adopt a reverse engineering point of view, namely we choose a desired dynamics for these functions and deduce from it the temporal evolution of the control parameters that needs to be enforced to perform the chosen dynamics. Functions A, B, D and control parameters U and T are linked via a set of equations that we obtain by plugging the ansatz form (10) into the Kramers equation (1) and sorting the v 3 , v 2 , v and constant terms where the dot stands for time derivative. This set of equations, general within the ansatz (3), must be obeyed to connect the two imposed equilibrium states. Of course, they have to be adapted to specific protocols with desired constraints, such as duration of the protocol, amplitude of the transition and number and nature of control parameters. The knowledge of the existence of at least one of these specific protocols is in general not a simple task, as we will see later. However, let us focus on these equations and try to extract as much information as we can. They describe Gaussian compressions and decompressions with transport when U (x, t) remains quadratic, but the x part of the particle density function can also stray from the Gaussian shape when arbitrary potentials are used. For clarity and simplicity purposes, we will restrict ourselves to harmonic potentials of angular frequency ω(t) (stiffness κ = mω 2 ) and center position x 0 (t) so as to carry out a compression or a decompression from the initial state characterized by (ω i , T i ) to the final state (ω f , T f ). With this harmonic potential, our Kramers equation is linear in x and v, ensuring that the particle density function keeps a Gaussian shape at all times. A first analysis of equations (4a) and (4b) shows that the kinetic temperature 1/B is time-dependent only and does not depend on space, and that D is linear in x Integrating equation (4c) with respect to x then yields a quadratic form in x for A(x, t) where the function A 0 (t) is related to the normalisation of the distribution. Its temporal evolution can be left aside from our study, as our Fokker-Planck equation conserves probability during the transformation. If we plug expressions (6) and (7) into equation (4d), and sort out the monomials in x, we obtain the equations controlling the time evolution of B(t) ...
as well as the time evolution of D 0 (t) Both equations are informative. The first one is not very tractable in itself, but highlights that the inverse kinetic temperature B is completely determined by the bath temperature T and the trap stiffness ω, and thus independent of the transport part x 0 (t) of the transformation. On the other hand, the second equation tells us that D 0 (t) is fully induced by transport. With no transport, x 0 (t) vanishes during the whole transformation, and so does D 0 (t) since initially D 0 (0) = 0. As a result, D(x, t) is simply proportional to x, referring to equation (6). Finally, coming back to equation (7), in a transport-free case, A(x, t) is itself simply proportional to x 2 .

A. Simplified formalism
We now restrict to transport-less transformations, carried out by a harmonic potential. A variant of this problem was numerically solved in [47], and the corresponding protocol displayed discontinuities at initial and final times.
Here, we will provide an exact explicit solution while imposing smooth boundary conditions, in order to create a protocol well-adapted to experiments. Our Kramers ansatz can be written in the lighter form with the correspondence A(x, t) = α(t)x 2 − ln N , where N (t) is a normalization factor, B(x, t) = β(t) and D(x, t) = δ(t)x, the function δ(t) being the amplitude of x − v correlations. The set of equations (4) then comes down to with the following initial conditions For the sake of simplicity, we now rescale all quantities and variables as follows We distinguish the rescaled quantities from the corresponding dimensioned ones by a tilde. From now on, the dot stands for derivative with respect to rescaled time s. This rescaling yields the following seṫ with initial conditions We introduced the two dimensionless parameters N γ = γt f and N ω = ω i t f that appear independently in Eqs. (14). They turn out to be two key parameters in the discussion within the underdamped framework. Indeed, most of the physics of the problem stems from the comparison between the characteristic timescales, namely the time associated with viscous friction 1/γ, the period of oscillation in the harmonic potential 1/ω and the duration t f of the protocol, or rather the position and velocity relaxation times t x = γ/ω 2 i and t v = 1/γ, as will be discussed later on. The two numbers N γ and N ω then simply compare the duration of the protocol with respectively the viscous time and the duration of an oscillation in the potential. Note that for clarity, we use the initial oscillation time to define N ω , as it would otherwise vary during the transformation, preventing us from performing a general analysis.
In an ESE approach, as highlighted above, we choose the dynamics of some of the parameters of the particle density function (here among α, β and δ, or combinations of them) and deduce from them the required temporal evolution of the control parameters: here κ (and possibly also the bath temperature [46]), that can be controlled experimentally. The three relations (16) include five unknowns ( α, β, δ, κ and T ). Therefore two of the three functions α, β and δ can be chosen freely, while the last function and the two control parameters are extracted from the equations. This choice is not unique, which is often an advantage, for it opens for a great versatility of STA or ESE protocols; we will next see two particular ways to proceed.

B. Specification for a fixed temperature bath
The equations that control the system couple non linearly the functions α, β and δ and the control parameters κ and T , turning out to be rather complex to solve explicitly. As building an explicit protocol is desirable for the theoretical study of its properties as well as for the experimental implementation of the transformation, we further restrict our investigation to isothermal transformations T = 1. In this case, our set of three equations becomeṡ and the initial conditions are where χ = κ f /κ i is the compression factor. Although the bath temperature T is constant, the kinetic temperature 1/ β is in general not. The overdamped limit of this isothermal case, where the kinetic temperature is supposed to stay at equilibrium at all times and consequently discarded from the treatment, was previously addressed in [28]. Here, in the underdamped regime, the problem is intrinsically more complex as both position and velocity distributions (via the functions α, β and δ) need to be engineered only through the stiffness κ.
In this particular isothermal case, we rephrase our set of three equations in a way that naturally leads to an explicit expression of the whole protocol. We first introduce the quantity a and its rescaled equivalent a It directly relates to the position-variance (width of the marginal distribution of the position of the particle) Consequently a is a measurable quantity, whereas the physical information borne by α is more elusive, and it is advantageous to work with a rather than α. Using equations (16), it is easy to obtain a differential equation on ȧ Then we extract δ from Eq. (16b) and get˙ It is convenient to define in order to recast equation (21) as˙ It turns out that all the used quantities can be expressed in terms of ∆ and its derivatives. Indeed, from equation (23), β reads and a follows from the definition of ∆ (equation (22)). Then the expression of δ stems from equation (20 and α from its definition (16a) together with equations (24) and (25). Finally, the control parameter ω(s) can be read on equation (16a) κ =˙ α C. Towards explicit protocols The above hierarchy of equations is convenient to devise explicit protocols, and they consequently fully qualify as STA or ESE in spirit. Indeed, the starting point provided by Eqs. (16) is not convenient, for these coupled equations involve the unknown forcing time-dependent term κ, a function that needs to obey subtle properties to be compatible with the boundary conditions (17). Thus, as such, the problem is not amenable to numerical solution, since it is of a functional-shooting type: find the proper family of κ(s) enforcing the desired final condition. On the other hand, equations (22), (24), (25), (16a) and (26) pave the way to a simple analytical solution: we first choose the shape of ∆(s) and subsequently deduce from it the evolution of β making use of (24), from which a is known invoking (22); δ then follows from (25), α from (16a) and the desired forcing is computed, at the end of the chain, with Eq. (26). That solution is referred to as Protocol A.
We emphasize that special attention has to be paid to the temporal boundary conditions, as in every shortcutto-adiabaticity-type procedure, in order to avoid excitations of the system when reaching equilibrium, at the end of the protocol. In addition, a protocol that is smooth enough at initial time, when launched, is more conducive to a successful experimental realization. We thus enforce the same boundary conditions at initial time s = 0 and final time s = 1, but it can be kept in mind that a different choice can be made at s = 0, under the proviso that the boundary conditions at s = 1 are as above. Equations (16) imply that the first derivative of α, β and δ vanishes at initial and final times, and so does the first derivative of a. This condition on the first derivative of δ also forces the second derivative of a to be zero, through equation (25). In turn, this imposes the same condition on the second derivative of β. Altogether, this requires at least that the first three derivatives of ∆(s) are zero at the final time. Finally, the values of ∆ at initial and final times are simply ∆(0) = 1 and ∆(1) = χ. For the sake of simplicity, we choose a polynomial for ∆. The lowest order admissible function reads This method is straightforward, and singles out ∆, over which "control" is exerted. It is a combination between the width of the position distribution and the kinetic temperature (related to the velocity distribution), and thus a rather "secondary" quantity. Our goal is next to present an explicit variant, where another quantity is controlled, with more direct physical meaning.
Another route towards an explicit protocol consists in choosing the inverse kinetic temperature β(s) and deducing from it the other functions, through the same aforementioned hierarchy. The temporal boundary conditions required for β(s) are again β(0) = β(1) = 1 and the first two derivatives vanish at initial and final times. Once β(s) is chosen, a(s) is obtained integrating equation (21) a(s) = 1 while the other quantities can be expressed in terms of functions a and β. The initial condition a(0) = 1 is fulfilled since β(0) = 1, but the final condition a(1) = χ imposes an additional integral constraint on the chosen β(s) For the previous procedure, only ∆ and its derivatives were employed to express all the other functions, so that specifying the boundary conditions of ∆ was enough to ensure the right initial and final states. Here, the solving procedure involves an integral of the chosen function β(s) in equation (28). As a result, the initial and final states cannot be both encoded in the temporal boundary condition of β and yield the integral constraint (29). Note that this relation is true whatever the protocol, and is in particular automatically verified for Protocol A when ∆ has smooth enough boundary conditions. We coin this variant "Protocol B". Finally, a more natural quantity to choose and control would be the inverse variance of position a(s). However, there is no straightforward solution of the equations in terms of a(s), as there was in terms of ∆ or β. Moreover, a numerical resolution of the hierarchy of equations in terms of a necessarily involves a numerical integration of a differential equation, for example equation (28). Even a careful choice of the boundary conditions of a fails to produce reliably the desired initial and final states for the transformation. We meet again a functional shooting problem, where the function a(s), which is a parameter of the following equatioṅ has to be tuned so that the solution ∆ fulfills the desired initial and final conditions. Equipped with protocols A and B, we are nevertheless in a position to discuss the robustness of our main findings.

IV. CHARTING OUT THE PHASE DIAGRAMS OF THE PROBLEM
We now study the characteristics of protocol A, as a function of the two dimensionless parameters N γ and N ω . We focus in particular on the very existence of the protocol, on the applicability of the overdamped limit, on the decoupling of the x and v degrees of freedom, and on the temporal evolution of the only control parameter of the problem, namely the stiffness κ of the trap. We represent this thorough study on two phase diagrams in Figures 1  and 2. They of course depend on the characteristics of the desired transformation, as will be discussed later on; we have chosen here a compression factor χ = 2 for Figure 1 and a decompression factor χ = 0.5 for Figure 2.
A. Zones of the phase diagrams

Existence of the protocol
We first investigate the existence of the underdamped protocol itself. For the underdamped ansatz (10) to be well-defined, the functions α and β, as well as a, must remain positive during the whole transformation; our ansatz is otherwise divergent. However, the expression (24) shows that β is positive only as long aṡ This criterion can always be satisfied in the case of a decompression, where ∆ can be chosen monotonically decreasing whatever the value of N γ , while ∆ is by construction always positive. On the other hand, for a compression, low values of N γ can make the function β become negative, resulting in a diverging ansatz. This happens typically when N γ is of order 1. The exact threshold of course depends on the chosen shape of ∆ and the compression factor χ, as will be discussed later. Moreover, since the product ∆ = a β is fixed and positive, a and β are always of the same sign. As for α, it is also positive when β (and thus a) is positive too. Therefore, the only region of the compression phase diagram where the ansatz is ill-defined is the left half-plane under a threshold N min γ of order unity. The ansatz is on the other hand well-defined in the whole decompression phase diagram.
We illustrate this on Figure 1, where a and β are plotted near the boundary of the existence region, on the red left inset curves. We notice that the minimum of β tends quite quickly to zero when N γ is decreased around one and Here, the compression factor is χ = 2, ∆ is given by equation (27), and the set of (Nγ, Nω) parameters is indicated by the position of the dots on the phase diagram. The straight thin black line represents a "trajectory" followed in the phase diagram when the duration t f of the protocol is decreased, all other parameters being fixed (see Fig. 5 for further details). The stars correspond to the different curves presented in Figure 5.
that in turn, a becomes very large in this region. When the inverse kinetic temperature β goes to zero, the velocity distribution of the particle becomes very broad and the kinetic temperature diverges. In order to keep a control on the particle, the stiffness of the confinement has to increase dramatically, yielding a very peaked a, as illustrated. In this region, where the protocol can be considered fast since it no longer allows for velocity equilibration, it is necessary to use large stiffness to obtain the desired compression, which provides work to the system. This results in a heating, marked by the drop of β.
Our approach fails when the duration t f of the compression protocol is of the order of the friction time 1/γ (see Figure 1, left area). This can be overcome by using a non-conservative potential that would penalize high velocities and then prevent this dramatic increase of the kinetic temperature, as suggested in [40,42]. On the other hand, for a decompression, the kinetic temperature cannot diverge during the transformation, and our approach remains valid at any friction in this case, without needing a non-conservative potential. We stress that realizing such forces in an experiment is a challenge, while the conservative case worked out here is routinely employed with optically confined colloids (see [45] and references therein).

Overdamped limit
Next, we consider the overdamped limit and show that we can recover the results obtained in [28] from the underdamped formalism developed here. We also determine the regime of validity of the overdamped approximation, which amounts to treating the velocity degrees of freedom as equilibrated at all times. As their distribution relaxes to the Gaussian distribution on a timescale 1/γ, the overdamped approximation requires this time to be much smaller than the other timescales of the problem; in our notations, this leads to N γ 1 and N γ N ω . Note that in the following, ∆(s) is kept unspecified for the sake of generality, but is always of "order 1", in that it is chosen independently of N γ and N ω . Moreover, the following discussion holds both for a compression and a decompression.
To investigate the overdamped limit, we focus on the temporal evolution of the stiffness of the potential κ. Starting from the overdamped formalism, it was found in [28] to be related to the function a(s) of equation (19) through Therefore, we concentrate on the quantity κ − a in the regime N γ 1 and N γ N ω . Using equation (16a), (18) and (20), we obtain This expression can be developed in power series of N γ . We emphasize that in the limit of large N γ , β remains close to one (as expected, the kinetic temperature is equilibrated at all times with the bath temperature), and a is then approximately ∆, which does neither depend on N γ nor on N ω . This yields and a β˙ a 2 δ˙ δ Finally under the assumption N γ N ω , we can compare the leading term of equations (34) and (35) a˙ ∆ and we obtain This matches exactly equation (32) in the overdamped limit, as expected when N γ N ω and N γ 1. This expression is also informative on the scalings with respect to N γ and N ω of the corrections to the overdamped asymptotics. Yet, a subtlety should be outlined here, and it is related to the polynomial shape chosen for ∆(s). When our parameters qualify the protocol as belonging in the overdamped regime, a inherits this functional form, which is slightly more complex than the one chosen in [28], where the smoothness requirement only concerns its first derivative. Indeed, the underdamped protocol requires more derivatives to vanish than its overdamped counterpart as discussed in paragraph III C, which results in a polynomial ∆(s) of higher order. Note here that comparing directly the stiffness κ computed within the full underdamped formalism to the overdamped one, instead of focusing on κ − a, yields a looser criterion: N γ 1 without any condition on N ω , as represented on the phase diagrams 1 and 2 in hatched green.
To conclude this discussion, we point out that our criterion for the overdamped limit differs from the one given in [42], which can be rephrased as N ω N γ . This discrepancy likely comes from the fact that the authors impose that the temperature of the particle (our inverse β) remains constant during the transformation. Though this choice offers much simpler calculations, it requires to implement a potential that is quadratic in impulsion, therefore creating a set-up that strongly strays from ours, where the protocol only involves conservative drivings.

x − v decorrelation
An interesting characteristics of the protocol is whether or not position and velocity degrees of freedom are correlated. This correlation is measured by comparing the (squared) cross term δ 2 to the product of position and velocity variances α β. In particular, these correlations vanish when a collapses onto α. Following the definition of a in (18) and of ∆ in (22), this amounts to the criterion | δ| 2 | ∆|. Going back to expression (25) for δ, we see that its scaling depends on that of β. At large N γ , β 1 so that δ scales as 1/N ω . On the contrary, for low N γ (that only concern decompressions as discussed in the paragraph about the existence of the protocol), β scales as 1/N γ so that δ scales as 1/ (N ω N γ ). The criterion | δ| 2 | ∆| then takes two different shapes in the limits of high and low N γ . For high N γ , x and v degrees of freedom are decorrelated for N 2 ω 1, which is confirmed by the two phase diagrams (see Figures 1 and 2). On the other hand, for low N γ and for decompressions, decorrelation arises when N ω 1/N γ . This is also confirmed by the decompression phase diagram (Figure 2), where the left frontier of the decorrelation area has a slope −1 in double log scale. The border of this region is determined numerically as the curve on which the relative difference between a and α is of one percent.  Figure 1, and the compression factor χ is 0.5. In a symmetric fashion compared to quick compressions where the stiffness exhibits strong overshoots, quick decompressions rely on a stiffness that involves a strong undershoot (see the three insets). Exceeding the target value of the stiffness (upwards for a compression, downwards for a decompression) is a means to accelerate the transformation. For very quick transformations, the shape of the stiffness can become very complex and display two undershoots as illustrated by the left two insets. Again, the straight thin black line represents a "trajectory" followed in the phase diagram when the duration t f of the protocol is decreased, all other parameters being fixed, and the stars correspond to the different curves of Fig. 5.

Implementation challenges and ESE relevance
Finally, a primordial aspect of a protocol such as devised here lies in the characteristics of the control parameter κ that needs to be enforced experimentally to achieve the desired transformation. The main experimental challenge arises when the stiffness of the trap becomes transiently negative. In this case, the potential switches from confining to repulsive, a feature that cannot be achieved simply with an optically trapped colloid. We determine numerically the part of the phase diagram where this change in the sign of the trap curvature takes place; it is represented in purple on Figs. 1 and 2.
The global behavior of the stiffness has a complex dependence on the physical parameters N γ and N ω . We capture this diversity with Figure 3, where we represent the shape of the stiffness as a function of rescaled time s for representative points of the phase space (N γ , N ω ). From the phase diagrams of Figs. 1 and 2, we only keep the features that concern the implementation of the protocol, namely the non-existence area for compressions, and the aforementioned zone where the stiffness is transiently negative. Insights into the zoology of behaviors of the stiffness can be gained by comparing the three timescales of the problem, namely the duration of the protocol t f , the timescale of position relaxation t x = γ/ω 2 i and the timescale of the velocity relaxation t v = 1/γ, which are combined in N γ and N ω . Note that the maximum of these last two timescales defines the global relaxation timescale of the problem. On Figure 3, the dashed lines indicate where these timescales are equal two by two, and the different colors stand for the six possible orders of the three timescales. We emphasize that for ESE purposes, the two red areas, where the protocol is slower than the global relaxation scale of the problem, are less interesting than the green areas where position or velocity (or both) relaxes more slowly than the protocol. The most interesting set of parameters (N γ , N ω ) ESE-wise are therefore the green zones of these two diagrams out of the grey-shadowed areas.
Determining the shape of the stiffness from general considerations on our explicit protocol is challenging, as it involves several quantities that themselves depend non trivially on N γ and N ω (see Eq. (26)). However, we can grasp some features of the shape of the stiffness with the help of the diagrams of Fig 3. First, the stiffness has a very smooth behavior in the red zones, as expected, where the protocol is slow compared to the intrinsic dynamics of the system. At large N γ , in the region where the stiffness is appropriately described by its overdamped expression (37) as discussed in section IV A 2, the dependence on the parameters N γ and N ω of the stiffness required to perform the transformation is proportional to N γ /N 2 ω (see Eq. (32)), i.e. t x /t f . Then the region where the amplitude of the stiffness is large compared to 1 corresponds to the part of the diagrams that is well below the line t f = t x , as shown on the bottom right insets for compression diagram and most of the bottom insets for the decompression diagram. The farther from the line t f = t x , the greater the amplitude of the stiffness. This tendency to have large variations in the stiffness is also amplified when t f becomes smaller than t v , as shown on the two bottom left insets of the decompression diagram. Finally, the regions in dark green, where ESE protocols are highly desirable because shorter than both position and velocity relaxation times, display extreme shapes of stiffness. These bizarre shapes result from the difficulty to control the evolution of both position and velocity, which cannot equilibrate themselves that quickly, with only one control parameter.

B. Robustness of the phase diagrams
We have analyzed above protocol A phase diagrams, in which some details can depend of the functional choice made for the chosen ∆(s). We now quickly address protocol B properties, as a means to put to the test the robustness of our findings. In this discussion, we choose to only address protocols where the kinetic temperature is temporarily increased in the case of a compression (resp. decreased in the case of a decompression). In other words, we suppose that β(s) is always smaller than one (resp. bigger than one). This is equivalent to restricting to monotonous ∆(s), as indicated by equation (23). As ∆(0) = 1 and ∆(1) = χ, this function remains bounded, irrespective of the values of the parameters N γ and N ω . The previous discussions about the overdamped limit and the decorrelation of position and velocity degrees of freedom are therefore still valid. These areas are then completely robust with respect to the exact shape of the protocol.
On the other hand, the exact position of the region where the ansatz no longer exists depends on the shape of the chosen function β for protocol B or ∆ for protocol A. The corresponding boundary is always a vertical line in the compression phase diagram, as equation (23) does not involve N ω . The most favorable case, in which the existence domain is the largest, corresponds to a situation where β(s) is flat during almost the whole compression (except near initial and final times to fulfill the temporal boundary conditions). Together with the fixed integral (29) for β, this indicates that the protocol-dependent existence threshold N min γ is always such that Therefore this non-existence zone is a generic feature of the process and cannot be eliminated by tinkering with the shape of the protocol. In particular, there is no way to devise a very fast isothermal compression protocol: t f > γ −1 (ln χ)/2. This is a strict lower bound, that can be significantly exceeded in cases where ∆ is non-monotonous (protocol A), or equivalently when β presents an overshoot (protocol B). The integral of β(s) being fixed by equation (29), such an overshoot necessarily needs to be compensated by lower values of β(s) in a different part of the transformation, therefore approaching the forbidden zero value. This may endanger the convergence of the ansatz and thus provide a further reason to dismiss non-monotonous ∆. As for the rather exotic shapes of the stiffness κ discussed earlier, we show in Figure 4 the differences between protocols A and B. As expected from the hierarchy of equations, a slight modification in the driving function yields significant changes, affecting all other quantities, including κ. However, Fig. 4 indicates that the essential features reported in the phase diagrams are robust with respect to a protocol change.

C. Consequences of accelerating the protocol
We finally address the effect of accelerating the protocol (diminishing t f , other parameters being fixed). How do the functions a(s), β(s), δ(s) and most importantly the stiffness κ(s) evolve? The results are reported in Figure 5 (a). For this "cut" across the phase diagrams (line of slope 1), we choose N γ = N ω , i.e. ω i = γ (straight thin black line in Figs. 1 and 2). This corresponds to the situation where t x = t v .
In Fig 5 (a), the purple curves represent a very slow transformation. As expected, the stiffness κ(s) interpolates smoothly between the boundary values, and is well followed by the inverse variance of position a(s): the dynamics is indeed slower than the timescales t x and t v . The velocity distribution always remains at equilibrium with the inverse kinetic temperature β(s) that stays close to 1. Velocity and position degrees of freedom are decoupled, as shown by a very low crossed term δ(s).
When the protocol duration t f decreases, the inverse kinetic temperature deviates, temporarily but more and more, from its unit equilibrium value, undergoing a transient heating for a compression and a transient cooling for decompression. As explained in paragraph IV A 1, the protocol no longer exists for too fast a compression, since β(s) becomes temporarily negative, preventing the ansatz for the particle density function from being well-defined. We also observe a change of sign in the crossed term δ during the transformation. As this quantity is proportional to the derivative of a(s) (see equation (25)), a change of sign transposes into a change of monotony of a(s), as illustrated in Figure 5 (b). This can be understood as follows. We first recast our ansatz (10) as so as to bring out a local mean velocity v loc (x) = −δx/(2β). The sign of δ then indicates the tendency for the local velocity, as sketched in Figure 5 (b). This yields a compression if δ is positive and a decompression if δ is negative. Indeed, the variance of the position decreases for a positive delta and increases in the opposite case (see Figure 5 (a)). . The left column is for a compression (χ = 2) and the right column is for a decompression (χ = 0.5). Both columns are for protocol A (with Eq. (27)). The color code is explained in the upper logarithmic time arrow. The red curves are thus for the faster protocol, that we omitted for decompression, as it does not bring a different information from the orange curve. (b) Illustration of the relation between the sign of δ and the monotony of a. If δ is positive, the particle density function tends to compress whereas when it is negative, it tends to expand.
When t f diminishes, a(s) steepens, leading to a pronounced overshoot for a fast compression (resp. undershoot for a fast decompression). Finally, the corresponding trap stiffness turns into a complicated non-monotonic function for fast protocols, with negative portions as well as peaked variations, as discussed in paragraph IV A 4. This rather unexpected behavior stems from the fact that controlling the dynamics of both the position and the velocity of the system with only one control function, the harmonic force depending on positional degrees of freedom only, is a delicate task. Fig. 5 (a) highlights the difficulty faced when devising an ESE protocol, since even in a case where t f exceeds the natural relaxation time by a factor three (yellow curves), the driving force and associated response significantly depart from their quasi-static counterpart.

V. CONCLUSION
In this article, we provide a general framework to study Engineered Swift Equilibration beyond the overdamped regime in which it was initially formulated. A Brownian particle is here confined in a harmonic potential, the stiffness of which can be changed in time as desired. In addition, the thermal bath is allowed to have a time dependent temperature T . As surprising as this situation might appear, the latter T -control is achievable in the laboratory with, for instance, optically confined colloids [45,46]. Trap stiffness and temperature are the two driving functions, that need in general to be carefully shaped to meet the desired goal: reaching the target state at the end of a chosen time t f . Yet, the formalism becomes cumbersome when T is time dependent, and explicit solutions become elusive. ESE techniques being designed especially for experiments and concrete applications, it is crucial to be able to exhibit such an explicit protocol. For this reason, we restricted our discussion to harmonic isothermal transport-free transformations, that is to say compression and decompression, where the formalism gets significantly simpler. Trap stiffness is thus the only quantity that is monitored by the experimentalist.
We discussed the explicit and analytical methods that can be employed, and analyzed the corresponding protocols as a function of the two key quantities of the problem, formed by the ratio between the relevant characteristic timescales. We summarized this analysis in two phase diagrams and discussed the range of applicability of our approach, that is "limitation-free" in the case of a decompression and limited to protocols longer than some lower bound ruled by the friction time in the case of compression. We also investigated core characteristics of our protocol, such as the crossover to the overdamped limit (where algebra is much simpler). Some attention was also paid to the relevance of the ESE protocol for each set of parameters, leading to investigate the influence of the protocol duration compared to the relaxation timescales on the shape of the trap stiffness. Finally, we discussed the robustness of the phase diagrams presented, by comparing the outcome of two distinct protocols in a parameter range where rather exotic drivings emerge.
Interesting venues for future work include extending our treatment to baths with time dependent temperature. This additional driving degree of freedom presumably leads to more regular protocols. Roughly speaking, the stiffness will ensure the compression or decompression in position space while temperature will take care of the velocity degrees of freedom. We thereby expect to overcome the limitations brought to the fore here, such as the non-existence of the underlying ansatz (which may become un-normalizable) and the odd shape of the stiffness. A first evidence of such an experimental achievement is presented in [46] where an effective modulation of the bath temperature allowed to perform a quick decompression without having to resort to a transiently negative stiffness. Our work also opens interesting perspectives for transformations including transport. Finally, the study of non-harmonic driving, energetics, and optimal features appears timely.

VI. ACKNOWLEDGMENTS
This work as been supported by the ERC contract OUTEFLUCOP. We acknowledge funding from the Investissement d'Avenir LabEx PALM program (Grant No. ANR-10-LABX-0039-PALM). We also thank A. Prados, C. Plata, A. Chepelianskii and O. Dulieu for insightful discussions.