A Feshbach engine in the Thomas-Fermi regime

Bose-Einstein condensates can be used to produce work by tuning the strength of the interparticle interactions with the help of Feshbach resonances. In inhomogeneous potentials, these interaction ramps change the volume of the trapped gas allowing one to create a thermodynamic cycle known as the Feshbach engine. However, in order to obtain a large power output, the engine strokes must be performed on a short timescale, which is in contrast with the fact that the efficiency of the engine is reduced by irreversible work if the strokes are done in a non-adiabatic fashion. Here we investigate how such an engine can be run in the Thomas-Fermi regime and present a shortcut to adiabaticity that minimizes the irreversible work and allows for efficient engine operation.


I. INTRODUCTION
Understanding and exploring concepts in quantum thermodynamics is currently a highly active topic with implications for the future development of quantum technologies [1].Within this area, the creation and operation of quantum engines which implement the Otto cycle and use cold atoms as their working medium, has received special attention, as they can be treated instructively and lend themselves to experimental realisation [2][3][4][5][6].
Similar to classical thermodynamical engines, quantum engines will achieve maximum efficiency if they are run without creating irreversible work.While this can be achieved by adiabatic evolution, this mode of operation has the drawback that it requires the external parameter changes to be very slow [7].As reliable, fast and simple control is needed e.g. for the development of new technologies [8], more recently the use of shortcuts to adiabaticity (STA) has received increased attention [9].Shortcut protocols provide a way to mimic adiabatic evolution in a finite time, mostly by requiring different parameter ramps or additional levels of control.While a large number of shortcuts have been found for single particle systems, shortcuts for interacting many-particle settings are still rare [10][11][12].However, first experiments have demonstrated the viability of shortcuts protocols for atomic Bose-Einstein condensates (BECs) in the repulsive mean-field limit [13,14] and for fermionic systems in the unitary limit [15].
Utilizing STA protocols to efficiently drive the dynamics through the control of the external potential parameters, such as trapping frequencies, can therefore increase the performance of finite time quantum heat engines [2,16].Furthermore, recent works have extended this idea to engines that are driven by changing internal parameters of the working medium, such as changes to the interparticle interaction strength [4,17].In ultracold atoms these are described by a scattering length, which can be controlled experimentally by varying an external magnetic field about a Feshbach resonance [18][19][20].While applying STA protocols to drive interactions is not a trivial task, Li et al. showed that it is possible in the case of a bright solitonic BEC which can be friction-lessly compressed and expanded using designed Feshbach pulses [4].Even though this can lead to an efficient Otto cycle, the operational range of the engine was very limited due to the possibility of BEC collapse in the presence of driven attractive interactions [21,22].
It is therefore interesting to extend the idea of the Feshbach engine to BECs in the stable Thomas-Fermi regime of large particle numbers and strong repulsive interactions [23].For this we derive in this work a novel interaction ramp that allows for the frictionless compression and expansion of such a Thomas-Fermi BEC in an almost arbitrarily short time, and show that these ramps can act as STAs in a Feshbach engine.The presentation is organized as follows: In Sec.II we introduce a scaling ansatz to derive an interaction ramp for a harmonically trapped d-dimensional BEC in the Thomas-Fermi limit which ensures that the system follows an adiabatic path at all times.We then verify that the ramp is working as intended, up to some minimum time, by numerically simulating the dynamics using the full Gross-Pitaveskii equation and comparing it to a non-optimized reference ramp.In Sec.III we show that the shortcut ramp can indeed be used to increase the power and efficiency of the engine and in Sec.IV we perform a stability analysis to derive the minimum time in which the interaction ramp can be performed before a modulational instability leads to a condensate collapse.

II. SHORTCUT TO ADIABATICITY
In this section we derive an interaction ramp for a BEC in the Thomas-Fermi limit that can act as an STA [9,24] and evaluate its performance in compressing and expanding a BEC compared to a smooth non-optimized reference ramp.

A. Interaction ramp
For simplicity we start by considering a onedimensional BEC trapped in an inhomogeneous potential, V (x), and generalise the results to higher dimen-arXiv:2005.06801v1[quant-ph] 14 May 2020 0 0.2 0.4 0.6 0.8 1 0.8 0.9 1. Interaction ramps obtained from the shortcut to adiabaticity according to Eq. ( 14) in three dimensions and for different values of T f .The black dashed line shows the timerescaled adiabatic reference (TRA).
sions later.The condensate can be described by a single wave function, ψ(x), whose dynamics is governed by the Gross-Pitaevskii equation [25] where m is the mass of an individual particle in the condensate, g(t) is the nonlinear interaction strength, which may vary in time, and the wave function is normalized to the number of particles, dx|ψ(x)| 2 = N .Since we are interested in compressing or expanding the width of the wave function without changing its general shape, we choose a scaling ansatz [26] of the form ψ(x, t) = 1 a(t) e iϕ(x,t) φ(y(x, t), τ (t)) , where the spatial coordinate is rescaled as y(x, t) = x/a(t) and we have also introduced a rescaled time τ (t).Inserting this ansatz into the Gross-Pitaevskii equation (1) and choosing the phase as to eliminate the term proportional to ∂φ/∂y, we get The choice of phase made in Eq. (3) can be interpreted as a gauge transformation Û = e iϕ(x,t) , which adjusts the momentum of the expanding or shrinking system as where ȧx/a is the local velocity [27,28].The same transformation is also commonly found in other scaling problems, e.g. as the optimal solution in a variational approach [29][30][31].Assuming the external potential is given by a harmonic trap, V (x) = 1 2 mω 2 x 2 , we obtain Choosing the rescaled time τ and the term ä + ω 2 a in such a way that it leads to a solvable Gross-Pitaevskii equation then allows one to design control pulses for a frictionless evolution of the BEC, e.g. by varying either the trap frequency ω, the interaction strength g or both [28,[32][33][34].
Many experiments involving repulsively interacting Bose-Einstein condensates are carried out in the so-called Thomas-Fermi (TF) regime, where the potential and the interaction energies are much larger than the kinetic energy, N g 3 ω/m [25].This allows one to neglect the kinetic energy term in the Gross-Pitaevskii equation ( 1) and obtain an analytical solution of the form and ψ(x, t) ≡ 0 otherwise.The chemical potential µ is determined via the normalization condition.Considering the TF limit in the scaling GPE Eq. ( 6) and choosing scaling functions as [35] ä + ω 2 a = ω 2 g(t) gi with some initial interaction strength g i , then gives This again has the aforementioned Thomas-Fermi solution with . Inserting everything back into the scaling ansatz gives us an analytic expression for the evolution of the wave function and choosing a suitable scaling function a(t) then allows us to reverse-engineer an interaction ramp that will take the system along this evolution.For this we rearrange Eq. ( 8) for g(t) and get Choosing appropriate boundary conditions for a(t) of the form we can drive the system from the Thomas-Fermi ground state at an initial interaction strength g i to the ground state at a final value g f in an almost arbitrarily short time T f , while mimicking an adiabatic evolution.It is worth noting that by doing this the system also acquires an additional, but irrelevant phase that depends on T f .The boundary conditions for a(t) can be fulfilled by a fifthorder polynomial a(t) = a i +(a f −a i ) 10s 3 − 15s 4 + 6s 5 with s = t/T f , which has the form of a smoother stepfunction [36].This shortcut to adiabaticity is easily generalised to a d-dimensional BEC in an isotropic harmonic trap in the Thomas-Fermi limit by driving the interaction strength according to and requiring a(t F ) = (g f /g i ) 1/(d+2) as well as replacing g i a(t ) with g i a d (t ) in the time scaling τ in Eq. ( 8).In the following we will use dimensionless units and scale lengths by x 0 = /mω, energies by ω, time by ω and interaction strengths by ωx d 0 .

B. Evaluation
As a reference for benchmarking the shortcut performance we define a time-rescaled adiabatic (TRA) stroke, which can be obtained by letting T f → ∞ or equivalently by setting ä(t) ≡ 0 in the interaction ramp g(t).A comparison between the TRA ramp for compressing a threedimensional Thomas-Fermi BEC from g i = 1 to g f = 0.8 with STA ramps for varying T f is shown in Fig. 1, and one can immediately notice that faster ramps require larger changes in the interaction strength over the duration of the STA.We can assess the performance of the shortcut by numerically simulating the Gross-Pitaevskii equation with the calculated interaction ramps and evaluating the irreversible work and fidelity, at the end of the stroke as a measure of how close the system's state |ψ(T f ) with energy E(T f ) after the evolution is to the desired target state |ψ target with energy E f .For a three-dimensional BEC consisting of N = 10 4 atoms and a compression going from g i = 1 to g f = 0.8 we show these quantities for the STA and TRA strokes in the upper row of Fig. 2. For comparison we also show them for the reverse expansion stroke going from g i = 0.8 to g f = 1 in the lower row of Fig. 2, but for N = 8000 instead.We note we use the actual ground states for the full Gross-Pitaevskii equation as the initial and target states for the evolution, which differ slightly from the Thomas-Fermi wave function in Eq. ( 7), especially around the condensate edges.Therefore their energy is slightly higher than the Thomas-Fermi value of Nevertheless, one can see that in each case the STA outperforms the TRA stroke by several orders of magnitude in the irreversible work for almost any stroke duration T f above some threshold T min f .Similarly, above this threshold the shortcut always achieves a perfect fidelity of F = 1 with the target state while the TRA falls short.The sharp dips in the irreversible work for some stroke times as well as the near-perfect fidelity for the TRA around T f ≈ π are accidental and can be attributed to the underlying dynamics in the harmonic trapping potential [37].
The sudden increase in irreversible work and the accompanying drop of the fidelity to basically zero if the shortcut is performed too fast is due to a modulational instability that exists for the nonlinear Schrödinger equation (see for example [38]).It is triggered by the shortcut ramp driving the BEC at attractive interactions for extended periods of time, resulting in a collapse of the condensate.This collapse can create trains of bright solitons [39][40][41] which lead to a clear deviation from the adiabatic path.In the particular example considered here and in the next section, the threshold below which the modulational instability appears is roughly T min f ≈ 0.05 and in 0.7 0.75 0.8 0.85 0.9 0.95 1 1.05 1.1 3. The energy of the system in the Thomas-Fermi regime according to Eq. ( 16), as a function of the interaction strength g for N = 10  Section IV we present a more detailed stability analysis to derive a general criterion for a given dimensionality, chemical potential, change in interaction and stroke time.

III. FESHBACH ENGINE
A Feshbach engine operates an Otto cycle which consists of two adiabatic interaction ramps that compress and expand a BEC, and which are connected by two iso-choric strokes that add or remove particles [4].The latter can be realized by cooling or heating the thermal cloud of atoms surrounding any BEC and thus prompting atoms to condense into or evaporate from the condensate, respectively.
In the following we will evaluate the Otto engine cycle driven in a harmonically trapped, three-dimensional BEC, where the interaction strength goes between g i = 1 and g f = 0.8 for the adiabatic strokes and the particle number between N i = 10 4 and N f = 8000 for the isochoric strokes (see Fig. 3).We quantify the performance of the engine by calculating its efficiency and power where the compression and expansion work always calculated from the actual values obtained after performing the strokes and not from the adiabatic values indicated in Fig. 3.For the cycle duration τ we assume that the isochoric strokes can be performed a lot faster than the work strokes and that we can therefore approximate it as τ ≈ 2T f with the work stroke duration T f .An analytical expression for the maximally attainable adiabatic efficiency of such a Thomas-Fermi Feshbach engine can be obtained by using the Thomas-Fermi wave functions to calculate the energy of the BEC at the engine cycle end points and is given as a function of the compression ratio g f /g i by with γ = 2/3 (1D BEC), γ = 1/2 (2D BEC) and γ = 2/5 (3D BEC).It is worth noting that for the one-dimensional case this is less efficient than the bright soliton Feshbach engine using the same compression ratio, which has an exponent of γ = 2 [4].
The efficiency and power of the engine cycle, using the adiabatic strokes shown in Fig. 2, are plotted in Fig. 4(a) and (b) respectively.One can immediately see that the STA enables the engine to reach its maximum adiabatic efficiency of η AD = 0.0854 already for cycle times four to five times shorter than in the TRA case and leads to a considerable increase in engine power for all cycle times considered, as long as the modulational instability is not triggered.Plotting the ratio P STA /P TRA in Fig. 4(d) shows that this increase can reach up to 60%.As expected, the advantage decreases and the STA and TRA perform equally well once the cycle times are increased to more and more adiabatic values.Finally, plotting the engine efficiency versus power in Fig. 4(c) shows that the shortcut enables the engine to run with high efficiency even at high power output, which is an important factor in the operation of any heat engine [42].
STA -T f = 0.5 STA -T f = 0.4 TRA FIG. 5. Interaction ramp given by the shortcut to adiabaticity for reducing the interactiong strength of a one-dimensional BEC from gi = 1 to g f = 0.8 in time T f according to Eq. ( 14).The time-rescaled adiabatic reference (TRA) obtained by setting ä ≡ 0 is used for comparison.For these parameters, the modulational instability seems to be triggered once the ramp's minimum goes below −gi (dashed line), which happens roughly around T f ≈ 0.45.x FIG. 6. Emergence of the modulational instability in the condensate density for compressing a one-dimensional BEC with N = 10 4 atom in the Thomas-Fermi regime from gi = 1 to g f = 0.8 by the interaction STA in Eq. ( 14) in a time T f = 0.4.

IV. MODULATIONAL INSTABILITY
To better understand the limit of the TF Feshbach engine, we perform in the following a stability analysis to determine the threshold times, T min f , below which the shortcuts to adiabaticity derived in Sec.II fail due to the appearance of a modulational instability.These threshold times act as an intrinsic quantum speed limit [43] for the manipulation of the Thomas-Fermi BECs considered here.For simplicity, we will again first perform the analysis for the compression of a one-dimensional BEC and show later that the obtained stability criterion can easily be extended to higher dimensions.This is confirmed by comparison with numerical simulations.
The general reason for an instability to occur is that for short manipulation times T f , the system needs to be driven into the regime of attractive interactions, g < 0, for a certain amount of time (see Fig. 5), and entering this regime too deeply will lead the condensate to collapse [44].Using the same parameters of N = 10 4 , g i = 1 and g f = 0.8 as in the previous sections, we find that the modulational instability for compressing a one-dimensional BEC via the shortcut is triggered around T min f ≈ 0.45.To understand this instability, let us first carefully examine its appearance in the Gross-Pitaevskii dynamics.
From the density distributions shown in Fig. 6 one can see that the instability manifests itself by the appearance of rapid oscillations that develop in the center of the condensate where the density is maximal.Inserting the analytical solution for the wave function dynamics from Eq. ( 11) into the GPE gives from which, or by directly setting x = 0 in Eq. ( 11), one can see that close to the trap center, where the density is approximately homogeneous and the kinetic and potential energies can be neglected, the wave function evolves according to where µ(t) = µ i (äa + a 2 ).It is worth noting that from Eq. ( 19) one can see that the interaction ramp accomplishes the condensate rescaling by creating a timevarying harmonic trapping potential with a time-varying ground state energy.We assume that the modulational instability can be described by a perturbation to the homogeneous solution of the form of ψ(x, t) = ψ hom (t) [1 + u(t) cos(kx)], with complex amplitude u(t) = u re (t) + iu im (t) and wave number k. Inserting this into the Gross-Pitaevskii equation and keeping only terms up to first order in u then leads to which is similar to Hill's equation [45].Note that the same equation can be obtained for the 2D and 3D case by replacing µ i and a f in µ(t) with the appropriate values.
Let us note that if µ(t) was periodic, i.e. if we had choosen a sinusoidal function to fulfil the boundary conditions (13) [46], Floquet theory would provide exact conditions for the stability of Eq. ( 22) [45,47].However, a ramp resulting from such an approach generally has maxima and minima with larger magnitude than the corresponding polynomial ramp and therefore triggers the instability even earlier.Furthermore, the instability would usually occur around t ≈ T f /4, when the ramp reaches its minimum, which means that the periodicity of the ramp would not come into play.
Since we are considering a complex amplitude u(t), it is helpful to do a change of coordinates.For this we rewrite 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.5 7. Minimum time for the compression of a onedimensional BEC in the Thomas-Fermi regime consisting of N = 10 4 atoms by changing its interaction from gi = 1 (blue curve) or gi = 2 (red curve) to a final value of g f via the STA in Eq. ( 14).The curves represent the lines ∆ = 10 14 for the stability criterion (27).The yellow curve shows T min f for a three-dimensional BEC with N = 10 4 , gi = 1 and the same ∆ = 10 14 .The circles show numerically obtained values for T min f from simulating the full Gross-Pitaevskii equation in each case and using the sharp increase in Wirr as an indicator for the instability.
Eq. ( 21) in polar coordinates using u(t) = r(t)e iϕ(t) and omitting the terms proportional to ȧ/2a to get In this form it is easy to see that for µ(t) < 0 and ϕ = arccos and where we have defined This allows us to define the total relative increase in the perturbation's amplitude after the interaction ramps as ∆ = exp The previous observation that the modulational instability starts forming at the center of the condensate indicates that it is triggered or seeded by noise [48,49].In our numerical simulations this is numerical noise, which appears on the level of 10 −14 .Since it grows exponentially, a good choice as a criterion for stability is to set ∆ = 10 14 , i.e. the point at which the small perturbation on the BEC's wavefunction becomes comparable to magnitude of the BEC wavefunction itself.The resulting curves for a compression stroke are shown in Fig. 7.
The minimal compression times T min f from an initial interaction g i to some final value g f < g i agree well with the actual times obtained from numerically simulating the GPE both for one-dimensional and three-dimensional BECs.While we found ∆ = 10 14 to be the best fit for our numerical data, changing the criterion even by several orders of magnitude only leads to slight deviations in the resulting stability curves.

V. CONCLUSION
By using a scaling ansatz we have exactly solved the dynamics of an atomic Bose-Einstein condensate subject to a specific interaction ramp in the repulsive Thomas-Fermi limit.This interaction ramp provides a shortcut to adiabaticity for driving the condensate from one interaction strength to another and thereby compressing or expanding it in a short amount of time while avoiding unwanted excitations.We have shown how these shortcuts can increase the efficiency of a Feshbach engine by using them as the adiabatic strokes of its Otto cycle.Using numerical simulations of the full condensate dynamics, we have shown that the speed-up of the condensate manipulation is limited by a modulational instability leading to a condensate collapse.The instability is caused by the need to drive the condensate at increasingly attractive interactions for longer periods of time as the ramp becomes shorter and shorter.Finally, we have performed a stability analysis and determined a criterion that provides an accurate limit T min f for a given initial chemical potential and final interaction strength.

FIG. 2 .
FIG. 2. Upper row:Irreversible work Wirr and fidelity F after compressing a 3D BEC consisting of N = 10 4 atoms from an initial interaction gi = 1 to g f = 0.8 in a time T f using the shortcut to adiabaticity (STA) and an adiabatic reference protocol (TRA).Lower row: Identical plots but for the reverse expansion stroke and N = 8000 atoms.

FIG. 4 .
FIG. 4. (a)Efficiency η and (b) power P of the threedimensional Feshbach engine as a function of the engine cycle duration τ = 2T f , comparing the performance of the shortcut to adiabaticity (STA) for the adiabatic strokes with an adiabatic reference (TRA).The maximally attainable efficiency for this specific engine cycle is ηAD = 0.0854.(c) Plot of efficiency against engine power for both the STA and TRA adiabatic strokes and (d) relative increase in engine power by using the STA strokes as a function of cycle duration τ .

k 2 √ 4 , ( 24 )
|µ(t)|we have φ = 0 and the amplitude of the perturbation can increase exponentially according toṙ = rk |µ(t)| − k 2 if k < 2 |µ(t)|.The increase is maximal for k = 2|µ(t)|, for which we have ṙ = μ(t)r or r(t) = e t 0 dτ μ(τ ) r(0) , 4(upper dashed line) and N = 8000 (lower dashed line).The solid lines indicate the Otto engine cycle consisting of two adiabatic strokes between gi = 1 and g f = 0.8, performing compression and expansion work, WC and WE respectively, and two isochoric strokes adding or removing heat QN + or QN − respectively by adding or removing particles to and from the condensate.