Introduction

The measurement’s outcome of a phase-sensitive sensor probing forces exerted on neutral atoms by inertial, material or electromagnetic sources depends dramatically on the initial conditions, i.e. on the position, velocity and size of the input matter-wave. A lack of knowledge or scattering in these initial properties inevitably leads to systematic effects or statistical errors harming the sensor’s performance. An example of the degree of control needed can be grasped if one considers making a test of the Universality of Free Fall (UFF) with two different atomic species to put bounds on a possible violation of the UFF at the femto-level in the Eötvös ratio1, level at which state-of-the-art experiments perform with material test masses2. Such a precise experiment requires that the initial positions, center-of-mass velocities and expansion rates are defined at a level better than 1 μm, 1 μm/s and 100 μm/s (35 pK in 3D), respectively3.

To meet these stringent requirements, the temperatures of the atomic ensembles have to be drastically reduced (down to a sub-nK level) and their size must remain compact (not exceeding a few mm after several seconds of free expansion) clearly indicating the necessity of using Bose-Einstein Condensates (BEC). Such a direction is taken by several metrology groups worldwide4,5,6,7,8,9,10, including the QUANTUS11 and MAIUS12 consortia, which reached important milestones in controlling quantum gases dynamics in microgravity conditions using atom chips13,14.

In a recent work15, we considered an approach based on Shortcut-To-Adiabaticity (STA) protocols to obtain analytic solutions for the transport of the BEC in an atom chip setup with realistic anharmonic and rotating trapping potentials. This approach based on the reverse engineering technique allows for a full control of the translational degrees of freedom of the BEC. It is, however, exciting several collective modes of the quantum gas, an effect which could eventually compromise the expected metrological gain if such a source is used without any precaution as an input of an atom interferometer. It is in this context that the use of optimal control theory (OCT) can reveal an unchallenged potential of targeting a given final state in timescales shorter than the trivial adiabatic manipulation, which is of no practical use in the metrology context since it is associated with poor cycling rates.

The aim of optimal control theory is to bring a dynamical system from one state to another, while minimizing a cost functional, such as the control time or the energy of the pulse used. The modern version of OCT is born with the Pontryagin’s Maximum Principle (PMP) in the late 1950s16,17. Originally applied to problems of space dynamics, OCT is nowadays a key tool to study a large spectrum of applications both in classical18,19 and quantum physics20,21,22. In the Pontryagin formulation, solving an optimal control problem is equivalent to finding extremal trajectories which are solutions of a generalized Hamiltonian system. These trajectories satisfy the maximization condition of the PMP as well as specific boundary conditions18,19,20. The implementation of the PMP is far from being trivial and numerical control algorithms have been developed to approximate the optimal solution23. Among others, we can mention the gradient19,24 and the Krotov23,25 algorithms, which are nowadays standard tools in physics.

OCT has been applied with success to quantum systems since the 1980s in domains extending from molecular physics and nuclear magnetic resonance to quantum information science (see refs22,26, for recent reviews, and references therein). The application of OCT to BEC dynamics has also been explored in different contexts. Using the Gross-Pitaevskii equation, the optimal coherent manipulation of an atomic BEC has been investigated in a series of studies (see refs27,28,29,30,31,32,33 to cite a few, and references therein). The transport of cold atoms has also been optimized for simple models in combination with invariant-based inverse methods15,34,35,36,37. It should be mentioned here that OCT and STA are usually compatible in the sense that an OCT methodology can be built on top of a basic STA frame of solutions34,38,39,40,41,42. One can also note that recently, new methods have been tested successfully to bridge the gap between an ideal STA and a realistic experimental implementation for the optical transfer of a degenerate gas, demonstrating fast highly non-adiabatic transfer with almost no residual sloshing using corrected STA trajectories43.

In this paper, we discuss the application of optimal control theory for the fast transport of Bose-Einstein condensates with atom chips while simultaneously controlling the quantum degrees of freedom of the problem to target the ground state of the final trap as the optimization result. The article is organized as follows: We describe in the next section the chip model used to transport the BEC. The chosen cost functional and the associated transport ramp are presented in the next section, which is followed by a comparison of our findings to the results of the STA technique applied in a similar context15. We finally illustrate the impact of the OCT ramp duration on the internal degrees of freedom of the final BEC state. We conclude by discussing the limits of the methodology we have developed, and by mentioning potential experimental implementations.

Chip model

We consider the case of a Z-shaped chip configuration used to trap and manipulate cold Rb atoms in micro-gravity (See reference15 for a detailed description of the numerical model and reference44 for the description of an experimental implementation). The three spatial dimensions are denoted by the three coordinates X, Y and z. z is the direction perpendicular to the chip. X and Y are two orthogonal directions in the plane of the chip. The trap is naturally rotating in the (XY) plane when the physical parameters governing the trap potential change. The diagonalization of the associated Hessian matrix allows to define two new eigen-coordinates x and y of the trap, rotated compared to the fixed X and Y coordinates15. The physical parameters which govern the trap potential are the chip intensity Iw and the bias magnetic field Bbias. For the present study, Iw is fixed at 5 A and the control parameter for the implementation of the transport ramp is the time-dependent bias magnetic field Bbias(t), which varies between Bbias(0) = Bi = 21.5 G at the initial time t = 0 and Bbias(tf) = Bf = 4.5 G at the end of the transport corresponding to t = tf  .

In such a configuration already described in our previous study15, the minimum of the trap is at the origin in x and y, and it is located at a distance z0(t) from the chip surface. At t = 0 we have \({z}_{0}(0)\simeq 0.45\) mm and at the end of the transport \({z}_{0}({t}_{f})\simeq 1.65\) mm. If we limit ourselves, in a first approximation, to the simplest case of a time-dependent harmonic trap, the center-of-mass of the condensate zA(t) in the direction normal to the surface follows Newton’s equations of motion

$${\dot{z}}_{A}(t)={v}_{A}(t)\,{\rm{and}}\,{\dot{v}}_{A}(t)=-\,{\omega }_{z}^{2}(t)[{z}_{A}(t)-{z}_{0}(t)],$$
(1)

where ωz(t) denotes the frequency of the trapping potential along z at time t. Moreover, in the Thomas-Fermi approximation45, the evolution of the size of the BEC is accurately described by a scaling approach46,47. The size of the BEC is defined by the three time-dependent radii rx(t), ry(t) and rz(t) of the paraboloid associated with the BEC wave function, with

$${r}_{x}(t)={r}_{x}(0)\,{\lambda }_{x}(t),\,{r}_{y}(t)={r}_{y}(0)\,{\lambda }_{y}(t)\,{\rm{and}}\,{r}_{z}(t)={r}_{z}(0)\,{\lambda }_{z}(t).$$
(2)

It was shown46,47 that the time-dependent scaling factors λx(t), λy(t) and λz(t) obey the three coupled second order differential equations

$$\begin{array}{ccc}{\ddot{\lambda }}_{x} & = & \frac{{\omega }_{x}^{2}(0)}{{\lambda }_{x}^{2}{\lambda }_{y}{\lambda }_{z}}-{\omega }_{x}^{2}(t)\,{\lambda }_{x},\\ {\ddot{\lambda }}_{y} & = & \frac{{\omega }_{y}^{2}(0)}{{\lambda }_{x}{\lambda }_{y}^{2}{\lambda }_{z}}-{\omega }_{y}^{2}(t)\,{\lambda }_{y},\\ {\ddot{\lambda }}_{z} & = & \frac{{\omega }_{z}^{2}(0)}{{\lambda }_{x}{\lambda }_{y}{\lambda }_{z}^{2}}-{\omega }_{z}^{2}(t)\,{\lambda }_{z},\end{array}$$
(3)

where ωx(t) and ωy(t) denote the frequencies of the trapping potential along x and y at time t. The full behavior of the trapping frequencies as a function of the control parameter Bbias can be found in15. Initially the trapping frequencies are \({\omega }_{x}\mathrm{(0)}\simeq 2\pi \cdot 15\,{\rm{Hz}}\) and \({\omega }_{y}(0)\simeq {\omega }_{z}(0)\simeq 2\pi \cdot 616\,{\rm{Hz}}\). At the end of the transport \({\omega }_{x}({t}_{f})\simeq 2\pi \cdot 10\,{\rm{Hz}}\) and \({\omega }_{y}({t}_{f})\simeq {\omega }_{z}({t}_{f})\simeq 2\pi \cdot 32\,{\rm{Hz}}\). The largest time scale associated with the trap is therefore of the order of 100 ms. An adiabatic transport would thus require transport durations larger than 1 s. Here we want to design a simple, fast and efficient transport ramp for the BEC. The OCT technique being very powerful, we have decided to optimize a single control parameter, Bbias(t), in order to control the final position of the BEC zA(tf), its final speed vA(tf), and its final size defined by the three final scaling factors λx(tf), λy(tf) and λz(tf). We also wish to control the final expansion rates given by \({\dot{\lambda }}_{x}({t}_{f}\,)\), \({\dot{\lambda }}_{y}({t}_{f}\,)\) and \({\dot{\lambda }}_{z}({t}_{f}\,)\). Finally, since we want the harmonic approximation to hold during the entire transport, we also wish to limit the time-dependent offset between the position of the center of mass of the BEC and the center of the trap |zA(t) − z0(t)| as well as the the time-dependent offset between their respective speeds \(|{v}_{A}(t)-{\dot{z}}_{0}(t)|\). To be compatible with metrology applications with an integration over thousands of experimental cycles, we want this transport to be realized quickly, i.e. in a duration of the order of the largest time scale associated with the trap, that is of the order of 100 ms with the present chip configuration.

Cost functional

To implement such an optimal control scheme, we first introduce the “classical” point-wise translational energy of the condensate in the reference frame of the trap

$${E}_{cl}(t)=\frac{m}{2}({\omega }_{z}^{2}{[{z}_{A}-{z}_{0}]}^{2}+{[{v}_{A}-{\dot{z}}_{0}]}^{2}),$$
(4)

as well as the “quantum” energy associated with the 3D Thomas-Fermi wave function

$${E}_{qu}(t)=\frac{m}{14}[{\omega }_{x}^{2}{r}_{x}^{2}+{\omega }_{y}^{2}{r}_{y}^{2}+{\omega }_{z}^{2}{r}_{z}^{2}]+\frac{m}{14}[{\dot{r}}_{x}^{2}+{\dot{r}}_{y}^{2}+{\dot{r}}_{z}^{2}]+\frac{15gN}{28\pi {r}_{x}{r}_{y}{r}_{z}},$$
(5)

where \(g=4\pi {\hslash }^{2}{a}_{s}/m\) is the scattering amplitude, as is the s-wave scattering length of Rb-87 and N = 105 denotes the number of condensed atoms. The first term in Eq. (5) describes the potential energy associated with the finite size of the condensate, the second term is the kinetic energy associated with the size dynamics, and the third and last term is the average mean-field interaction energy between the atoms of the condensate. The numerical factors (1/14) and (15/28) seen in Eq. (5) come from the specific definition given in Eq. (2) of the size of the condensate using a Thomas-Fermi expression for the probability density.

The goal we want to achieve is the minimization of a total cost functional Ctot, defined by the sum

$${C}_{tot}={C}_{term}+{C}_{run}$$
(6)

of a terminal cost

$${C}_{term}={\lambda }_{1}\,{E}_{cl}({t}_{f})+{\lambda }_{2}\,{E}_{qu}({t}_{f})$$
(7)

and a running cost

$${C}_{run}={\lambda }_{3}\,(\frac{1}{{t}_{f}}{\int }_{0}^{{t}_{f}}{E}_{cl}(t)\,dt).$$
(8)

The terminal cost was designed to insure the formation of the ground state of the trap at time tf  . It imposes the minimization of the total energy of the condensate at the end of the transport. The running cost is introduced in order to limit the transient excitation of the condensate in the moving harmonic trap. Here we fix λ1 = 1 and the two other dimensionless parameters λ2 and λ3 are chosen to express the relative weights between the three terms of the cost functional. Changing the values of λ2 and λ3 affects the progress of the optimization procedure by changing the path it takes during optimization. This can lead in practice to different final transport ramps, which will take into account the relative weight assigned to each of the terms of the cost functional.

Transport ramp

The initial and final traps are defined by the initial and final values Bi = 21.5 G and Bf = 4.5 G of the bias magnetic field Bbias(t). Since in experiments one can be limited by the switch on/off speed of the magnetic field we circumvent this problem by insuring a smooth variation of Bbias(t) at t = 0 and at t = tf  . For this reason we have chosen to impose

$$B(t)={B}_{i}+({B}_{f}-{B}_{i})(10\,{[\frac{u(t)-{u}_{0}}{{u}_{f}-{u}_{0}}]}^{3}-15{[\frac{u(t)-{u}_{0}}{{u}_{f}-{u}_{0}}]}^{4}+6{[\frac{u(t)-{u}_{0}}{{u}_{f}-{u}_{0}}]}^{5}),$$
(9)

where u(t) is a continuous function of time, with u0 = u(0) and uf = u(tf). This definition allows to impose the following boundary conditions for the bias magnetic field

$$\begin{array}{lll}{B}_{bias}(0)={B}_{i}, & & {B}_{bias}({t}_{f})={B}_{f},\\ {\dot{B}}_{bias}\mathrm{(0)}=0, & {\rm{and}} & {\dot{B}}_{bias}({t}_{f})=0,\\ {\ddot{B}}_{bias}\mathrm{(0)}=0, & & {\ddot{B}}_{bias}({t}_{f})=0\,.\end{array}$$
(10)

Note that a consequence of these boundary conditions imposed on Bbias(t) is that similar relations hold for all trap parameters such has the trap position z0(t) and the trap frequencies in all directions ωx(t), ωy(t) and ωz(t). The optimization procedure we have adopted is therefore using the dimensionless control function u(t), from which we can calculate the optimal bias magnetic field using Eq. (9).

Optimal control

We now reformulate our optimization problem in the framework of optimal control theory. We refer the interested reader to standard textbooks for details18,19,20,21. The state of the system is described by a state vector x, with

$$\begin{array}{llll}{x}_{1}={z}_{A}(t), & {x}_{3}={\lambda }_{x}(t), & {x}_{5}={\lambda }_{y}(t), & {x}_{7}={\lambda }_{z}(t)\\ {x}_{2}={v}_{A}(t), & {x}_{4}={\dot{\lambda }}_{x}(t), & {x}_{6}={\dot{\lambda }}_{y}(t), & {x}_{8}={\dot{\lambda }}_{z}(t)\end{array}$$
(11)

As suggested by Eqs (1) and (3), the time evolution of all components of the state vector x is governed by a set of coupled first order differential equations controlled by u(t) through the time dependence of the trap position and frequencies. Once u(t) is chosen and for well defined initial conditions at t = 0, these equations are easily solved using a Runge-Kutta algorithm48,49 or the Verlet method50, for instance.

According to the Pontryagin maximum principle16,17, the extremal solutions of the problem, candidates to be optimal, satisfy the equations of Hamiltonian

$${\dot{x}}_{i}=+(\frac{\partial {H}_{p}}{\partial {p}_{i}})\,{\rm{and}}\,{\dot{p}}_{i}=-(\frac{\partial {H}_{p}}{\partial {x}_{i}}),\,{\rm{for}}\,i=1,2,\ldots ,8$$
(12)

where p is the adjoint state vector and where the Pontryagin Hamiltonian of the system is defined by

$${H}_{p}({\bf{x}},{\bf{p}},t,u)={\bf{p}}\cdot \dot{{\bf{x}}}-{\lambda }_{3}(\frac{{E}_{cl}(t)}{{t}_{f}}).$$
(13)

From Eq. (12), it can be easily shown that the dynamics of the adjoint state is governed by the following set of coupled first order differential equations

$${\dot{p}}_{1}={\omega }_{z}^{2}(t)[{p}_{2}+{\lambda }_{3}\frac{m}{{t}_{f}}({x}_{1}-{z}_{0})]$$
(14a)
$${\dot{p}}_{2}=-\,{p}_{1}+{\lambda }_{3}\frac{m}{{t}_{f}}({x}_{2}-{\dot{z}}_{0})$$
(14b)
$${\dot{p}}_{3}={p}_{4}[{\omega }_{x}^{2}(t)+\frac{2{\omega }_{x}^{2}\mathrm{(0)}}{{x}_{3}^{3}{x}_{5}{x}_{7}}]+{p}_{6}[\frac{{\omega }_{y}^{2}(0)}{{x}_{3}^{2}{x}_{5}^{2}{x}_{7}}]+{p}_{8}[\frac{{\omega }_{z}^{2}(0)}{{x}_{3}^{2}{x}_{5}{x}_{7}^{2}}]$$
(14c)
$${\dot{p}}_{4}=-\,{p}_{3}$$
(14d)
$${\dot{p}}_{5}={p}_{4}[\frac{{\omega }_{x}^{2}(0)}{{x}_{3}^{2}{x}_{5}^{2}{x}_{7}}]+{p}_{6}[{\omega }_{y}^{2}(t)+\frac{2{\omega }_{y}^{2}(0)}{{x}_{3}{x}_{5}^{3}{x}_{7}}]+{p}_{8}[\frac{{\omega }_{z}^{2}(0)}{{x}_{3}{x}_{5}^{2}{x}_{7}^{2}}]$$
(14e)
$${\dot{p}}_{6}=-\,{p}_{5}$$
(14f)
$${\dot{p}}_{7}={p}_{4}[\frac{{\omega }_{x}^{2}\mathrm{(0)}}{{x}_{3}^{2}{x}_{5}{x}_{7}^{2}}]+{p}_{6}[\frac{{\omega }_{y}^{2}\mathrm{(0)}}{{x}_{3}{x}_{5}^{2}{x}_{7}^{2}}]+{p}_{8}[{\omega }_{z}^{2}(t)+\frac{2\,{\omega }_{z}^{2}\mathrm{(0)}}{{x}_{3}{x}_{5}{x}_{7}^{3}}]$$
(14g)
$${\dot{p}}_{8}=-\,{p}_{7}.$$
(14h)

In addition, the transversality conditions for the adjoint state read

$${p}_{n}({t}_{f})=-\,{\lambda }_{1}{(\frac{\partial {E}_{cl}}{\partial {x}_{n}})}_{t={t}_{f}}-{\lambda }_{2}{(\frac{\partial {E}_{qu}}{\partial {x}_{n}})}_{t={t}_{f}}$$
(15)

thus leading to the following boundary conditions at time t = tf

$${p}_{1}({t}_{f})=-\,{\lambda }_{1}\,m{\omega }_{z}^{2}({t}_{f})\,[{z}_{A}({t}_{f})-{z}_{0}({t}_{f})]$$
(16a)
$${p}_{2}({t}_{f})=-\,{\lambda }_{1}\,m{v}_{A}({t}_{f})$$
(16b)
$${p}_{3}({t}_{f})=-\,{\lambda }_{2}[\frac{m}{7}{\omega }_{x}^{2}({t}_{f}){r}_{x}(0){r}_{x}({t}_{f})-\frac{15gN{r}_{x}(0)}{28\pi \,{r}_{x}^{2}({t}_{f}){r}_{y}({t}_{f}){r}_{z}({t}_{f})}]$$
(16c)
$${p}_{4}({t}_{f})=-\,{\lambda }_{2}\frac{m}{7}{r}_{x}(0){\dot{r}}_{x}({t}_{f})$$
(16d)
$${p}_{5}({t}_{f})=-\,{\lambda }_{2}[\frac{m}{7}{\omega }_{y}^{2}({t}_{f}){r}_{y}(0){r}_{y}({t}_{f})-\frac{15gN{r}_{y}(0)}{28\pi \,{r}_{x}({t}_{f}){r}_{y}^{2}({t}_{f}){r}_{z}({t}_{f})}]$$
(16e)
$${p}_{6}({t}_{f})=-\,{\lambda }_{2}\frac{m}{7}{r}_{y}(0){\dot{r}}_{y}({t}_{f})$$
(16f)
$${p}_{7}({t}_{f})=-\,{\lambda }_{2}[\frac{m}{7}{\omega }_{z}^{2}({t}_{f}){r}_{z}(0){r}_{z}({t}_{f})-\frac{15gN{r}_{z}(0)}{28\pi \,{r}_{x}({t}_{f}){r}_{y}({t}_{f}){r}_{z}^{2}({t}_{f})}]$$
(16g)
$${p}_{8}({t}_{f})=-\,{\lambda }_{2}\frac{m}{7}{r}_{z}(0){\dot{r}}_{z}({t}_{f}).$$
(16h)

We use a standard first-order gradient algorithm which is adapted to the control problem under study. The optimization procedure proceeds as follows:

  1. (i)

    First we fix an initial control ramp u(t) arbitrarily, such as the linear ramp u(t) = t/tf for instance, or the STA ramp obtained from ref.15.

  2. (ii)

    We then compute the magnetic field Bbias(t) using Eq. (9) and we deduce the trap dynamics by calculating the trap motion z0(t) and the trap frequencies ωx(t), ωy(t) and ωz(t);

  3. (iii)

    Using the Verlet method50, we then solve Eqs (1) and (3) to simulate the condensate dynamics in the Thomas-Fermi regime from the initial time t = 0 to the final time t = tf  ;

  4. (iv)

    We calculate the adjoint state p(tf) at the end of the transport using Eq. (16) and we propagate p(t) backward in time until t = 0 using Eq. (14);

  5. (v)

    Finally, we add a first order correction to the control ramp by replacing the control function u(t) by [u(t) + δu(t)], where δu(t) = ε(∂Hp/∂u), ε being a small positive constant.

This procedure is repeated until convergence is reached.

Convergence

Figure 1 shows a typical example of convergence of this algorithm. The condensate is assumed to be initially at rest in the ground state of the initial trap. The initial control ramp is the shortcut-to-adiabaticity solution (see ref.15 for details). The weight parameters are λ1 = 1, λ2 = 5.105 and λ3 = 0.001. We have chosen in this example a large value for λ2 in order to impose a fast convergence for the control of the final size of the condensate. In practice the correction parameter ε has to be chosen small enough to insure the convergence of the optimization algorithm. Since the correction to the control ramp is introduced at first order only, decreasing the value of ε beyond a reasonable limit does not improve the accuracy of the optimization procedure but it slows down the convergence. In the present example we have chosen ε = 10−11. In Fig. 1, panel (a) shows the classical energy Ecl(tf) of the condensate at the end of the transport (tf = 150 ms in this case) as a function of the optimal control theory iteration number (logarithmic scaling). Panel (b) shows the quantum energy Equ(tf) of the condensate at the end of the transport as a function of the iteration number. Panel (c) shows the average classical energy \(\frac{1}{{t}_{f}}{\int }_{0}^{{t}_{f}}{E}_{cl}(t)dt\) of the condensate during the transport as a function of the iteration number.

Figure 1
figure 1

Example of convergence of the different cost functionals as a function of the optimal control theory iteration number: (a) Final classical energy in nK, (b) Final quantum energy in nK, (c) Average classical energy in nK. See text for details.

Full size image

Since the total cost functional given in Eq. (6) is characterized by a very large weight λ2 associated with the final quantum energy, we see that Equ(tf) is very quickly minimized, in about 1,000 iterations. This limit of 1,000 iterations is emphasized in Fig. 1 with a vertical dashed red line. Once this convergence is reached, the final 3D size of the condensate adopts the size of the ground state of the final trap and the size dynamics of the BEC is frozen. This convergence was obtained at the cost of a transient degradation of the final classical energy, which reaches a maximum of about 20 nK after about 60 iterations, but the final classical energy is then minimized very quickly to reach a near-zero value in about 1000 iterations. It is only when this first stage of convergence is reached (iteration number >1000) that the last cost functional, associated with a smaller weight λ3, starts to decrease. One can note that the convergence of the average classical energy during the transport [in panel (c)] is rather slow since it requires more than 107 iterations before it starts to stabilize at values close to 30 nK. This value can be compared with the energy of the condensate in the initial trap, which is close to 120 nK, and with the energy of the condensate in the final trap, close to 10 nK. The transient excitation during the transport is therefore relatively limited.

Comparison of different optimization procedures

In Fig. 2 the Shortcut-To-Adiabaticity (STA) transport ramp obtained in ref.15 (dotted blue line) is compared with two results obtained with the present optimal control technique (OCT). Note that strictly speaking, these two methods assume slightly different constraints and that, in principle, STA can be combined with the minimization of a cost functional. However, a quantitative study on this specific point is beyond the scope of this article where we concentrate mainly on developing the transport method with the OCT protocol. The correction parameter is ε = 10−10. The dashed green line labeled as “cl-OCT” shows the result obtained for the weight factors λ1 = 1, λ2 = 0 and λ3 = 5.510−4. The solid red line labeled as “qu-OCT” is for λ1 = 1, λ2 = 3.3 and λ3 = 5.510−4. The difference between these two OCT results lies in the fact that qu-OCT takes into account the influence of the finite size of the BEC in the cost functional, while cl-OCT considers the BEC as a classical point-wise particle. The BEC model used for cl-OCT is therefore similar to the model used in STA and these two approaches can be compared directly. The optimized time variation of the bias magnetic field Bbias(t) is shown as a function of time in the first panel (a). The duration of the transport is tf = 150 ms, and all results are plotted from t = 0 to t = 250 ms i.e. up to 100 ms after the end of the transport. This time interval was chosen in order to detect the eventual presence of a residual excitation at the end of the transport. The position [zA(t) − z0(t)] and velocity \([{v}_{A}(t)-{\dot{z}}_{0}(t)]\) offsets are shown in panels (b) and (c). Finally, Panels (d), (e) and (f) present the condensate size dynamics Δα(t) along the three coordinates α ≡ x, y or z, where \({\rm{\Delta }}\alpha (t)={r}_{\alpha }(t)/\sqrt{7}\) represents the width (standard deviation) of the Thomas-Fermi condensate wave function in the directions α ≡ x, y or z.

Figure 2
figure 2

Comparison of different optimization procedures. Shortcut-to-adiabaticity (STA): dotted blue line, classical optimal control (cl-OCT): dashed green line, quantum optimal control (qu-OCT): solid red line. (a) Bias magnetic field in Gauss as a function of time, (b) Position offset [zA(t) − z0(t)] in μm as a function of time, (c) Velocity offset \([{v}_{A}(t)-{\dot{z}}_{0}(t)]\) in μm/ms as a function of time, (d–f) Size dynamics of the condensate along the three coordinates x, y and z in μm as a function of time. The duration of the transport is tf = 150 ms. See text for details.

Full size image

We see in panels (b) and (c) that the three methods are very efficient for the control of the final average position and velocity of the BEC since the condensate is fully at rest in the center of the trap at the end of the transport and for all times t > tf = 150 ms. In addition, the transient position and velocity offsets during the transport reach similar values using these three different optimization methods. One can note in panels (b) and (c) that in terms of maximum transient offset in position and speed, from the two methods that we can compare directly, cl-OCT is a little better than STA (maximum offsets of 4.5 μm vs. 5.3 μm in position and 14 μm/ms vs. 22 μm/ms in speed) but this difference is not very significant in practice. The transient offsets of the qu-OCT approach are slightly larger than those of the cl-OCT method (with maximum offsets of 6.2 μm in position and 15 μm/ms in speed). Again this increase would be very benign in a practical implementation. Note finally that the three control fields Bbias(t) shown in panel (a) are relatively similar, with a fast initial decrease during the first half of the ramp, before 75 ms, followed by a much slower decrease afterward. A first conclusion of this study is therefore that, if one is mainly interested in the control of the average translational degree of freedom of the BEC, the STA approach, whose numerical implementation is much simpler than OCT, is sufficient.

It is in the size dynamics shown in panels (d), (e) and (f) that there is a striking difference between qu-OCT and the two other optimization methods. In terms of size dynamics, cl-OCT and STA give very similar results which consist in a persistent size excitation of the condensate after the transport. This result was already seen in ref.15 where it was shown that it was mainly the first quadrupole mode Q1 which was excited, thus explaining that the size oscillation along x, y and z is almost periodic after the transport. The qu-OCT approach is able to suppress efficiently this quadrupole-mode excitation and, at the end of the transport, the sizes Δx, Δy and Δz remain constant. We can therefore conclude that the introduction of a minimization goal for the quantum energy associated with the finite size dynamics of the condensate allows the qu-OCT transport ramp to prepare the true ground state of the final trap at t = tf  . When the size dynamics is not accounted for, as in the STA and cl-OCT approaches, it is impossible to insure the preparation of the lowest energy state in the final trap using short transport ramps.

The optimized OCT transport ramps were obtained using a Thomas-Fermi approximation in a 3D harmonic trap. We have therefore verified, by solving the 3D mean-field time-dependent Gross-Pitaevskii equation for the evolution of the time-dependent macroscopic condensate wave function ψ(x, y, z, t), that this control is robust when taking into account the anharmonicities and the rotation of the trap. The numerical method used for this calculation is described in ref.15. This result is illustrated in Fig. 3, showing the time evolution of the average atomic densities

$${P}_{x}(x,t)={\int }_{-{\rm{\infty }}}^{{\rm{\infty }}}dy{\int }_{-{\rm{\infty }}}^{{\rm{\infty }}}dz\,|\psi (x,y,z,t){|}^{2}$$
(17a)
$${P}_{y}(y,t)={\int }_{-{\rm{\infty }}}^{{\rm{\infty }}}dx{\int }_{-{\rm{\infty }}}^{{\rm{\infty }}}dz\,|\psi (x,y,z,t){|}^{2}$$
(17b)
$${P}_{z}(z,t)={\int }_{-{\rm{\infty }}}^{{\rm{\infty }}}dx{\int }_{-{\rm{\infty }}}^{{\rm{\infty }}}dy\,|\psi (x,y,z,t){|}^{2}$$
(17c)

along x, y and z. In each panel the black dashed line shows the expected center of mass trajectory obtained by solving Newton’s equations of motion (1). Similarly, the dotted blue lines highlight the expected widths obtained by solving the scaling equation (3). The condensate wave function follows clearly these predicted positions and widths. It is therefore clear from Fig. 3 that the controls predicted by the cl-OCT and qu-OCT methods are robust with respect to the anharmonicities and with respect to the inherent rotation of the trap in this realistic atom chip setup. In addition, the control of collective excitations using the qu-OCT approach appears clearly when comparing the lower line (qu-OCT ramp) of Fig. 3 with the upper line (cl-OCT ramp) of the same Figure.

Figure 3
figure 3

Condensate dynamics in the x, y and z directions using the cl-OCT ramp (upper line) and the qu-OCT ramp (lower line) shown in Fig. 2. The transport duration is tf = 150 ms. The average atomic density, solution of the time-dependent Gross-Pitaevskii equation, is shown as a function of time and position: (a) and (d) Px(x, t), (b) and (e) Py(y, t), (c) and (f) Pz(z, t). The black dashed lines show the expected center of mass trajectory. The dotted blue lines highlight the expected width dynamics according to the scaling approach. The dotted vertical white lines mark the time of the end of the transport. The total atom number is N = 105. See text for details.

Full size image

Influence of the transport duration

What remains to be seen is the efficiency of these various optimization procedures for different transport durations. Figure 4 shows in panel (a), for the three optimized ramps, the variation of the average translational energy

$$\langle {E}_{cl}\rangle =\frac{1}{{t}_{f}}{\int }_{0}^{{t}_{f}}{E}_{cl}(t)\,dt$$
(18)

as a function of the ramp duration tf  . Panel (b) shows, in the same conditions, the maximum position offset Max|zA(t) − z0(t)| during the transport. We see here that whatever the transport duration STA and cl-OCT are characterized by a very similar performance in terms of transient excitations. This confirms the advantage of the STA protocol in practical implementations, due to its overall simplicity when compared to cl-OCT. We also see that, on one hand, when the transport duration is larger than 140 ms (i.e. about 1.4 times the largest time scale associated with the trap), the transient excitations realized by the improved qu-OCT procedure are very close to the ones of cl-OCT and STA. On the other hand, for transport durations smaller than 140 ms larger transient excitations are obtained when using qu-OCT.

Figure 4
figure 4

Influence of the transport duration tf on: (a) the average translational energy 〈Ecl〉 of the condensate and (b) the maximum position offset |zA − z0| during the transport. Shortcut-to-adiabaticity (STA): dotted blue line, classical optimal control (cl-OCT): dashed green line, quantum optimal control (qu-OCT): solid red line. The weight parameters λ1, λ2 and λ3 are the same as those used in Fig. 2. See text for details.

Full size image

We could however verify that for all transport durations in the range 100 ms ≤ tf ≤ 200 ms, the qu-OCT method is able to minimize very efficiently the residual size excitations after the transport, a goal which is not achievable with the STA or cl-OCT procedures. This can be seen in Fig. 5, which shows the residual oscillation amplitudes

$${\rm{\Delta }}{\alpha }_{{\rm{res}}}=\frac{1}{2}[\mathop{{\rm{Max}}}\limits_{t > {t}_{f}}({\rm{\Delta }}\alpha )-\mathop{{\rm{Min}}}\limits_{t > {t}_{f}}({\rm{\Delta }}\alpha )]$$
(19)

of the size of the condensate Δα (standard deviation of the Thomas-Fermi condensate wave function) after the transport, for α ≡ x [panel (a)], α ≡ y [panel (b)], and α ≡ z [panel (c)]. The results shown in Figs 4 and 5 demonstrate that for tf ≥ 140 ms the residual size excitations of the condensate can be limited efficiently by optimal control and that this limitation does not introduce any detrimental effect on the transient excitation of the BEC. The same result can also be obtained by optimal control with shorter transport ramps, but at the cost of an increased transient excitation of the condensate.

Figure 5
figure 5

Residual oscillation amplitudes in the size dynamics after transport in the (a) x, (b) y and (c) z directions, as a function of the transport duration tf  . Shortcut-to-adiabaticity (STA): dotted blue line, classical optimal control (cl-OCT): dashed green line, quantum optimal control (qu-OCT): solid red line. The weight parameters λ1, λ2 and λ3 are the same as those used in Fig. 2. See text for details.

Full size image

Conclusion

In conclusion, we engineered optimal control theory protocols allowing for the fast, excitation-less transport of BECs over large distances compatible with a precision atom interferometric use. The ramps presented in this work relied on a single-parameter (bias magnetic field) optimization to shift the trap minimum position of the atom chip, promising a straightforward experimental implementation. The results of the OCT procedure relied on a scaling approach assuming a harmonic trapping. Real-life implementations on atom chips come with anharmonic corrections, mainly cubic in the direction of the transport, that scale with the position offset between the atoms and the trap minimum during the transport and with an inherent rotation of the trap. We demonstrated in this study, by solving 3D Gross-Pitaevskii equations for typical anharmonic and rotating chip traps, that the proposed OCT protocol does not compromise the target state solution even for very competitive ramp times of 150 ms. This also suggests a successful transfer to experiments. Moreover, we indicated by studying the impact of different transport durations, the methodology to follow in order to device the shortest ramps possible. Indeed, by quantifying the maximum offset induced by each ramp duration, every experimental implementation would be characterized by an anharmonicity range explored according to the specific trap configuration considered. This range determines, ultimately, the success of the ramp in reaching the ground state of the final trap. The positive outcome of this study suggests a natural generalization to the dual-species transport case. For this latter, no analytic neither intuitive solutions do exist. The STA approach generally fails since the two species experience different potential frequencies due their mass difference. A comparable OCT approach to the one adopted in this study, based on a pair of coupled mean-field equations, would allow to find trap trajectories that bring a quantum mixture to a target position in its ground state. Such a source will allow for precision interferometric measurements such as equivalence principle tests.