Non-equilibrium quantum thermodynamics of a particle trapped in a controllable time-varying potential

Many advanced quantum techniques feature non-Gaussian dynamics, and the ability to manipulate the system in that domain is the next-stage in many experiments. One example of meaningful non-Gaussian dynamics is that of a double-well potential. Here we study the dynamics of a levitated nanoparticle undergoing the transition from an harmonic potential to a double-well in a realistic setting, subjecting to both thermalisation and localisation. We characterise the dynamics of the nanoparticle from a thermodynamic point-of-view, investigating the dynamics with the Wehrl entropy production and its rates. Furthermore, we investigate coupling regimes where the the quantum effect and thermal effect are of the same magnitude, and look at suitable squeezing of the initial state that provides the maximum coherence. The effects and the competitions of the unitary and the dissipative parts onto the system are demonstrated. We quantify the requirements to relate our results to a bonafide experiment with the presence of the environment, and discuss the experimental interpretations of our results in the end.


I. INTRODUCTION
The potential for enhanced performances above and beyond the possibilities offered by classical devices underpins the development of quantum technologies for applications, in information and communication technology, sensing, and computation [1][2][3][4][5][6]. Recently, it has been realized that quantum laws could also be exploited to manage more efficiently the energetics of nanoscale technologies, thus contributing to the emergence of a framework for the thermodynamics of quantum processes [7][8][9][10][11][12][13]. However, such advantages do not come easily due to unavoidable necessity to protect the mechanisms responsible for the quantum process being studied from detrimental environmental effects [14][15][16][17]. Such deleterious influences have a tendency to scale with the size of the system being used, and become -in principleprohibitively severe in the mesoscopic and macroscopic regimes.
Fortunately, impressive progress on the experimental control of quantum systems and dynamics have been achieved in the course of the past twenty years. This has allowed the design and implementation of effective strategies for the control of mesoscopic systems such as cold atoms [18], large arrays of superconducting systems [19,20], and (electro-/opto-/magneto-)mechanical structures [21,22]. Levito-dynamics [22], i.e., the levitation of nano-and micro-objects in vacuum, holds the potential to become a key experimental platform for the demonstration of quantum features at the mesoscopic scale. It enables dynamical control over nearly arbitrary poten- * Email: qwu03@qub.ac.uk tials as well as high degrees of isolation from environmental influences. Today, experiments based on opticaltweezer technology [23] have allowed near Heisenberglimited position readout, real-time control of motion and preparation in the ground state of motion [24][25][26]. Also, dynamical shaping of an optical potential for a levitated nanoparticle has been demonstrated, with the scope to implement a logically irreversible transformation [27]. These developments pave the way to the exploration of non-equilibrium phenomena in open mesoscopic systems and thus the consolidation of a framework of controllable quantum thermodynamics of large-scale systems [28].
In this work, we make significant theoretical steps towards the characterization of such framework by addressing thermodynamic irreversibility stemming from the non-equilibrium process involving a levitated nanoparticle -subjected to both thermalisation and decoherence mechanisms -in a potential landscape that is changing from a quadratic to a double-well configuration, similar to the experimental situation in [27]. We focus our attention on the quantification of irreversible entropy production [29], a key quantity that characterizes logical and thermodynamic irreversibility, constraints the performances of quantum and classical engines, and appears to be strongly related to the occurrence of nonequilibrium quantum critical phenomena [29,30]. Our investigation sets the methodological toolbox for the successful simulation of non-equilibrium processes subjected to real-time potential-shaping transformation, as envisioned in particular for levitated systems. It also establishes the context for the quantification of potential thermodynamics-based limitations to the efficiency of quantum memories.
The remainder of this paper is organized as follows. In Sec. II we introduce the model and introduce a suit-able angular-momentum algebra that allows for an agile simulation of the dynamics of the system. We illustrate such capabilities through significant numerical examples. In Sec. III, we introduce the thermodynamic quantities of interest, namely the Wehrl entropy, its production rate and the entropy flux rate [29]. These form a toolbox of figures of merit for the characterization of thermodynamic irreversibility stemming from the nonequilibrium process undergone by the system, and the energy-exchanging mechanisms with the environment resulting from its open-system dynamics. In Sec. IV, we focus on a quantitative characterization of such quantities for a system undergoing temporal shaping of its trapping potential from harmonic to double-well shape. An experimental perspective is given in Sec. V, while we draw our conclusions in Sec. VI.

II. THE MODEL
We consider a mechanical system under the influences of the surrounding environment. The dynamics of the system can be accounted by the following master equation for the density matrix ρ, The second term in Eq. (1) is the thermalisation dissipator D th [ρ], which describes the effect of heat exchange between system and environment and would eventually lead to thermal equilibrium at the environmental temperature. Such term takes the form of the Lindblad superoperator [7,11] Here, a and a † are bosonic annihilation and creation operators of a harmonic oscillator of frequency ω that will work as a basis of our problem, γ is the coupling strength between system and its thermal environment, which has a mean number of excitationsn = (e β ω − 1) −1 with β = 1/k B T being the inverse temperature of the environment. The last term in Eq. (1) quantifies the decoherence effects of the collision of the residual gas in the vacuum, and it can be described through a localisation term [17,31,32], which reads where Λ is the coupling strength and x = /2mω(a † + a) is the position operator. Finally, H s (t) in Eq. (1) denotes the system Hamiltonian, whose time dependence can be used to dynamically modify the potential. In particular, we focus on the following form where we have introduced the momentum operator p = i mω/2(a † −a). Eq. (4) describes a harmonic potential modified by an additional time-dependent term whose form we take as Here, E is a suitably chosen energy scale. The protocol embodied by the time-variation of H add (t) creates a double-well potential from one of the form of an inverted Gaussian via a suitable change of α(t) andᾱ(t).
Here, we define a protocol that linearly switches the potential in time from 0 to τ , such thatᾱ(t) = 1 − α(t), and α(t) = 1 − t/τ , with τ the characteristic time of the protocol, and we fix E = 10 ω so that the ground state energy is below the central peak of the final double-well potential. The initial potential (at t = 0) well approximates the low energy-sector of the harmonic oscillator. The total potential is shown in Fig. 1a and we refer to Sec. V for experimental considerations on the plausiblity of the choices made here. In order to analyze the dynamics of the system as the shape of the trapping potential varies in time, we resort to a spin-coherent state picture that provides an accurate effective description of the dynamics of the system in the low-energy sector of the Hilbert space. This is done through the Holstein-Primakoff (HP) transformation, shown as follows. The discretization of the quantum oscillator system is helpful for numerical simulations, as well as the assessment of the non-equilibrium thermodynamics of the process, as we explain in Sec. III.

Transformation from harmonic to angular momentum description
Computing dynamics of a quantum system with a varying potential under the action of an environment is in general difficult due to the infinite dimensions of the Hilbert space. The general non-Gaussian nature of the time-dependent potential prevents the use of techniques typical of Gaussian-state analysis, requiring the consideration of high-order moments of the position and momentum operators for a faithful account of the features of the system.
In order to bypass these bottleneck and provide an agile and physically intuitive picture of the dynamics at hand and its consequences, we aim to describe the system and the potential in Eqs. (4) and (5) through an effective picture based on the physics of angular momenta. With this aim in mind, we proceed as follows: 1. we discretize the system, and 2. we write the discretized system in terms of angular momentum operators.
Step 1: Discretization of the Hamiltonian can be achieved through the Holstein Primakoff (HP) transformation [33][34][35], which intorduces an effective bosonic system with a fixed dimension N . Such HP bosonic system through which one can define the HP position and momentum operators mω/2. Note that the commutation relation between such operator reads [b, b † ] = I−N |N N |, where I is the identity matrix. We are now able to discretize the system by recasting Eqs. (4) and (5) through the HP bosonic operators as x → x HP and y → y HP . The mapping well approximates the low energy sector as far as the value of N is sufficiently large. For N → ∞ the mapping is exact.
Step 2: Now, we consider a spin-J system with spin operators {J x , J y , J z , J 2 } (the same approach works for general angular momentum operators), we denote with j and j z the quantum numbers of J 2 and J z , respectively, whose simultaneous eigenstates |j, j z will be used as computational basis of our problem. The action of J z and J ± = J x ± iJ y is defined as usual through Moreover, the standard commutation relations are satisfied [J i , J k ] = i ikm J m , where {i, j, k} = {x, y, z} and ikm is the Levi-Civita symbol. To connect this algebra to that of the discrete system, we fix the value of j = (N − 1)/2, while j z = −j, −j + 1, . . . , j − 1, j. Then, the HP transformation imposes the relations between the HP bosonic operators and the spin operators as We notice that J + and J − are respectively proportionally mapped to b and b † . To express the HP bosonic operators with spin operators, we perform the Taylor expansion on the non-linear term in Eq. (8) up to the order κ, thus obtaining where M κ is a κ-th order polynomial in b † b. To make an explicit example, with an expansion up to the second order, we have For any order κ, M κ is a real and diagonal matrix, whose inverse matrix M −1 κ can be calculated. Then, we can approximate and correspondingly the position and momentum operators become This approximation becomes exact for κ → ∞. Finally, Step 2 is performed by the following mapping x hp → J x and p hp → J y . Thus, we can rewrite the dynamics in terms of the spin operators. In particular, the system Hamiltonian defined in Eq. (4) and Eq. (5) becomes with H sp add (t) = −E α(t) +ᾱ(t) where the label sp stands for spin. Similarly, one can easily recast the two dissipators with the spin ladder operators to and As we already pointed out, we are focusing on the lowenergy sector, thus the value of N can be rather low. To have an accurate description of the first 8 energy levels, we can take the spin system with N = 25 levels (namely j = 12) and the Taylor expansion to the order κ = 15. In Fig. 1c, we compare the first 8 energy levels of the continuous Hamiltonian defined in Eq. (4) with those of the Hamiltonian in terms of the discrete spin degrees of freedom as constructed in Eq. (13), which underlines the good approximation of our approach throughout the time of the process.

III. WEHRL ENTROPY IN TERMS OF Q-REPRESENTATION
To characterise the non-equilibrium dynamics of the system, we focus on the entropy production of the process under study. Here, we consider the Wehrl entropy, which is qualitatively similar to the von Neumann entropy, although better at capturing non-equilibrium features beyond the Gaussian regime. Generally, the Wehrl entropy provides an upper bound to the von Neumann one, and the two coincide for coherent states. The Wehrl entropy can be calculated through spin-coherent states, which are defined as [29,36,37] |Ω = e −iφJz e −iθJy |j, j .
Here |j, j is the angular momentum state with largest quantum number of J z , and Ω = {θ, φ} is the set of Euler angles identifying the direction along which the coherent state points. Thus, the Wehrl entropy for a system N = 2j + 1 degrees of freedom, as it is the one under investigation, reads [28,36,38] where dΩ = sin(θ) dθ dφ, while Q(Ω) is the Husimi-Q function associated with the state ρ, which is defined as Using the spin-coherent state picture, the key terms in the dynamical equations of motion formulated in the density-matrix representation can be mapped to the phase-spaces representation, such that With such correspondences at hand, we are allowed to access thermodynamic aspects of the system, namely irreversible entropy production and entropy flux rates, as stated in the following. Eq. (1) can now be reformulated in terms of the Husimi function as with the terms related to the non-unitary mechanisms being D i (Q) = Ω| D i [ρ] |Ω (i = th, lc), and U Q the transformed version of the term − i [H s , ρ] inducing the unitary evolution of the system. Correspondingly, one obtains three contributions to the rate of the Wehrl entropy production, which can be identified by merging the time derivative of Eq. (18) with Eq. (21), thus finding The dissipative Wehrl entropy rates connected to the dissipators can be further decomposed into the corresponding irreversible entropy production rate Π and entropy flux rate Φ as [29,36] dS The second thermodynamic law states an ever-increasing irreversible entropy production rate Π, while the whole entropy production rate of the system dS D / dt can be negative if one has a sufficiently large positive entropy flux rate Φ. The explicit forms of these two terms for the thermalisation dissipation read [36] Φ th = γ j(2j and For the localization process one has only the irreversible entropy production with no entropy flux. Thus, which was derived in Ref. [36] for the localization in J z basis while we derived that in J x basis in Appendix A. Then, the entropy associated to this process, strictly connected to the loss of coherence in J x basis, is always increasing.

IV. ENTROPY PRODUCTION FOR SIMULATED DYNAMICS
This section presents the main results of our study, which are the numerical simulations of the dynamics of the system and the corresponding analysis of its thermodynamics. Throughout the simulations, we set the parameters = ω = m = k B = 1, and the processing length time τ ω = 10.
The system is subject to the thermalisation and localisation dissipators, and next we show the effects of these dissipators in terms of their thermodynamic aspects, i.e. the Wehrl entropy, the irreversible entropy production rate and the entropy flux rate. To fully understand their respective actions, we look at the thermodynamics in different coupling regimes. When both the dissipators are considered, the system is driven to different non-equilibrium steady states depending on the relative coupling strengths.
Before digging into the problem at hand, we start the simulations with a study case of a spin system, which is simulated under three different conditions: when it is under the action of thermalization; that of localization; and finally under the combined action of the two dissipators. Then, we move back to the quantum oscillator case, where we instead focus on two cases: when the system is isolated, where we highlight the effects ot the dynamical change in the potential, and when the system is subject to the combined action of the two dissipators.

A. Spin-J systems
To understand the influences of the two dissipators, we consider the study case of a 4-level spin system, whose Hamiltonian H s is given by where J z is the spin operator along z direction. The initial state is prepared to be the thermal state ρ i = e −βHs Z(ω,β) with Z(ω, β) = i e −βei being the partition function, and {e i } are the eigenvalues of the system Hamiltonian. The system is initially prepared with inverse temperature β = 2/ ω. In the following, we consider the processes of thermalisation and localisation separately, then we study their combined action.

Thermalisation
The system is attached to a thermal environment at the inverse temperature β B = 1/ ω. This leads to a heating process where the system is asymptotically led to the thermal state ρ th = e −β B Hs Z(ω,βB) with respect to the environmental temperature. The thermal dissipator is given in Eq. (15) and we set γ/ω = 0.5. The corresponding evolution of the populations of the system is shown in Fig. 2a. The increase of temperature leads to a system with less distinctive population distribution. Such process also changes the Wehrl entropy S th Q of the system, which is calculated in Eq. (18) and its evolution -always increasing -is shown in Fig. 2d. The decomposition of the Wehrl entropy rate into the entropy production rate Π and the entropy flux rate Φ is shown in Fig. 2e. They can be respectively calculated from Eq. 24, and Eq. (25).
The trend of Π th shows that the production of irreversible entropy of the system is slowing down, while that of Φ indicates that the flux of entropy from the environment to the system is decreasing as the system is getting thermalised. One can see that the irreversible entropy and the entropy flux pace differently during the growth of the Wehrl entropy. Indeed, the irreversible entropy rate Π th is quickly suppressed while the entropy flux rate Φ th lasts for a longer time. For longer times, both the rates go to zero, indicating that the state is fully thermalised.

Localisation
Given any initial state, the localisation dissipator in Eq. (16) destroys the coherence of the system in positions, namely in the J x basis. Since the system Hamiltonian in Eq. (27) does not commute with J x , their combined action drive the system into a mixed state. We start with the same initial state considered in the previous case, and set the localisation parameter to Λ/ω = 0.5. The corresponding evolution of the populations of the system is shown in Fig. 2b, where the final state is found to be the completely mixed state ρ f = I/4, which denotes that any coherence in the system is destroyed. The evolution of the Wehrl entropy S lc Q for this process is shown in Fig. 2d, and it is increasing asymptotically towards the value that corresponds to the thermal state at infinite temperature. For this process, we have that the entropy flux rate Φ lc is null. Thus, the Wehrl entropy only rate consists in the irreversible entropy production rate Π lc . Its evolution in time is shown in Fig. 2e. Moreover, we notice that, although we took γ = Λ, the irreversible entropy production rates of the localization and thermalization process are not equal but following relation holds Π lc > Π th .

Combined dynamics
Here we consider the dynamics of the system that is subject to the combined action of the thermalisation and the localisation dissipators. We consider the previous settings, namely the couplings are γ/ω = Λ/ω = 0.5 and the inverse temperature for the environment is β B = 1/ ω. The corresponding evolution of the populations of the system is shown in Fig. 2c. Since two dissipators drive asymptotically the system to different states, this will eventually end in a non-equilibrium steady state. This is manifested by the evolution of its Wehrl entropy and production rates. One can see in Fig. 2d that the Wehrl entropy grows quickly at the beginning of the process, due to collective action of both dissipators, and then it stabilises to an average value between those of the two single processes. The entropy production rates of the combined dynamics are shown in Fig. 2f. In contrast to the previous two cases with only one dissipator acting, here the rates no longer approaching zero but a stable positive value, underling that the system has reached a non-equilibrium steady state. Moreover, since the locali-sation dissipator is driving the system toward an infinite temperature state, the system is overheated with respect to the temperature of the environment. Thus, the environment plays the role of cooling the system, which is indicated by the transition of the entropy flux rate Φ th from negative (when is still heating the system) to positive values. On the other hand, the irreversible entropy rates Π th and Π lc remain positive throughout of the process. Once a stable situation is attained, the following relation holds which indicates that the system has stabilised at a nonequilibrium steady state. Here, the system is overheated, and feeds to the environment the excessive entropy produced by the positive irreversible entropy productions of the two dissipators. We now deepen the study of the non-equilibrium steady state under the action of the combined dynamics by changing the relative strength of the coupling constants of the two dissipators. Without loss of generality, we set again γ/ω = 0.5 and let the value of Λ change. The first quantity of interest is the coherence, which is quantified with a l 1 measure in the basis of J z as [39] C l1 (ρ) = j =k |ρ jk |.
We observe that at the end of the process there is still present some coherence in the final non-equilibrium steady state as it is shown in Fig. 3a. Here, the two dissipators drive the system in different bases: the system is prepared in basis of J z and the D lc [ρ] localizes the system in basis of J x ; due to their non-commutativity, the system is driven out of equilibrium during the process, thus maintaining some coherence. Such a coherence would disappear if only the localisation dissipator D lc [ρ] acts on the system. However, when both dissipators are present, the thermalisation dissipator D th [ρ] prohibits the system from being fully localised in basis J x , yet fails to fully thermalise the system, resulting in a final non-equilibrium steady state containing non-zero coherence. We see the coherence peaks at Λ/γ ≈ 2, and the coherence disappears when Λ = 0 or ∞, where the dynamics is dominated by one of the two dissipators. We proceed by quantifying the influences of the thermalisation and localisation dissipators in driving the system towards the fully thermalised state ρ th and the fully localised state ρ lc , respectively. To quantify this driving, we employ fidelity, which is defined as In particular, Fig. 3b shows the fidelity F (ρ f , ρ th ) between the final and the thermal state, with referencing lines that represent ρ th (the red dot line) and ρ lc (the gray dot-dashed line). On the other hand, we also demonstrate the Wehrl entropy of the final state w.r.t different coupling strengths, compared with Wehrl entropies of ρ th and ρ lc . From both plots, one can observe a non-linear transition as Λ increases. Finally, in Fig. 3c we compare the entropy production rates for different values of the ratio of the couplings. We analyse them at the end of the simulations, when the system reaches the non-equilibrium final state and the relation in Eq. (28) holds true. As the ratio Λ/γ increases, one can observe an increase of the irreversible entropy production rate Π th and the flux rate Φ th from the thermalisation dissipator. This is due to the localization dissipator, which pushes the system toward an "infinite temperature" state. On the other hand, one has a peak of the irreversible entropy production rate Π lc from the localisation dissipator at Λ/γ ≈ 2, for whose value the final state displays the maximum coherence, as it is shown in Fig. 3a.

B. Quantum oscillator with non-Gaussian potential
We now focus on the quantum oscillator with the timevarying Hamiltonian redefined in Eqs. (13) and (14), with dimension N = 25 and Taylor expansion up to the order κ = 15. The initial state is prepared as the thermal state with respect to H sp s (0) at β = 2/ ω (we remand below for the analysis with different initial states). Here we analyze two different situations: the case of the isolated system, which is not as trivial as in the case of the spin system, and the open dynamics under the influence of both the dissipators.

Isolated dynamics
In the case of the isolated dynamics, we are interested in verifying the effective splitting of the wavefunction distribution when going from the Gaussian-harmonic potential to the double-well. Thus, we first verify, in the J x basis, the proper separation of the initial population in the two peaks as the dynamics takes place. The plot in Fig. 4a shows the effects of the time-varying potential in the population distribution. Since at t = 0 the system is prepared at a low temperature, the population initially resides mostly in the lowest energy level, which is localized in a region being small with respect to the distance of the double well [cf. Fig. 1a]. Then, the population divides in two peaks corresponding to the two minima of the double-well potential at t = τ . Since the double-well potential has a shallower depth and larger eigenenergies [cf. Fig. 1c], the energy of the system increases during the transition, as described by Fig. 4b. Also the coherence of the system, which is quantified by Eq. (29) in basis of J x , grows in time. Indeed, the double-well potential, by splitting the population in two localized regions, generates coherence between these two that results in its increasing behaviour reproduced in Fig. 4b. Finally, we show the Wehrl entropy of the system and its rate respectively in Fig. 4c and Fig. 4e. They both capture the non-equilibrium features caused by the change in the potential. Initially, the small deformation of the potential from its quadratic behaviour pushes the state sightly out of the equilibrium. Consequently, the Wehrl entropy rate presents small oscillations, which are captured in the first half of the evolution presented in Fig. 4e. In the second half of such evolution, the potential starts to present the central peak characterising the double-well. Subse- shows the mean energy and the coherence, while the Wehrl entropy is displayed in Panel (c). Panel (e) shows the unitary contribution to Wehrl entropy production rate, while Panel (f ) exhibits the corresponding irreversible entropy production rates Π th , Π lc and the entropy flux rate Φ th for the combined dynamics at the beginning of the dynamics. Their entire evolution down to time t = τ is illustrated in its inset.
quently, the non-equilibrium driving becomes stronger and leads to a stronger increases of the Wehrl entropy of the system and its rate.

Combined dynamics
We now consider the combined dynamics where the system is under the action of both the thermalisation and the localisation dissipators. Here, we set the coupling rates to γ/ω = Λ/ω = 0.5 and the inverse temperature for the environment β B = 1/ ω, which are the same values considered in the study case considered in Sec. IV A. The main effect of dissipators is to wash away all the quantum features of the system, and thus it tries to nullify the action of the transition to the double-well. In Fig. 4d, we show the evolution of the population under this combined dynamics. As one can see, the population is modified only in the first moments of the evolution, where the distribution slightly broadens without splitting in the expected two peaks. The first part of the evolution corresponds to the initial increase of the energy driving the system in a higher temperature state by the two dissipators (cf. Fig. 4b). After this initial change, the distribution in position appears as "frozen" for the rest of the dynamics. It is shown in Fig. 4b, the energy of system increases with an initial surge driven by the two dissipators until the corresponding inverse temperature β B is reached. Then, the energy continues to increase but at a lower rate only due to the change of the potential. Indeed, the combined action of the two dissipators does not contribute to the energy change: although the localisation heats the system, the thermalisation cools it back to β B . Alongside with this, the coherence becomes smaller, which reflects the illegible effect of the potential on the population distribution, until it reaches a constant value as a result of the action of the two dissipators driving the system in two different basis (the same was noticed in the spin study case, cf. Fig. 3a). The Wehrl entropy of the system and its rates are shown respectively in Fig. 4c and Fig. 4f. Similarly as the spin situation, the Wehrl entropy increases due to the thermalisation and the localisation. The latter overheats the system with respect to the thermal dissipator, so that the entropy flux rate Φ th changes from being negative to positive. Eventually, the system reaches a non-equilibrium steady state for which the decomposed rates then follow Eq. (28) since the unitary rate dS U / dt is negligible when compared to the other rates.
a. Different coupling regimes The action of the two dissipators has prevented the generation of a population in position characterized by two peaks relative to the minima of the double-well potential. We thus proceed  to analyse which can be the suitable coupling strengths for which such a two-peak distribution can still be observed even though the presence of the two dissipators.
Here, we consider three scenarios: the one where one has only thermalisation, that with only localisation, and that with both the dissipators being present. Today, optical levitation experiments performed in low vacuum allow to reach a condition where photon recoil is significantly stronger then damping. We thus explore this regime fixing Λ/γ = 10 for the rest of these simulations. See Sec. V for detailed discussions. We consider two figures of merit to make a full comparison of these scenarios. The first one is the fidelity F (ρ f , ρ is f ) between the non-equilibrium steady state for the considered dynamics and the final state for the isolated dynamics, and it is shown in Fig. 5a. One can infer that the previous simulation (where Λ/ω = γ/ω = 0.5) reproduced a low-fidelity final state. To improve such a situation, one is required to reduce the values of Λ and γ. In particular, for being able to reach a fidelity of F > 80%, one needs to set the coupling strength γ or Λ to approximately 10 −3 ω in the case of a scenario with a single dissipator. The requirement becomes instead more stringent, with γ/ω ≈ 10 −4 and Λ/ω ≈ 10 −3 , when both the dissipators are involved in the dynamics. These values are close to the current state-of-art experimental setting (See Sec. V). Notably, the coherence C l1 , which is defined in Eq. (29), provides approximately the same behaviour as the fidelity, up to a proportionality constant. We display its scaling in the right vertical axis of Fig. 5a. Finally, we consider the Wehrl entropy of the final state as the second figure of merit, which is shown in Fig. 5c. When considering the two single-dissipator cases, one notices that the thermalisation has a milder action on the entropy compared to that of the localisation. Indeed, the thermalisation dissipator heats the system up to the inverse temperature β B , but not beyond that. On the other side, the localisation dissipator, which can be treated as a bath at infinite temperature, keeps heating the system. Consequently, also its action in terms of entropy is stronger. In the case where both the dissipators are considered with the fixed relation of Λ = 10γ, one can identify two different regimes: up to γ/ω ∼ 10 −2 both dissipators slowly heat the system; beyond that value, the thermalisation works as a cooling mechanism, in contrast to the localisation that keeps heating the system. This results in a increasing behaviour of the Wehrl entropy up to γ/ω ∼ 10 −2 , after which it is slowly descending.
b. Low coupling regime. Given the results of the analysis above, we consider a low coupling regime that can potentially showcase stronger quantum features at the end of the protocol, and maintain them also after it. In particular, we set Λ = 10γ and consider two values for γ: case 1. γ/ω = 10 −3 , and case 2. γ/ω = 10 −2 . For both cases, we extend the dynamics after the transformation protocol to t = 4τ . This gives an insight on the environmental effects on the system also after the change of the potential.
In the low-coupling regimes, the first thing to notice is that we observe again the double-well effect in the population. See for example the population of case 1. in Fig. 6a where one can see the double-peak distribution also for t = 4τ . In both cases, the thermodynamics of the system gets comparable influences both from the unitary dynamics due to the potential and the action of the environment. Indeed, the magnitude of the unitary contribution dS U dt and the dissipative one dS D dt to the Wehrl entropy rate are compared in Fig. 6b and Fig. 6c respectively for the case 1. and 2. From these figures, one observes that, for larger the coupling strength, the difference between the unitary and the dissipative contributions gets larger. This is the reason why, in the previous simulation with γ/ω = 0.5, it was possible to neglect the unitary part and only investigate the dissipative contribution. Moreover, after the end of the transition to the double-well potential, the system is still in contact with the environment and it tends to the steady state with a rate depending on the values of the coupling strengths. Such a process suggests that the total rates of Wehrl en- tropy production is going to zero asymptotically. However, with smaller coupling strengths, one can no longer ignore the unitary contribution dS U dt . This causes an interesting situation, which is shown in Fig. 6c. The system tends to a non-equilibrium steady state that is neither in equilibrium with the system Hamiltonian nor the environment. Indeed, one gets dS U dt > 0 and dS D dt < 0, which means that the system keeps producing entropy due to non-equilibrium with the the Hamiltonian, while in meantime it dissipates the entropy in the environment.
c. Different initial states Due to the Gaussian term in the Hamiltonian in Eq. (5), the initial state of the single-well potential is naturally slightly squeezed in position, as it is shown in Fig. 6d. Here, we want to quantify which are the effects of such a squeezing on the performance of the protocol under study. To this end, we consider other initial states by varying the amount of initial squeezing and study if this can provide a better outcome. The squeezing operator is given by [40] and the squeezed initial state is prepare by ρ sq in = S(ζ)ρ in S(ζ) † with ζ ∈ [−1, 1]. In particular, we consider a similar simulation as above with the values of the parameters fixed at γ = 10 −3 , Λ = 10γ. Firstly, we show in Fig. 6e the energy and coherence of the system at t = τ at different values of initial squeezing ζ. Regardless of the quadrature being squeezed, the initial squeezing adds energy to the system. On the other hand, a suitable squeezing can increase the coherence. The maximum in coherence is given by ζ = −0.2, which corresponds to a small squeezing in momentum. Taking such a squeezed state as the initial state, we simulate the system and compute the dynamics of unitary and dissipative Wehrl entropy rates, which are shown in Fig. 6f. Comparing them qualitatively to the non-squeezed case (cf. Fig. 6b), the squeezing induces huge oscillations to the unitary Wehrl entropy rate, while it does not affect significantly the dissipative rate. After the deformation of the potential, both rates approach to zero as in the non-squeezed case.

V. EXPERIMENTAL PERSPECTIVE
The experimental implementation of time-controlled double-wells in the quantum regime is challenging. Yet, recent experiments with optically levitated nanoparticles provide the necessary ingredients. Quantum limited position readout and the preparation of nearly pure quantum states of motion [24][25][26] as well as the dynamical shaping of double-well potentials [27] have been demonstrated. Thus, the protocol discussed here is in principle implementable, once these two methods are combined. This requires the implementation of existing state detection methods [25,26] to the detection of the transversal motion [41], where the double-well potential landscape is realized. Yet, the downside in optical levitation is the unavoidable decoherence due to recoil of optical tweezer photons, resulting predominantly in a localization term. Additionally, thermalization with the environment may be induced by collisions with the surrounding gas molecules. In the following, we base the discussion on the currently most favorable parameters in this respect, achieved by the recent cryogenic implementation of an optical harmonic trap by Tebbenjohanns et al. [26].
a. System and potential. We consider a silica nanoparticle of radius R = 50 nm trapped in an dynamically controlled optical double well potential, created from a combination of a TEM 00 and TEM 01 optical mode [27]. The trapping laser has a wavelength of λ = 1550 nm. We assume a frequency of ω = 2π × 54.9 kHz in Eq. (4), which corresponds to an effective optical trapping frequency of 2π × 77.6 kHz when we add also the contribution from H add with α = 1. The timescale of the protocol τ is limited by the switching times of the light field, which can be several orders of magnitude faster than the mechanical frequency when done with electrooptic modulators. It can therefore be considered near instantaneous. Using a numerical aperture of the optical tweezer N = 0.75 the beam waist is W ≈ 660 nm. Assuming equal restoring forces contributing to H s and H add we find E = 8 × 10 5 ω in Eq. (5).
b. Localisation dissipator. The major effect of localization is due to the scattering of tweezer photons from the levitated nanoparticle. We assume that this effect in the optical trap exceeds localization due to blackbody radiation, which is expected to be fulfilled in current experiments. The corresponding decoherence parameter at the maximum intensity of a Gaussian beam is given by [42]: where c = −1 3( +2) and is the relative dielectric constant of the nanoparticle. The above formula results in Λ ≈ 2.4 × 10 −4 ω from photon recoil.
The surrounding gas pressure of 3 × 10 −9 mbar results in a damping rate of γ exp = 5.9 × 10 −5 Hz at an environmental temperature of 60 K (thermal occupation n th = 10 7 ) resulting in a thermalisation rate γ = n th γ exp ≈ 2 × 10 −3 ω. As a further reduction in pressure by 3 orders of magnitude seems feasible [26], the thermalization rate can still be significantly suppressed below the localization rate as discussed in Appendix B.
Accordingly, the effect of the environment is already reduced to the relevant range considered in Fig. 5. Note, however, that the experimental parameter for the strength of the double-well Hamiltonian H add of E = 8 × 10 5 ω significantly exceed the value chosen for the simulations of E = 10 ω. On an experimental level, reducing the power sufficiently to match that constraint corresponds to performing the experiment at the level of ω ≈ 1 Hz. While the normalized localization rate would not alter, a corresponding reduction of the pressure by 5 orders of magnitude would be required and also the feasibility of such optical traps is unclear. From a different perspective, this is however not necessary, as it may be expected that similar effects as discussed here will emerge also for large E. Theoretical analysis of this behaviour however requires treatment of a significantly larger number of energy levels, which is hard to treat numerically.

VI. CONCLUSIONS
In this paper, we constructed an approach to study the thermodynamics a nanoparticle under the combined action of a time-varying potential, which reproduce the transition from a harmonic to double-well landscape, and external dissipative influences. To simulate numerically its dynamics, we needed to discretise the system without loosing its dynamical features, a considerably complex task due to the presence of the non-trivial timevarying potential. Here, we provide a general approach to implement such a discretisation at a suitable level of approximation, and that can be also used to characterise the non-equilibrium thermodynamics of the system. The latter is here provided in terms of the Wehrl entropy, which well captures the non-Gaussian features of the system, and its production rates. The analysis showed that, starting from a thermal initial state relative to the initial potential, a non-equilibrium steady-state is eventually achieved. Moreover, after suitably setting the coupling strengths of the thermalisation and localisation processes, one can arrive at a superposition state in position, namely a non-equilibrium high-coherence steadystate exhibiting two delocalised peaks in its position distribution. We demonstrated that this state survives the end of the variation of the potential and maintains its quantum features. By suitably squeezing the initial thermal state, one can further improve the final coherence value.
Finally, we investigated the feasibility of implementing the proposed scheme in optical-levitation experiments. While state-of-the-art experiments operate at the relevant levels of decoherence, our study addresses a change in the trapping potential -from quadratic to double-well -involving low-energy states that is currently experimentally unexplored. It would be interesting to further our analysis to extend it to the context of electrostatic or magnetic levitation, and nanomechanical resonators (cf., for example, Ref. [43]), which might offer more favourable working regimes.