Shortcuts to Squeezed Thermal States

Squeezed state in harmonic systems can be generated through a variety of techniques, in-cluding varying the oscillator frequency or using nonlinear two-photon Raman interaction. We focus on these two techniques to drive an initial thermal state into a ﬁnal squeezed thermal state with controlled squeezing parameters— amplitude and phase—in arbitrary time. The protocols are designed through reverse engineering for both unitary and open dynamics. Control of the dissipation is achieved using stochastic processes, readily implementable via, e.g., continuous quantum measurements. Importantly, this allows controlling the state entropy and can be used for fast thermalization. The developed protocols are thus suited to generate squeezed thermal states at controlled temperature in arbitrary time.


Introduction
Squeezing is a paradigmatic quantum effect that allows reducing fluctuations of one variable beneath the standard quantum limit. This is achieved at the expenses of increasing the variance of the conjugated variable, such that Heisenberg uncertainty principle still holds true for the product of the variances. Squeezed states have kept their promise in improving measurement accuracy beyond quantum noise [1,2] and have become central in quantum optics [3] through demonstrated applications in quantum metrology and sensing [4,5]. Advanced techniques to generate squeezed light [6][7][8] facilitated the detection of gravitational wave [9][10][11]. Theoretical works have proposed applications in quantum information, where coupling a qubit to a squeezed reservoir allows erasing information below the Landauer's limit [12]. In the context of quantum thermodynamics, the proposed theories of coupling the working medium of a Aurélia Chenu : aurelia.chenu@uni.lu nanoscale heat engine to a squeezed reservoir to generate work beyond the Carnot's limit [13][14][15][16][17][18] have been experimentally demonstrated using a vibrating nano-beam driven by squeezed electronic noise [19].
Renewed interest in squeezing has come with progress in quantum optomechanics [20][21][22][23][24]. In a simple parametric interaction where the spring constant of a mechanical oscillator is controlled with radiation pressure forces, the emergence of mechanical instabilities prevents reducing fluctuations to at most 50% (the 3 dB limit) in the steady state [25,26]. Schemes to generate squeezing below this limit have been developed by combining parametric driving and weak measurements [26,27]. Continuous measurements, where measuring one variable-e.g. position-precisely reduces its associated variance, present an intuitive technique that has been broadly explored [28][29][30][31][32][33][34][35][36][37]. Alternatively, a simple but powerful scheme that lifts the requirement of explicit measurement or feedback has been put forward using a dissipative mechanism where the driven cavity acts Figure 1: Schematic representation of the control processes considered in this work. Starting from an initial thermal state with isotropic density in phase-space (top), we design dynamical protocols to generate a squeezed thermal state (bottom) at controlled temperature in arbitrary time using two different experimental implementations: (left) in a harmonic oscillator with controlled frequency, or (right) using two-photon Raman interaction. Thermalization is achieved by engineering a dissipator in position.
as an engineering reservoir [38]. This theoretical protocol shows the similarities between processes relying on coherent feedback or reservoir engineering [14,39], and has been experimentally demonstrated in [40,41].
While much progress has been achieved for increasing the squeezing parameter, the protocols have so far been restricted to unitary dynamics and do not allow for the control of entropy. We here lift this limitation and provide protocols to generate squeezed thermal states at controlled temperature in a fixed time. Controlling the temperature of squeezed thermal states is all the more relevant since the variance of such states not only depends on the squeezing amplitude but also on the average thermal phonon number [42][43][44], as further discussed below.
Specifically, we focus on two different methods, already known as useful to generate squeezing, that we extend to open setups: (i) squeezing from nonadiabatic driving of the oscillator frequency-how is the trap control frequency modified by the dissipative dynamics; (ii) squeezing through the use of twophoton Raman interaction. In both cases, knowledge of the analytical dynamics allows finding the control processes through reverse engineering. Thermalization is achieved through the use of white noise with controlled amplitude, that generates the open dynamics.
The presented protocols aim at the dynamical control of states without relying on adiabatic evolution. In this sense, they fall under the umbrella of Shortcuts to Adiabaticity (STA) [45,46]. In essence, STA provide the control Hamiltonian to generate in a fixed time the state otherwise reached through a reference adiabatic trajectory. Extension of STA to open quantum systems requires, in addition to the control Hamiltonian, a control dissipator. This was first explored for Markovian dynamics [47] and a general scheme has recently been put forward for arbitrary dynamics [48]. Other works have shown how to control the thermalization of a harmonic oscillator [49][50][51]. Here, we provide protocols combining squeezing and thermalization. Note that squeezed thermal states have been experimentally achieved in, e.g., a massive mechanical harmonic oscillator using sudden quenches [20]. Our protocols allow to operate in the quantum regime and avoid creation of excitations thanks to the use STA techniques. Also, they allow us to control the state squeezing parameters as well as its temperature.
The paper is organized as follow: Section 2 presents the general evolution for a squeezed thermal state with time-dependent parameters. Section 3 focuses on the harmonic oscillator setup, clarifying the relations of known STA with squeezing. Section 4 presents the results using two-photon Raman interaction that allows controlling the squeezing amplitude, phase, and state temperature.

Squeezed thermal states
We consider the thermal state σ 0 = e −β 0 H 0 /Tr(e −β 0 H 0 ) of a harmonic oscillator (HO) with Hamiltonian H 0 = ω 0 (a † 0 a 0 + 1/2). This initial thermal state is driven to a target squeezed thermal state σ f at an arbitrary final time t f following where we allow for changes in entropy through the dimensionless time-dependent parameter ε t = ω t β t and its initial value ε 0 = ω 0 β 0 . The partition function Z t normalizes the state. The squeezing operator is defined from the squeezing parameter z t = rt 2 e −iφt with time-dependent amplitude r t and phase φ t , which defines the correlations between position and momentum. The annihilation operator a 0 is defined from The static properties of state (1) have been thoroughly studied and described in e.g. [52]. We focus on its dynamics to design control protocols generating a squeezed thermal state at arbitrary final temperature β −1 f and target parameters {r f , φ f } at the end of the control protocol, t f . Let us first rewrite (1) as the thermal state of the squeezed harmonic oscillator. Using S r,φ S † where H gho = ωt ω 0 S r,φ H 0 S † r,φ and the partition explicitly reads Z t = Tr(e −βtHgho ). We split define new creation and annihilation operators, explicitly found using the Baker-Campbell-Hausdorf (BCH) formula [53] , yielding These 'A' operators are bosonic operators fulfilling For each A † t boson created, there is both creation and annihilation of some 'a' bosons 1 In this basis, the squeezed harmonic oscillator simply reads We proceed to describe the evolution of its thermal state. Under unitary evolution, the state dynamicsσ t = known as the counter-diabatic (CD) Hamiltonian [56][57][58], which ensures that each eigenstate remains instantaneous eigenstate of the time-dependent Hamiltonian and evolves as i |ṅ t = H cd |n t . The explicit form of this Hamiltonian is obtained from the kth time-derivatives of (e −iφt a 2 0 −e iφt a †2 0 )-see App. A-yielding Let us comment on the possible implementations set by the phase. Case φ t = 0: (Section 3) Squeezing with no final correlation betweenx andp, i.e. using S r,0 , reduces the Hamiltonian (7) to and the evolution readsσ We elaborate in Sec. 3.1 how this Hamiltonian can be implemented in a harmonic oscillator by varying the trap frequency and relates to common STA techniques. Indeed, since −i (a 2 0 − a †2 0 ) = {x,p}, the counterdiabatic term recovers the known squeezing term [56] such that, for a proper choice of r t , it can be recast into the form that, in first quantization, reads H Case φ t 0: (Section 4) The generalized HO 2 ) commutes with the dynamical state 1 Squeezing then appears similar to the physical setup of the independent-boson model [54], which is best dealt with using two different basis for the bosons [55]. and can be removed from the Hamiltonian. The remaining Hamiltonian in Eq. (7) can be implemented in a rotated frame, as we detail later in Sec.

For now, note that the unitary rotation
This unitary evolution is dictated by the two-photon Raman Hamiltonian. We further detail how this allows implementing (7) provided thatφ t = −2ω 0 . Importantly, the phase linearly depends on the process time. So with this setup, a process of fixed time t f generates a fixed squeezing phase, φ f = −2ω 0 t f . We show in Sec. 4 how to lift this constraint and generate a squeezed state with arbitrary position-momemtum correlation in arbitrary time.

Fast Squeezing and Thermalization through trap and dephasing control
On one hand, it is quite established that squeezing can be achieved with a change of the trap frequency, as proposed in trapped ions already decades ago [59][60][61] and demonstrated experimentally [62]. On the other hand, creating a thermal state from the thermal state of a different system requires rearrangement of the initial distribution of eigenstates so as to match the Gibbs distribution of the final system. STA protocols generate, in a finite time, the adiabatic evolution of a reference Hamiltonian [46], thus preserving the initial eigenvalue distribution. Although these two results are well established, the connection between squeezing and STA on a HO seems to not always be made. We explicit it here and then consider dissipative dynamics to extend the technique to generate squeezed thermal state at arbitrary temperature. To do so, this Section focuses on generating squeezed thermal state with no correlation betweenx andp, i.e. φ = 0.

Squeezing through trap control
We first consider unitary, isentropic dynamics, i.e. ε t = ε 0 constant. In the case φ t = 0, the squeezed thermal state directly maps to the instantaneous thermal state of a HO with time-dependent frequency ω t , . (10) Indeed, for φ t = 0, the A t operator (5) becomes a t = S r,0 a 0 S † r,0 = a 0 cosh r t + a † 0 sinh r t . It evolves aṡ a t =ṙ t a † t which maps, for r t = ln ω t /ω 0 , to the annihilation operator factorizing the time-dependent HO, H . Note that the operation S r,0 is also known as a dilatation [63], T w = exp − i log(w) 2 (xp +px) with w ≡ ω 0 /ω t , that transforms position and momentum as S r,0 f (x)S † r,0 = f (x/w) and S r,0 f (p)S † r,0 = f (wp), respectively. The time-dependent HO itself is thus equivalently a squeezed or dilated HO The dynamics directly follows from (8) For the purpose of experimental implementation, let us consider the state t = U Ω 0 σ t U † Ω 0 in a frame rotated by the unitary where Ω 0 so far is an arbitrary, time-dependent frequency. The evolution of this density matrix˙ t = − i [H c , t ] is governed by Taking Ω 0 = −ω t /(2ω t ) = −ṙ t =ẇ/w removes the correlations in position and momentum. Under this condition, the control Hamiltonian H c = p 2 2m + 1 2 mω 2 cx 2 is a HO with time-dependent control frequency (13) 3 We used So for a given reference trajectory with frequency ω t , ω c is the frequency to be implemented to do the shortcut. This control frequency can be recovered inserting w in the Ermakov equation [64],ẅ + ω 2 eff w = ω 2 0 /w 3 , and corresponds to known results from local counterdiabatic driving [57,65]. The state implemented in the lab evolves as and corresponds to the target squeezed state at the end of the process only-at which time U Ω 0 becomes the identity and f = σ f (φ = 0).

Extending the range of accessible squeezed states with a control dissipator
The thermal state (10) is diagonal in the instantaneous Fock state basis and reads σ t = is conserved during any unitary evolution. The product βω, that characterizes the eigenvalue distribution of the thermal state, is constant under unitary evolution. We refer to 'cooling' for processes decreasing the von Neumann entropy [66], which contrast with phase-space preserving processes-with constant entropy. To extend the range of accessible states, we allow for changes in temperature and entropy during the dynamics, taking ε t time dependent.
Another motivation to design open protocols for squeezing is to enhance the variance. Indeed, the variance in position of a squeezed thermal state depends on both the average phonon numbern = 1 e ε t −1 and the squeezing parameter r t [42]. Specifically, using the operatorx = mω 0 (a 0 + a † 0 ) on the state σ t (φ = 0) = (1 +n) −1 ∞ n=0 n n+1 n |n t n t | has a variance in position Using the link to the time-dependent harmonic os- tails. So open protocols allow controlling the variance with two parameters, (ω t , ε t ). Figure 2 shows the variance ratio between final and initial states as function of the inverse temperature and trap frequency relative changes. While the variances at reach with unitary protocols are limited to those on isentropic lines, the map is extended to arbitrary values (including beyond the 3dB limit) thanks to changes in entropy. Unitary dynamics are along these isentropic lines and restrict the target state to β f ω f = β 0 ω 0 , thus also restricting the accessible variances for given initial conditions. By contrast, open dynamics and engineered dissipation extend the variances at reach to the full map. Blue background correspond to cooling processes, orange is for heating.

Let us now detail the open control protocols.
When entropy is allowed to change-ε t is time- with H (0) cd given in Eq. (8). The system is initialized in the thermal state (10) and evolves according to (15). As in the unitary case, for the sake of experimental implementation, we consider the state matrix t = U Ω σ t U † Ω rotated by the unitary where all terms accounted for population changes are in the counter-diabiatic dissipator, D cd ( t ) = nṗ n,t U Ω |n t n t |U † Ω . Alternatively, part of the counter-diabatic Hamiltonian can be written as a 'control' harmonic oscillator, H c , with a 'control' frequency chosen of the form of the closed results (13) if the dynamics were unitary (cf. Eq. 12)-and an additional frequency Ω 1 that accounts for changes due to the open dynamics. The master equation (16) thus becomeṡ The 'control' dissipator can be written in a compact form by defining the annihilation so the position operator x is equivalently written in one basis or the other.
This dissipator can be further written in a more 'experimentally-friendly' form. For this, notice that the relations the dissipator (19) can be recast as The control dissipator thus becomes the well-known form of localization in the position eigenbasis, often referred to as Joos-Zeh term [68,69], and is easily implementable in current experimental platforms. In turn, the modulation of γ t can be engineered, e.g., by postselection measurement of the position or via stochastic parametric driving, as proposed in [51] and also used below in Section 4.3.
The designed master equation of interest thus reads  The Wigner function is plotted at times t = 0, t f /2, and t f . It starts as a symmetric Gaussian and rotates in phase-space during the process to reach a squeezed state along thex quadrature, as expected for compression. The inset shows the dissipation rate γ t .
with the control parameters It allows generating the target state f = σ f at final time-since thenω f =ε f = 0. Implementation easily follows from knowledge of the control parameters (23): given the boundary conditions β 0 , β f and ω 0 , ω f , one can fixe the time evolution as e.g. a fifthorder polynomial Ansatz, p(τ ) = 10τ 3 −15τ 4 + 6τ 5 , on β t and ω t .
With this choice, we illustrate the dynamics in Figure 3, which shows the control parameters and the state Wigner function along the open dynamics. As expected for a compression process (ω f = 3ω 0 ), the Wigner representation of the final state evidences a thermal state squeezed in position.
To summarize this part, we showed how opening the dynamics influences the control frequencycompare Eq. (23a) with (13). The particular choice of the correction Ω 1 made in Eq. (20) leads to a dissipator controled in position, as proposed in [48,51]. Should this dissipator not be the most experimentally suited, Eq. (19) provides the more general result for the control dissipator. Also note that the dynamics is resolved at the level of operators and not restricted to the coordinate representation. As such, the results here are more general and complements these previous works.

Squeezing and Thermalization with two-photon Raman interaction
A squeezing protocol alternative to controlling the trap frequency is based on two-photon Raman interaction, that was successfully used to squeeze the ground vibrational of trapped ions [70][71][72]. We detail below the experimental setup used for implementation of equation (9), and extend the known technique to allow for (i) squeezing in arbitrary time thanks to reverse engineering, and (ii) at arbitrary temperature with engineered dephasing. We notably explain how the modulation of the lasers amplitudes allow to modify the variance (14) in arbitrary time.

Experimental setup
Consider a trapped ion interacting with two monochromatic laser beams-see Fig. 4. In the experimental situation of interest, the electronic structure of the ion is reduced to a two-level system described by the atomic Hamiltonian H a = ω 2 σ z , with σ z = |e e| − |g g|. The motion of the trapped atom can be considered harmonic in all three dimensions, as obtained either from a classical or quantummechanical treatment [73,74], and thus described by H m = ω 0 (a † a + 1/2). With suitable electromagnetic fields, the electronic levels can be coupled to each other and to the vibrational motional degrees. Each of the electromagnetic field is treated as a classical plane wave of the form, in the direction x of interest, E l (x, t) · x = A l (t)(e i(k lx −ω l t−Φ l ) + c.c.)/2 with time-dependent amplitude A l (t), wave vector k l = k l x, and detuning δ l from the atomic transition, ω l − δ l = ω. The interaction Hamiltonian resulting from the applied two laser fields can be described as [74] H int (t) = l={1, 2} 2 with σ x = |g e| + |e g|. The Rabi frequency describing dipole coupling to a single charge q is given by Ω l /2 = q g|x|e A l (t).
We aim at preparing a squeezed thermal state on the vibrational levels of the system with total Hamil- tonian starting from an initial vibrational state that is thermal.

Closed dynamics
We first consider the unitary dynamics and denote |ψ t the solution of the Schrödinger equation. It is useful to change the energy scale [75] and look at the evolution of |Ψ t = U r,t |ψ t rotated by a unitary transformation U r,t ≡ e i Hrt . The rescaling Hamiltonian H r = H a + H m + ∆ 2 σ z effectively shifts the electronic energy of the gap ω into an energy defined by the average laser detuning, ∆ = (δ 1 + δ 2 )/2, and yields to an interaction picture. The rotated state evolves as i |Ψ t = H|Ψ t , with H ≡ U r,t H tot U † r,t + i U r,t U † r,t explicitly reading The position being quantized, we usedx = x 0 (a + a † ) to write e ik lx = e iη l (a+a † ) in the expression of the electromagnetic field [72]. The interaction picture leads to using the time-dependent operators a t ≡ U r,t aU r,t = ae −iω 0 t andx t ≡ U r,tx U r,t = x 0 (a t + a † t ). The Lamb-Dicke parameter η l = k l x 0 is defined from the extension of the ground-state wave function of the reference oscillator, x 0 = /(2mω 0 ). The evolving wave function is a superposition of the electronic ground and excited states dressed with the vibrational levels |n , and we look for a solution in the form |Ψ t = n g n (t)|g, n + e n (t)|e, n . The electronic and vibrational degrees of freedom can be decoupled through an adiabatic elimination [75], which assumes constant excited-state population. We follow [56,76] and set ω 1 − ω 2 = 2ω 0 . Keeping only the resonant, second blue sideband, which effectively is a vibrational form of the RWA, and neglecting the Lamb term shifting, the evolution of the atom state density is dictated by the effective squeezing Hamiltonian as computed in App. E (27) One can rewrite the evolution of |Ψ t in the Liouville-von Neumann form, so that ρ t = |Ψ t Ψ t | evolves as with and Φ t = (Φ 1 − Φ 2 ). Choosing the dephasing between the lasers to be Φ 1 − Φ 2 = π 2 recovers the 2-photon Raman Hamiltonian.
Importantly, this recovers the dynamics of the squeezed state given by (9) in the rotated frame, thus generating the general squeezing Hamiltonian (7) in the original frame. In this setup, the squeezing parameter is directly related to the experimental parameters asṙ This simple relation is experimentally very relevant. It implies that using a constant amplitude for the lasers leads to a squeezing parameter linear in time. So to achieve a fixed squeezing parameter, one must wait a given time. By contrast, the same squeezing parameter can be achieved in an arbitrary time through reverse engineering. For example, consider a fifth-order polynomial interpolating between the initial and final squeezing parameter, r t = r 0 + (r f − r 0 )p(t/t f ). Its derivative gives, through Eq. (30), the laser amplitudes needed to reach the target squeezing in arbitrary time. This is further illustrated in Fig.  5, that shows the benefit of reserve engineering the laser amplitudes over using a constant amplitude in order to reach a target squeezing parameter in a desired time. Furthermore, as mentioned in Sec. 2, the time of the process determines the phase. It is then possible to choose the time of the process accordingly to the desired target phase, t f = − φ f 2ω 0 . We now extend the dynamics to open processes and solve the state evolution. In addition to provide control over the temperature, the derived solution allows generating an arbitrary squeezing phase.

Open dynamics: Engineering the master equation in a stochastically-shaken trapped ion
In order to generate a controlled dissipation and extend the map of final states at reach, we process as in the HO, and use stochastic processes. Specifically, we add to the total Hamiltonian a stochastic component, and consider This system is now characterized by the Wiener process W t = W 0 + t 0 ξ t dt defined in terms of the normally distributed random variable ξ t taken as white noise-with zero mean and vanishing correlation, ξ t ξ t = δ(t − t ). The stochastic term allows to create the control dissipator and has the advantage of being readily implementable via continuous quantum measurement or in a stochastically-shaken trap [77]. Note that the stochastic term in (31) is taken as acting only on the electronic ground state. While this yields a rigorous analytical derivation, it can be experimentally challenging, so an alternative scheme is presented in App. G that, instead of shaking the trap, relies on two additional laser beams, one having a stochastic amplitude.
Let |ψ t denote the solution of the Schrödinger equation, |ψ t+dt = e − i hstdt |ψ t . Following the uni-tary results, we look at the evolution of the state vector |Ψ t = U r,t |ψ t . The influence of the additional stochastic term is detailed in App. F. The density matrix of interest is the ensemble one, obtained from averaging over the realizations of the noise and denoted ρ t = |Ψ t Ψ t | . We find its evolution dictated by the master equation the second line following from the RWA. It corresponds to a master equation of Lindblad form, where the dissipators, defined from D(a) = aρ t a † − 1 2 {a † a, ρ t }, are modulated with an amplitude κ t = γ t x 2 0 . The parameter α t is the same as in the unitary case-Eq. (29). We solve the dynamics and show how this setup can be used to generate a target squeezed thermal state.

Solving the dynamics
The system is initialized in a state with density matrix |g g| ⊗ σ 0 , where σ t ≡ e −βtHm /Z t denotes the thermal state on the vibrational manifold at initial inverse temperature β 0 . We next solve the dynamics to find the dynamical control parameters {α t , κ t } for which the squeezed thermal state is solution of the dynamics (32). The time-dependent parameter λ t ≡ −β t ω 0 allows varying the temperature. It is useful to work with the factorized form, that we derive in normal ordering following McCoy [78] as (see App. D) where the parameters are defined as J t = j(r t , λ t )e iφt with the real functions j(r, λ) ≡ 1 2 sinh(2r)(e 2λ −1) 2(cosh 2 r−sinh 2 re 2λ ) , B t ≡ − ln 1 + (e λ t −1)(cosh 2 (rt)+sinh 2 (rt)e λ t ) cosh 2 (rt)−sinh 2 (rt)e 2λ t , and the normalizing constant K t -given explicitly in Eq.
(S24). The master equation (32) gives Using the adjoint representation, detailed in App. F, yields to the simple system linking the vector of the control parameters v c = (κ α R α I ) T to the vector of the squeezing parameters v sq = (J RJIḂ ) T with T the transposition. Specifically, we obtain v c = M −1 t v sq (36) with the transfer matrix This is the main result of this section. Eqs. (36 -37) give the control parameters α t = |α t |e iΦt = α R + iα I and κ t = γ t x 2 0 to engineer the squeezed state characterized by J t = J R + iJ I and B t at the desired temperature through λ t ≡ −β t ω 0 .
We show numerical applications of the controlled parameters to be implemented to drive an initial (possibly squeezed) thermal state characterized by The state parameters J R , J I and B t are assumed to follow a smooth evolution taken as a fifth-order polynomial, with additional boundary conditions taken as null first and second derivatives at initial and final times. The relative detuning between the lasers is fixed to ω 2 − ω 1 = 2ω 0 , as required to generate the squeezing Hamiltonian (27). The control parameters are obtained by solving Eqs. (36,37). The dynamics can thus be implemented through the controlled dephasing strength κ t = γ t x 2 0 , the controlled laser amplitudes, and their Rabi frequencies. The latter are directly related to the control parameters α = |α|e i(Φ 1 −Φ 2 ) that gives the relative laser phases Φ 1 − Φ 2 = arctan (α I /α R ) and Rabi frequencies through |α| = (η 2 − η 1 ) 2 Ω 1 Ω 2 4∆ . Figure 6 shows the control parameters for squeezing with different temperature conditions, namely cooling, isothermal, and heating. The normalized laser amplitude appears to be quite similar for all squeezing processes, which can be expected as it mainly controls the squeezing amplitude. In turn, its maximum is influenced by the variation of squeezing, as shown in Fig. 7. When the state retains an isentropic density (∆r = 0), no squeezing term is needed, as intuitively expected. Figures 7 and 8 show the influence of changing the squeezing amplitude and temperature, respectively. We verify that the dephasing strength is 'symmetric in squeezing' in the sense that squeezing by a positive or negative variation |∆r = r f − r i | only changes the sign of the dephasing, not its strength.
As mentioned above, implementation of the stochastic Hamiltonian (31) assumes a spindependent term on the position of the trap, that could be developed following the techniques proposed in e.g. [79]. This allowed for a rigorous derivation of the effective Hamiltonian through the adiabatic elimination. Shaking the full trap (ground and excited electronic states) would require further approximations of the excited state populations, although the adiabatic elimination might still hold at large detunings. Further work could be done using the recently developed adiabatic elimination for open bipartite systems [80][81][82]. An experimental  alternative is to install a feedback loop that enforces the qubit to remain in its ground state [83]. Should the proposed model still be experimentally limiting, we provide in App. G an alternative scheme where the dissipator is engineered with two additional laser field instead of shaking the trap.
Finally, note that we have here focused on the trapped-ion setup for the sake of proposing a scheme that can be directly implemented experimentally. However, the dissipator need not be of the form [x, [x, ρ]]. For example, using a dissipator of the type D sq = γ (n + 1) L(a)ρ+γnL(a † )ρ in Eq. (32), with L(a) = aρa † − 1 2 a † aρ + ρa † a , could also be used to generate a squeezed thermal state in a photonic platform, where the squeezing Hamiltonian could be obtained with, e.g., parametric downconversion.

Conclusion
Starting from the general evolution for a squeezed thermal, we first clarified how squeezing without phase control can be achieved in arbitrary time by modulating the trap-frequency of a harmonic oscillator and as such, relates to known STA techniques. In turn, control of the phase can be implemented with a two-photon Raman Hamiltonian. Importantly, the two approaches presented here include dissipative dynamics in order to control the state entropy, that is engineered using stochastic fields. We provided a detailed analysis in a trapped-ion setup, giving the control laser amplitude, relative phase, and dephasing strength suited to generate a target squeezed state in arbitrary time. The general formalism could also capture, e.g., photonic thermal states squeezed by parametric downconversion in a lossy cavity [84] and is thus adaptable to other experimental platforms. Among possible applications, the generated squeezed states can be used for trapped-ion transport [85], which is relevant to quantum computing architectures.

Acknowledgements
It is a pleasure to thank Adolfo del Campo, Kihwan Kim, and Mauro Paternostro for insightful discussions and comments on the manuscript.

A Finding an expression for S r,φṠ † r,φ
To find an explicit expression of S r,φṠ † r,φ for any squeezing angle φ t , it is useful to define the operatorã t = e −i φ t 2 a 0 . It fulfills the bosonic commutator relation [ã t ,ã † t ], and gives S r,φ = e r t 2 (ã 2 t −ã †2 t ) . We can expand the exponential in Taylor series and look for S r,φṠ † r,φ in terms of the kth derivatives ofã 2 t −ã †2 t . From the time derivativeȧ t = −iφ t 2ã t , it follows that tã t ] to obtain the expression for k = 1 in the form of commutators, explicitly, This form allows to generalize the results and obtain Eventually, we get which corresponds to Eq. (7) given in the main text.

B Wigner function for squeezed thermal state
We use the coordinate representation of the state evolving in a time-dependent harmonic oscillator under dephasing in position [51] y|ρ t |x = with k t = mωt and ε t = β t ω t . This readily gives the Wigner function as Integrating over momenta, one can find the marginal distribution for the position and identify the variance in position ∆x as This result recovers the known variance for a squeezed thermal state [42].

C Control dissipator D c (19)
The 'control' dissipator is defined using Eqs. (16)(17)(18) as In order to find a compact form, it is useful to define the operator (S10) First note that since a t +a † t = b t +b † t , the position operator is equally represented in both operator basis, namelŷ we can recast the control dissipator (S9) into the compact form given in Eq. (19) of the main text.

D Factorization of the squeezed thermal state
It is useful to write the squeezed thermal state in a product form in order to solve its dynamics. We show how to obtain (S11) The full analytical demonstration we propose here is alternative to the one provided in [86].
The factorized form of a function can be obtained following the use of differential equations, as first proposed by McCoy [78]. For any function g(a † , a) of the non-commuting operators [a, a † ] = c, the partial derivative can be defined as [87] c ∂g ∂a † = [a, g], and c Note that c is a constant that will be taken equal to unity at the end, but is useful in the general derivation for the purpose of normalization. The expression of the function g in normal ordering (annihilation operators a to the right, creation a † to the left) is obtained by integration of a system of partial derivatives. We consider the particular function ρ = e λA † A = n λ n n! (A † A) n , which is quadratic in a and a † since, for (S13) To obtain differential equations, we start from the obvious observation that Using the relations between the 'A' and 'a' operators (5) this readily gives With the future integration in mind, we write the aρ term on the r.h.s as ρa + ∂ρ ∂a † and obtain the first differential equation A similar equation can be obtained starting from the observation that A(A † A) n = (AA † ) n A. This gives [A, ρ] = (e cλ − 1)ρA and yields to the differential equation where again, we have chosen to have terms on the r.h.s in the ordering ρa and a † ρ. So we now have the system of differential equation with the constants (S19) It is now easy to verify that the factorized form (S11) is solution of the system (S18). So we obtain the following factorized form, in normal ordering For the squeezed thermal state, f + = cosh(r t ) and f − = sinh(r t )e iφt as follows from (5), which allows to relate directly the squeezing parameters with the factorized form as given in (34). In order to compute the constant K, we further follow the derivation proposed by McCoy [78]. During the factorization, only the commutation relation is important. So we choose to replace a † →x and a → cp = c d dx , with [p,x] = c. We look at how λA † A acts on 1, and denote this action λA † A{1}. With the change of The constant term in (λA † A) n can be found in applying n times the operator (λA † A) on the identity, and is of the form a n c n . Let us denote the constant in the first term of the serie P 1 (c) such that P 1 (c) = ∞ n=0 a n c n /n!. We denote P 2 (c) the constant when acting twice (λA † A), namely the constant term in e λA † A {λ |f ∂c . We further know from Eq. (S11) that . At the limit for which the operators commute, K(c → 0) tends to unity. Hence, One can compute the integral Inserting this in (S22) and using the constant c = 1, we obtain the normalization constant in the factorized state (34, S11, S21) as (S24)

E Effective Hamiltonian in the Unitary case
We are interested in the evolution of the state |Ψ t characterized by the Hamiltonian Let us first give the explicit form of the interaction Hamiltonian in the rotated frame. The atomic part evolves as e i 2 (ω+∆)σz σ x e − i 2 (ω+∆)σz = e −i(ω+∆)t |g e| + h.c.. The bosonic part is obtained from e iω 0 ta † a a † = a † e iω 0 t(a † a+1) that gives e iω 0 ta † a e iη l (a † +a) e −iω 0 ta † a = e iη l (a † e iω 0 t +ae −iω 0 t ) . Keeping only the terms with the lowest frequency (RWA), we thus have We are looking for a solution of the wave function as a linear combination of the dressed basis |Ψ t = n (e n (t)|e, n + g n (t)|g, n ) .
For large detuning, |∆| |Ω l |, ω 0 , a state initially in the electronic ground state mainly remains in this electronic level. The small population of the electronic excited state can be eliminated abiabatically. We thus seṫ e n (t) = 0, and the evolution follows as with the effective Hamiltonian equal to Eq. (S31). After the adiabatic elimination, the increment reads It is now easy to characterize the evolution of the density matrix ρ st = |Ψ t Ψ t |. The Leibnitz chain rule that, in the Itô calculus, generalizes to d(AB) = (A + dA)(B + dB) − AB = (dA)B + A(dB) + dAdB [88], yields which preserves the norm at the level of each individual stochastic realization. The density matrix of interest here in the ensemble one, obtained from averaging over the realizations of the noise and denoted ρ t = ρ st . Since the average of any function F t of the stochastic process vanishes, F t dW t = 0 [88], we find that the evolution of the ensemble density matrix ρ t is dictated by the master equation (32) given in the main text. We now solve this equation and find the control parameters for which the squeezed thermal state (33) is a solution. To do so, we look at dρ dt ρ −1 and use the factorized form of the squeezed thermal state that allows recasting all needed terms of the master equation (35) We then explicit the transformation for each element of the basis. The action of the adjoint is thus a linear transformation that can be represented in matrix form in the basis B, the needed terms being explicitly Finally, dρ dt ρ −1 = Figure 9: Experimental setup: 2-photon Raman interaction is generated by the (blue) laser pair with ω 2 − ω 3 = 2ν, while dephasing is generated with the (red) laser pair, ω 1 − ω 0 = ν, one amplitude being taken as stochastic.
at the evolution of the rotated vector |Ψ t ≡ U r,t |Ψ t . The unitary U r,t ≡ e i Hrt is defined from the rotation Hamiltonian H r = H a + H m + ∆ So this dynamics creates the squeezed thermal state (33) provided that the control parameters fulfill  Figure 10 presents the control parameters for implementation of the dynamics for cooling, isothermal and heating processes. Interestingly, in the case of simple cooling and heating (with no squeezing), the squeezing hamiltonian is not zero anymore, which is different from the former setup (cf. Fig. 6). Adding squeezing (dashed curves) leads to similar results. In turn, the parameter controlling the dephasing, κ t , is positive for heating and a negative for cooling, which matches with intuition. The influence of temperature and squeezing variations on the maxima of control parameters are presented in Figures 11 and 12.   : Maxima of the control parameters κ max and |α| max as function of the initial or final squeezing amplitudes, for different variation ∆r = r f − r i . The two control parameters are 'symmetric in squeezing', i.e. their maxima only depend on the absolute value |∆r| of squeezing variation. In other words, a unique value of |κ max | is associated to a given couple of values (r i , r f ). Note that, while a high variation of the squeezing parameter is hard to engineer, only small values are needed since the variance exponentially depends on the squeezing amplitude-Eq. (14). For instance ∆r = 2 drastically reduces the variance by seven times. Plots are for t f = 1.