Single-atom heat engine as a sensitive thermal probe

As depicted in figure 1(c), an atomic ion is trapped in a Paul trap consisting of four blade-shaped electrodes and two endcap-electrodes [27]. In contrast to a conventional Paul trap, the blade electrodes are tilted with respect to the axial symmetry axis by an angle θ. The four blade-electrodes are driven by radio-frequency signals and the two endcap-electrodes are biased with positive DC voltages. The ponderomotive potential formed by this trap can be approximated as harmonic in the radial (x, y) and axial (z) directions with a variation of the radial trapping frequencies in the axial direction:

$$\begin{equation}{\omega }_{x,y}\left(z\right)=\frac{{\omega }_{x0,y0}}{{\left(1+\mathrm{tan}\enspace \theta \cdot z/{r}_{0}\right)}^{2}},\end{equation} \tag{ 1 }$$

with r₀ the radial distance of the ion to the blade electrodes and ω_x0,y0 the radial trapping frequencies at z = 0. In an experiment, a small anisotropy of the radial potential is normally present, which lifts the degeneracy of ω_x0,y0. The confinement in the axial direction is much weaker compared to that of the radial directions, signified by a smaller axial trapping frequency ω_z ≈ ω_x0,y0/10.

The Hamiltonian describing the motional state of an ion in this potential can be written as

$$\begin{equation}\hat{H}=\frac{{\hat{p}}_{x}^{2}}{2m}+\frac{1}{2}m{\omega }_{x}^{2}{\hat{x}}^{2}+\frac{{\hat{p}}_{y}^{2}}{2m}+\frac{1}{2}m{\omega }_{y}^{2}{\hat{y}}^{2}+\frac{{\hat{p}}_{z}^{2}}{2m}+\frac{1}{2}m{\omega }_{z}^{2}{\hat{z}}^{2},\end{equation} \tag{ 2 }$$

with m the mass of the ion, $\hat{x}$, , and ${\hat{p}}_{x}$, ${\hat{p}}_{y}$, ${\hat{p}}_{z}$ the position and momentum operators in the corresponding directions. Viewed as a heat engine, the radial states represent the working medium which thermalizes with external hot and cold baths periodically. The engine drives and amplifies a coherent oscillation of the ion in the axial direction. The axial oscillator thus serves as a flywheel that stores the energy output. In open-cycle operations, the flywheel is not actively damped. As we will show in later sections, the energy stored in the flywheel grows quadratically with the number of engine cycles. In our analysis, we will treat the large-amplitude coherent oscillation in the axial direction as a classical oscillator, while keeping a quantum mechanical description for the radial states. The Hamiltonian of the system can then be rewritten as

$$\begin{align}\hfill \hat{H}& =\hslash {\omega }_{x0}\left({\hat{a}}_{x}^{{\dagger}}{\hat{a}}_{x}+\frac{1}{2}\right)+\hslash {\omega }_{y0}\left({\hat{a}}_{y}^{{\dagger}}{\hat{a}}_{y}+\frac{1}{2}\right)+\frac{1}{2}m\left({\omega }_{x0}^{2}{\hat{x}}^{2}+{\omega }_{y0}^{2}{\hat{y}}^{2}\right)\left(\frac{1}{{\left(1+\gamma z\right)}^{4}}-1\right)+\frac{{p}_{z}^{2}}{2m}+\frac{1}{2}m{\omega }_{z}^{2}{z}^{2},\hfill \end{align} \tag{ 3 }$$

where ${\hat{a}}_{x}^{{\dagger}}$ (${\hat{a}}_{y}^{{\dagger}}$) and â_x (â_y ) represent the bosonic creation and annihilation operators of the phonons in x (y) direction, respectively. The first two terms in equation (3) describe the energy stored in the working medium. The third term represents the coupling of the working medium to the flywheel, where we define γ = tan θ/r₀. The last two terms denote the kinetic and potential energy in the flywheel, where p_z and z are the classical momentum and position of the ion in the axial direction.

As indicated by equation (3), the two radial directions contribute to the Hamiltonian independently. For the sake of brevity, we consider only the x direction in our further analysis, while assuming that the oscillator in the y direction is maintained at a low and constant temperature. We denote the creation and annihilation operators in the x direction as â^† and â. The equations of motion for the quantum harmonic oscillator in the x direction read

$$\begin{align}\hfill \begin{aligned}\hfill \frac{\mathrm{d}}{\mathrm{d}t}\hat{X}\left(t\right)& =2\left(\hslash {\omega }_{x0}+2g\left(t\right)\right)\hat{Y}\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}}{\mathrm{d}t}\hat{Y}\left(t\right)& =-2\left(\hslash {\omega }_{x0}+2g\left(t\right)\right)\hat{X}\left(t\right)-8g\left(t\right)\hat{N}\left(t\right)-4g\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}}{\mathrm{d}t}\hat{N}\left(t\right)& =-2g\left(t\right)\hat{Y}\left(t\right),\hfill \end{aligned}\end{align} \tag{ 4 }$$

where we define

$$\begin{align}\hfill \begin{aligned}\hfill \hat{X}\left(t\right)& =\left({\hat{a}}^{{\dagger}2}\left(t\right)+{\hat{a}}^{2}\left(t\right)\right)\hfill \\ \hfill \hat{Y}\left(t\right)& =i\left({\hat{a}}^{{\dagger}2}\left(t\right)-{\hat{a}}^{2}\left(t\right)\right)\hfill \\ \hfill \hat{N}\left(t\right)& ={\hat{a}}^{{\dagger}}\left(t\right)\hat{a}\left(t\right)\hfill \\ \hfill g\left(t\right)& =\frac{\hslash {\omega }_{x0}}{4}\left(\frac{1}{{\left[1+\gamma z\left(t\right)\right]}^{4}}-1\right).\hfill \end{aligned}\end{align} \tag{ 5 }$$

The classical equations of motion for the flywheel can be written as

$$\begin{align}\hfill \dot{z}\left(t\right)& ={p}_{z}\left(t\right)/m\hfill \\ \hfill {\dot {p}}_{z}\left(t\right)& =\mathcal{F}\left(t\right),\hfill \end{align} \tag{ 6 }$$

with $\mathcal{F}$ the force acting on the axial oscillator. In the mean-field approximation, $\mathcal{F}$ can be expressed as

$$\begin{equation}\mathcal{F}\left(t\right)=-\frac{\partial \langle \hat{H}\rangle }{\partial z}=-m{\omega }_{z}^{2}z\left(t\right)+F\left(t\right),\end{equation} \tag{ 7 }$$

which is the sum of the restoring force of the harmonic potential and the force F resulting from the radial–axial coupling

$$\begin{equation}F\left(t\right)=\frac{\gamma \hslash {\omega }_{x0}R\left(t\right)}{{\left(1+\gamma z\left(t\right)\right)}^{5}}-\frac{\hslash {\omega }_{x0}}{4}\frac{\partial R}{\partial z}\left(\frac{1}{{\left(1+\gamma z\left(t\right)\right)}^{4}}-1\right).\end{equation} \tag{ 8 }$$

Here, we define $R\left(t\right)=\langle {\left({\hat{a}}^{{\dagger}}\left(t\right)+\hat{a}\left(t\right)\right)}^{2}\rangle$. We note that, one could alternatively assume the same temperature of the two radial modes, which would result in doubling of the force F. In typical experimental scenarios [27], the small axial displacement z(t) ≪ r₀ holds, leading to γz(t) ≪ 1. This allows the approximation

$$\begin{equation}\mathcal{F}\left(t\right)=-m{\omega }_{z}^{2}z\left(t\right)+\gamma \hslash {\omega }_{x0}R\left(t\right),\end{equation} \tag{ 9 }$$

with

$$\begin{align}\hfill R\left(t\right)& =\langle \hat{X}\left(t\right)+2\hat{N\left(t\right)}+1\rangle \hfill \\ \hfill & =1+2{N}_{0}+{X}_{0}\enspace \mathrm{cos}\left(2{\omega }_{x0}t\right)+{Y}_{0}\enspace \mathrm{sin}\left(2{\omega }_{x0}t\right),\hfill \end{align} \tag{ 10 }$$

and X₀, Y₀, N₀ the expectation values of $\hat{X}\left(t\right)$, Ŷ(t) and $\hat{N}\left(t\right)$ at t = 0.

The model can be cast to describe the dynamics of a four-stroke Otto engine [31–34]. One engine cycle, as shown in figures 2(a) and (b), can be described as follows: hot isochore [A]: the radial trapping frequency is kept constant at ω_h and the working medium is in contact with a hot bath at temperature T_h. After a time of τ_h, the working medium thermalizes at this temperature. The time scale of the interaction is much shorter than the half period of the axial oscillation τ_z = π/ω_z . During this stroke, the axial displacement of the ion and the change of the radial trapping frequency are negligible. Isentropic expansion [B]: the ion is isolated from the baths and evolves in the trapping potential. After a time of τ_z , the radial trapping frequency changes from ω_h to ω_c due to the displacement in the axial direction. This trapping frequency change may be either adiabatic or engineered using optimal control techniques [35] that also allow for the minimization of different sources of error in the variation of ω_x [36, 37]. The dynamics describing this process is unitary and the entropy of the working medium remains constant. Cold isochore [C]: the radial trapping frequency remains constant at ω_c. The ion is coupled to a cold bath and its radial state thermalizes to the temperature T_c (T_c < T_h) after time τ_c (τ_c ≈ τ_h ≪ τ_z ). Isentropic compression [D]: in the last stroke, the ion is again isolated from the baths and evolves isentropically in the trapping potential for another time of τ_z . In reference [27], the isochoric strokes were realized by coupling the ion to a laser or external electric fields, which emulate the cold and hot baths, respectively. The cold bath is realized by exposing the ion to a laser beam red detuned from the ² S_1/2–² P_1/2 transition of the calcium ion. The hot bath is emulated by applying noises to the external DC electrodes, while simultaneously exposing the ion to the cooling laser. The noise is generated and timed by arbitrary waveform generators. The two baths provide an accessible temperature range of 1 mK to 4 K for the working medium.

Zoom In Zoom Out Reset image size

In the following, we analyze the dynamics of the engine, presenting both the approximated analytical and the full numerical solutions of the Hamiltonian. We establish a numerical routine using a combination of propagation of the radial wavefunctions following the Liouville–von-Neumann equation [38, 39] and classical trajectory simulation in the axial direction with partitioned Runge–Kutta method [40]. This numerical platform allows us to simulate the dynamics of the engine driven by thermal and non-thermal baths. Details about this numerical routine are given in appendix A.

We consider typical experimental parameters of a ⁴⁰Ca⁺ ion with m = 40 amu confined in a tapered Paul trap with θ = π/6, r₀ = 1 mm and trapping frequencies of ω_x0 = 2π × 1 MHz and ω_z = 2π × 0.1 MHz. Before starting the engine, all the motional states are initialized using Doppler cooling to a thermal state at 1 mK. This results in more than 200 average phonon occupation in the axial oscillator, which justifies a classical treatment in this direction. After the initialization and at time t = 0, the ion is at an axial position z₀ with axial velocity v_z0 = p_z0/m. At this time, a hot thermal bath of temperature β_h = 1/k_B T_h is switched on. The working medium thermalizes to this temperature and has ${X}_{0}={Y}_{0}=0,\enspace {N}_{0}={\left({\mathrm{e}}^{{\beta }_{\mathrm{h}}\hslash {\omega }_{x0}}-1\right)}^{-1}$. The force due to the radial–axial coupling in the small axial displacement limit becomes

$$\begin{equation}{F}_{\mathrm{h}}=\gamma \hslash {\omega }_{x0}{R}_{\mathrm{h}},\end{equation} \tag{ 11 }$$

with

$$\begin{equation}{R}_{\mathrm{h}}=\mathrm{coth}\left(\frac{{\beta }_{\mathrm{h}}\hslash {\omega }_{x0}}{2}\right).\end{equation} \tag{ 12 }$$

This force leads to a displacement of the axial potential toward the open end of the taper. After the interaction with the bath, the radial state expands isentropically in the potential. The classical trajectory of the ion in the z direction can be solved by integrating the equations of motion (6). After τ_z , the ion reaches the position z₁ at t = τ_z with velocity v_z1. At this point, a cold bath with temperature β_c = 1/k_B T_c is switched on and cools the radial state to this temperature. This process is again isochoric. After the cooling process, the force is reduced to

$$\begin{equation}{F}_{\mathrm{c}}=\gamma \hslash {\omega }_{x0}{R}_{\mathrm{c}},\end{equation} \tag{ 13 }$$

where R_c is given by equation (12) and exchanging β_h → β_c. The ion experiences a displaced axial potential toward the narrower end of the trap. For the next half axial oscillation period, the ion evolves again isentropically in the displaced potential. At the end of a completed engine cycle t = 2τ_z , the flywheel does not restore the original (z₀, v_z0) point but ends at (z₂, v_z2) due to the work done by the forces F_h and F_c. The energy–frequency diagram for the radial states undergoing one engine cycle is illustrated in figure 2(b).

In figures 2(c) and (d) we present the results of numerical simulations with T_h = 1.2 mK and T_c = 1.0 mK. The blue line in figure 2(c) shows a close-up of the axial trajectory of the ion as a function of time. The oscillation amplitude grows linearly with the number of engine cycles. The inset displays the full axial trajectory over four engine cycles, where the blue line shows the numerical simulation of the exact Hamiltonian in equation (3) and the dashed red line represents the approximated analytical solution using equations (11) and (13). The analytical and numerical results show excellent agreement. Figure 2(d) displays the evolution of the energy in the flywheel (red curve) and the working medium (blue curve). As a result of the linear growth of the axial oscillation amplitude, the energy in the flywheel increases quadratically with the number of engine cycles.

The measurement protocol is illustrated in figure 3. To begin with, all the motional states are initialized by Doppler cooling to a thermal distribution of temperature T₀ = 1.0 mK. The initial axial position z₀ and velocity v_z0 of the ion follow the Boltzmann distribution $f\left({z}_{0},{v}_{z0}\right)\propto \mathrm{exp}\left[-m\left({v}_{z0}^{2}+{\omega }_{z}^{2}{z}_{0}^{2}\right)/2{k}_{\mathrm{B}}{T}_{0}\right]$. After the initialization, the engine is set to the four-stroke operation driven by two baths of temperature T₁ and T₂, respectively. The position of the ion upon its nth interaction with the baths can be obtained by integrating the equations of motion in equation (6) stroboscopically. In a typical experimental scenario, where the axial oscillation amplitude is much smaller compared to the radial dimension of the trap (z ≪ r₀), we arrive at

$$\begin{align}\hfill \begin{aligned}\hfill {z}_{n}& ={z}_{0}+\frac{n\hslash {\omega }_{x0}\gamma \left({R}_{2}-{R}_{1}\right)}{m{\omega }_{z}^{2}}\quad \text{for}\enspace n=0,2,4\dots \hfill \\ \hfill {z}_{n}& =-{z}_{0}+\frac{2\hslash {\omega }_{x0}\gamma {R}_{1}-\left(n-1\right)\hslash {\omega }_{x0}\gamma \left({R}_{2}-{R}_{1}\right)}{m{\omega }_{z}^{2}}\quad \text{for}\enspace n=1,3,5\dots ,\hfill \end{aligned}\end{align} \tag{ 14 }$$

where R₁ and R₂ are given by equation (12) with the corresponding temperatures T_1,2. The change in the axial oscillation amplitude of two consecutive peaks (see figure 2(c))

$$\begin{equation}{\Delta}z={z}_{n+2}-{z}_{n}=\frac{2\hslash {\omega }_{x0}\gamma }{m{\omega }_{z}^{2}}\left({R}_{2}-{R}_{1}\right),\end{equation} \tag{ 15 }$$

which in the limit of ℏω_x0 ≪ k_B T simplifies to

$$\begin{equation}{\Delta}z=\frac{4{k}_{\mathrm{B}}\gamma }{m{\omega }_{z}^{2}}{\Delta}T,\end{equation} \tag{ 16 }$$

where ΔT = T₂ − T₁. We note that the assumptions on the axial oscillation amplitude and the bath temperatures are only made to arrive at the analytical solutions of the equations of motion. The function principle of the thermal probe remains valid out of these limits and the equations of motion can be solved numerically.

Zoom In Zoom Out Reset image size

As shown in equation (16), Δz carries the same sign with T₂ − T₁. After N(N ≫ z₀/Δz) engine cycles, the energy fed by the engine dominates over the initial energy in the flywheel. The axial oscillation will thus be brought in phase with the bath interactions. The ion travels to the open side of the tapered potential (z > 0) when put in contact with the bath of lower temperature, and vice versa. The temperature difference between the two baths can be determined by measuring the axial oscillation amplitude. We propose two sets of triggered measurements as follows. The first set of measurements is performed by illuminating the ion with a short laser pulse of duration τ_z /10 and at a time τ_z after the 2Nth bath interaction (see upper panel of figure 3(a)). Fluorescent light from the ion is collected by a microscope objective and imaged on a camera (see figure 3(b)). The duration of the laser pulse is chosen to be much shorter than the axial oscillation period to minimally disturb the position of the ion. To obtain a significant signal-to-noise ratio, the measurement has to be reproduced M times with a restart of the engine each time. The mean position of the ion is determined by performing single-particle localization analysis [41] on the integrated fluorescent image. We note that each restart of the engine samples a random z₀ and v_z0 from the Boltzmann distribution. As shown in figure 4(a), the influence of the distribution of z₀ averages out due to a large number of sampling. The measured mean axial position of the ion is thus

$$\begin{equation}{\overline{z}}_{2N}=\frac{2N\hslash {\omega }_{x0}\gamma \left({R}_{2}-{R}_{1}\right)}{m{\omega }_{z}^{2}}.\end{equation} \tag{ 17 }$$

The second set of measurements is performed following the same procedure but at a time 2τ_z after the 2Nth bath interaction, while the (2N + 1)th bath interaction is skipped (see lower panel of figure 3(a)). The measured mean axial position is then

$$\begin{equation}{\overline{z}}_{2N}^{\prime }=-\frac{2N\hslash {\omega }_{x0}\gamma \left({R}_{2}-{R}_{1}\right)}{m{\omega }_{z}^{2}},\end{equation} \tag{ 18 }$$

leading to the difference between the two positions

$$\begin{equation}{\overline{z}}_{2N}-{\overline{z}}_{2N}^{\prime }=2N{\Delta}z,\end{equation} \tag{ 19 }$$

which is a linear function of the temperature differences ΔT.

Zoom In Zoom Out Reset image size

Figure 4(a) displays a simulated distribution of the axial positions of the ion upon the two sets of measurements with ΔT = 0.1 mK, T₀ = 1.0 mK, N = 10⁵ and M = 2 × 10⁵. Figure 4(b) shows the linear dependence of ${\overline{z}}_{2N}-{\overline{z}}_{2N}^{\prime }$ with the temperature difference ΔT for N = 10⁵. The black line displays the results of the approximated analytic solution. The green circles and red crosses display the same value deduced from numerical simulations performed at T_c = 1.0 mK and 0.2 mK, respectively. The good agreement of the numerical simulations confirms the validity of equation (16) at different base temperatures.

The precision in determining ΔT is limited by experimental uncertainties in measuring the amplitude 2NΔz. In single-particle localization analysis, the uncertainty in the center position is given by the signal-to-noise ratio of the camera image. This is determined by the illumination time of the laser, the efficiency of photon collection, the quantum efficiency, and the background noise of the camera. In reference [27], an objective with a numerical aperture of 0.26 was used to collect the fluorescent photons at 397 nm and to form an image on an intensified charge-coupled-device sensor. At each oscillation phase, 2 × 10⁵ repeated measurements were necessary to obtain a localization precision of ±250 nm [42]. Assuming the same measurement conditions and using equations (16) and (19), we arrive at an uncertainty of ±2 μK in determining ΔT.

We note that in the limit of ℏω_x ≪ k_B T_1,2, the outcome of the measurement depends only on the temperature difference between the two baths and is independent of their absolute temperature. The protocol is particularly suitable for detecting small temperature differences at high base temperatures. On the other hand, when using a bath of well-characterized temperature as a reference, the absolute temperature of an unknown bath can be determined.

The amplitude of the measured signal can be amplified by exploiting quantum resources such as squeezing of the working medium. Here, we show that by squeezing the working medium after the isochoric strokes, a significant amplification of the oscillation amplitude can be obtained. Using this method, temperature differences smaller than ±2 μK can be detectable in experiments.

We assume no prior knowledge of the bath temperatures. The squeezing operations described by the operator $\widehat{S}\left(\xi \right)=\mathrm{exp}\left(\frac{1}{2}\left({\xi }^{{\ast}}{\hat{a}}^{2}-\xi {\hat{a}}^{{\dagger}2}\right)\right)$ are applied to the radial state right after each of the two isochoric strokes, as such, the temperatures of the baths remain unchanged. After the squeezing operation, the state of the working medium is described by a squeezed thermal state, Ŝ(ξ)ρ_th Ŝ^†(ξ). In this case, as X₀ ≠ 0, Y₀ ≠ 0, R(t) becomes time-dependent (see equation (10)). The growth of the axial oscillation amplitude between two consecutive peaks reads

$$\begin{equation}{\Delta}{z}^{\prime }=\frac{4\kappa \gamma \hslash }{m{\omega }_{z}}\left(\mathrm{cosh}\left(2r\right)+\mathrm{sinh}\left(2r\right)\mathrm{cos}\left(\alpha \right)/\left(4{\kappa }^{2}-1\right)\right)\left({n}_{\mathrm{t}\mathrm{h}}^{1}-{n}_{\mathrm{t}\mathrm{h}}^{2}\right),\end{equation} \tag{ 20 }$$

with κ = ω_x0/ω_z , ${n}_{\mathrm{t}\mathrm{h}}^{1,2}={\left({\mathrm{e}}^{{\beta }_{1,2}\hslash {\omega }_{x0}}-1\right)}^{-1}$, r and α the amplitude and phase of the squeezing parameter following ξ = r e^iα . Note that Δz' is maximized when α = 0. Considering the limit of ℏω_x0 ≪ k_B T, we obtain

$$\begin{equation}{\Delta}{z}^{\prime }=\frac{4\gamma \hslash }{m{\omega }_{z}^{2}}\left(\mathrm{cosh}\left(2r\right)+\mathrm{sinh}\left(2r\right)/\left(4{\kappa }^{2}-1\right)\right){\Delta}T.\end{equation} \tag{ 21 }$$

To evaluate the impact of squeezing the working medium on the temperature difference resolution, we define the squeezing amplification factor

$$\begin{equation}\mathcal{A}=\frac{{\Delta}{z}^{\prime }}{{\Delta}z}=\mathrm{cosh}\left(2r\right)+\mathrm{sinh}\left(2r\right)/\left(4{\kappa }^{2}-1\right),\end{equation} \tag{ 22 }$$

where Δz denotes the growth of the amplitude without squeezing. The dependence of the amplification factor on the squeezing parameter is depicted in figure 5, where the system was numerically simulated for different squeezing parameters with T₁ = 0.11 mK, T₂ = 0.1 mK and κ = 10. The squeezing operations amplify the oscillation growth by an order of magnitude favoring the detection of small temperature differences. More details about the numerical simulations are presented in appendix B.

Zoom In Zoom Out Reset image size

The quantum regime is attained when the squeezing operation generates a state in which the fluctuations are reduced below the symmetric quantum limit in one of the quadrature components. This comes at the expense of increasing fluctuations in the canonical conjugate component, such that the uncertainty principle is still satisfied. In other words, the quantum limit is reached when the variance of one of the quadratures is smaller than 1/4, which implies that there is a phase in which the Glauber–Sudarshan distribution turns negative [43, 44]. This condition is given by

$$\begin{equation}\frac{1}{4}\left(2{n}_{\mathrm{t}\mathrm{h}}+1\right){\mathrm{e}}^{{\pm}2r}{< }\frac{1}{4}.\end{equation} \tag{ 23 }$$

In the example considered in figure 5, violation of equation (23) corresponds to a value of r > 0.77, as indicated by the vertical dashed line and the shaded region. The white region, on the other hand, points out that even when the fluctuations are not reduced below the symmetric quantum limit, the amplification factor is greater than one. However, to obtain a significant amplification one is forced to enter the quantum regime.

In trapped ion experiments, squeezing of the motional states can be realized by fast trap voltage control [45–47] or dynamic optical forces [48–50]. Both methods lead to a modulation of the motional state at twice the harmonic oscillation frequency. In a recent experimental demonstration of a quantum absorption refrigerator [50], squeezing operations with r up to 2 were realized with two detuned laser fields.

One should note that the optimal performance of a squeezed Otto cycle is obtained when squeezing is performed after the interaction of the working medium with the hot bath [28]. In this case, the amplitude growth is typically greater than applying squeezing after both isochoric strokes, enhancing the amplitude resolution of the flywheel. Although the simple linear relation Δz ∝ ΔT will not be satisfied, this scheme can be applied to evaluate the hot bath temperature given that the temperature of the cold bath is known.

We employ the Strömer–Verlet method [40] for propagating the classical oscillator in the axial direction. The classical oscillator is described by the Hamiltonian $H\left(z,v\right)=\frac{{v}^{2}}{2m}+{\Phi}\left(z\right)$, where Φ(z) denotes a potential which is solely dependent on the position of the particle. The canonical equations of motion can be written as z_n+1 − 2z_n + z_n+1 = −Δt(q/m)∂_z Φ(z) and v_n = (z_n+1 − x_n−1)/2Δt, with the time step Δt = t_n+1 − t_n . The following recursion relations are used [40]

$$\begin{equation}{v}_{n+\frac{1}{2}}={v}_{n}+\frac{{\Delta}t}{2}\frac{q}{m}{\partial }_{z}{\Phi}\left({z}_{n},{t}_{n}\right)\end{equation} \tag{ A.1 }$$

$$\begin{equation}{z}_{n+1}={z}_{n}+{\Delta}t{v}_{n+\frac{1}{2}}\end{equation} \tag{ A.2 }$$

$$\begin{equation}{z}_{n+1}={v}_{n+\frac{1}{2}}+\frac{{\Delta}t}{2}\frac{q}{m}{\partial }_{z}{\Phi}\left({z}_{n},{t}_{n}\right).\end{equation} \tag{ A.3 }$$

The time evolution of the radial quantum oscillator is described by the Liouville–von Neumann equation

$$\begin{equation}\mathcal{L}\left[\hat{\rho }\left(t\right)\right]=\mathrm{d}\hat{\rho }\left(t\right)/\mathrm{d}t=-\frac{\mathrm{i}}{\hslash }\left[\hat{H}\left(t\right),\hat{\rho }\left(t\right)\right]+\mathrm{i}{\mathcal{L}}_{D}\left[\hat{\rho }\right],\end{equation} \tag{ A.4 }$$

with the Liouvillian $\mathcal{L}$, the Hamiltonian Ĥ(t) and the Lindbladian ${\mathcal{L}}_{D}$. The Liouvillian is expanded in Newton polynomials in order to numerically integrate the Liouville–von Neuman equation. An arbitrary function f(x) with $x\in \mathbb{C}$ can be represented in terms of Newton polynomials R_n (x) with a set of sampling points $\left\{{x}_{\ell }\right\}$

$$\begin{equation}f\left(x\right)\approx \sum _{n=0}^{N-1}{a}_{n}{R}_{n}\left(x\right)\quad \text{with}\quad {R}_{n}\left(x\right)=\prod _{\ell =0}^{n-1}\left(x-{x}_{\ell }\right).\end{equation} \tag{ A.5 }$$

The coefficients a_n are computed by the recursion relation

$$\begin{equation}{a}_{n}=\frac{f\left({x}_{n}\right)-{\sum }_{\ell =0}^{n-1}{a}_{\ell }{\prod }_{m}^{\ell -1}\left({x}_{n}-{x}_{m}\right)}{\ell {\prod }_{m}^{\ell -1}\left({x}_{n}-{x}_{m}\right)},\end{equation} \tag{ A.6 }$$

which is referred to as divided differences [60]. The first two coefficients are obtained by imposing the interpolation condition

$$\begin{equation}{a}_{0}=f\left({x}_{0}\right)\quad \text{and}\quad {a}_{1}=f\left({x}_{1}\right)-f\left({x}_{0}\right).\end{equation} \tag{ A.7 }$$

With f(x) = e^−ixdt and $x=\mathcal{L}/\hslash$, the time evolution of the density matrix can be expressed as

$$\begin{equation}\hat{\rho }\left(\mathrm{d}t\right)={\mathrm{e}}^{-\frac{\mathrm{i}}{\hslash }\mathcal{L}\mathrm{d}t}\hat{\rho }\left(0\right)\approx \sum _{n}^{N-1}{a}_{n}\left(\mathcal{L}-{x}_{n}\mathbb{1}\right)\hat{\rho }\left(0\right).\end{equation} \tag{ A.8 }$$

For the simulations presented in this paper, we used the python implementation of the Newton propagator [39]. Details about the numerical methods can be found in reference [38].

Description of the engine dynamics requires a combination of the classical and quantum equations of motion. At every step of the propagation, the density matrix of the radial state, the position z and velocity v of the axial oscillator, and the force $\mathcal{F}$ needs to be updated. For each iteration, R(t) is computed and the force $\mathcal{F}$ is updated according to equation (9). The axial oscillator is then propagated by Δt using the Strömer–Verlet method. The newly obtained position z(t + Δt) is used to update the Hamiltonian of our system described by equation (3) and the Liouvillian is computed. The Newton propagator [39] is called and the state $\hat{\rho }\left(t\right)$ is propagated for the finite time step Δt. This finishes one step of the propagation. After $t=\frac{\pi }{{\omega }_{z}}$ the state is coupled to the hot or cold bath. We describe the thermalization process by updating the density matrix to the bath temperature in one step of propagation. This algorithm is illustrated in figure A1.

Zoom In Zoom Out Reset image size