Optimal control of Bose–Einstein condensates in three dimensions

The harmonic potential in this example is given by

$$\begin{eqnarray*}&&{V}_{{\boldsymbol{\lambda }}}(x,y,z)=\displaystyle \frac{m}{2}\left({\left[{\omega }_{x}({\lambda }_{1})\right]}^{2}{x}^{2}+{\left[{\omega }_{y}({\lambda }_{2})\right]}^{2}{y}^{2}+{\omega }_{z}^{2}{z}^{2}\right),\end{eqnarray*}$$

wherein the frequencies ω_x and ω_y can be set independently via the control parameters λ₁ and λ₂. More precisely, we transform the external potential from an initial configuration with ω_x = ω_xⁱ and ω_y = ω_yⁱ at time t = 0 to a final configuration with ω_x = ω_x^f and ω_y = ω_y^f at the final time t = T. To this end, we parametrize ${\omega }_{x}$ and ω_y as

$$\begin{eqnarray*}{\omega }_{x}({\lambda }_{1}) & = & {\omega }_{x}^{i}+{\lambda }_{1}({\omega }_{x}^{f}-{\omega }_{x}^{i}),\\ {\omega }_{y}({\lambda }_{2}) & = & {\omega }_{y}^{i}+{\lambda }_{2}({\omega }_{y}^{f}-{\omega }_{y}^{i}),\end{eqnarray*}$$

with

$$\begin{eqnarray*}{\lambda }_{1}(0) & = & 0,\quad {\lambda }_{1}(T)=1,\\ {\lambda }_{2}(0) & = & 0,\quad {\lambda }_{2}(T)=1.\end{eqnarray*}$$

We note that these parametrizations, as all others discussed below, are chosen as an example and can easily be adjusted to the parameters accessible in a specific experimental realization.

In the following simulations the number of atoms is N = 5000, the final time is set to T = 9 ms and ${\omega }_{z}=5\;{\rm{kHz}}.$ The initial configuration of the trapping potential is given by ${\omega }_{x}^{i}=5\;{\rm{kHz}}$ and ${\omega }_{y}^{i}=0.75\;{\rm{kHz}},$ the final configuration by ${\omega }_{x}^{f}=0.75\;{\rm{kHz}}$ and ${\omega }_{y}^{f}=5\;{\rm{kHz}}.$

Before we discuss the result of the optimal control algorithm we first consider a numerical simulation as a benchmark, in which the control parameters λ₁ and λ₂ are varied linearly. The corresponding time-evolution of the trap frequencies ω_x and ω_y is depicted in figure 2(a).

Zoom In Zoom Out Reset image size

In order to investigate the overlap of ψ with ψ_d beyond the end of the control we continue the time-evolution with ${\boldsymbol{\lambda }}(t)={\boldsymbol{\lambda }}(T)$ for $t\gt T.$ We proceed analogously in the other examples. As can be seen from figure 2(b) the infidelity decreases only slightly until t = T and shows a strong oscillation for $t\gt T.$ This behavior of the infidelity indicates that the final state differs significantly from the desired state ψ_d. This is also strikingly visualized by example snapshots of the density at time ${t}^{*}=22$ ms in figures 2(c)–(e).

Next, we consider the result of the optimal control algorithm. Using ${{\boldsymbol{\lambda }}}^{0}(t)\;=\;$ $[0.25\mathrm{sin}(\pi t/T)+t/T,-0.25\mathrm{sin}(\pi t/T)+t/T]$ for $t\in [0,T]$ as a starting point, the algorithm converges to a solution that reduces the cost functional by four orders of magnitude. The time-evolution of the frequencies ω_x and ω_y is shown in figure 2(f), the time-evolution of the corresponding infidelity in figure 2(g). It can clearly be seen that the infidelity strongly decreases until the end of the control at t = T. Moreover, the infidelity remains on a very low level for $t\gt T,$ indicating that the desired final state has been reached with high precision. Consequently, the deviations of the density to the density of the desired state at time $t={t}^{*}$ are very close to zero as can be seen from figures 2(h)–(j). We note at this point that the evolution of the 3D wave functions can naturally only be described here in limited detail. A supplementary video that visualizes these dynamics in greater detail is available online [33].

The trapping potential is given by a slightly elongated harmonic potential and a Gaussian function centered at the origin of our coordinate system [40]

$$\begin{eqnarray*}{V}_{{\boldsymbol{\lambda }}}(x,y,z) & = & \displaystyle \frac{m}{2}\left({\left[{\omega }_{x}({\lambda }_{1})\right]}^{2}{x}^{2}+{\omega }_{y}^{2}{y}^{2}+{\omega }_{z}^{2}{z}^{2}\right)\\ & & +{V}_{0}({\lambda }_{2})\mathrm{exp}(-2({x}^{2}+{y}^{2})/{w}_{0}^{2}).\end{eqnarray*}$$

In an experiment this Gaussian function could for example correspond to a red-detuned laser beam realizing a repulsive dipole potential.

As illustrated in figure 3 we consider the transformation of the potential from an initial harmonic configuration with ω_x = ω_xⁱ and V₀ = 0 at time t = 0 to a toroidal configuration with ω_x = ω_x^f and ${V}_{0}={V}_{0}^{*}$ at the final time t = T. Hence, a suitable parameterization of ω_x and V₀ is given by

$$\begin{eqnarray}&&{\omega }_{x}({\lambda }_{1})={\omega }_{x}^{i}+{\lambda }_{1}({\omega }_{x}^{f}-{\omega }_{x}^{i}),\end{eqnarray} \tag{ 12a }$$

$$\begin{eqnarray}&&{V}_{0}({\lambda }_{2})={V}_{0}^{*}\;\chi ({\lambda }_{2}),\end{eqnarray} \tag{ 12b }$$

where

$$\begin{eqnarray*}{\lambda }_{1}(0) & = & 0,\quad {\lambda }_{1}(T)=1,\\ {\lambda }_{2}(0) & = & 0,\quad {\lambda }_{2}(T)=1.\end{eqnarray*}$$

In equation (12b), χ plays the role of a saturation function. The use of the saturation function ensures that V₀ remains positive—and thus experimentally realizable—for any possible choice of λ₂. This does not restrict the original control problem, as every experimentally realizable trajectory ${V}_{0}(t)\geqslant 0$ can be parametrized through a suitable ${\lambda }_{2}(t)$ in ${V}_{0}^{*}\chi ({\lambda }_{2}(t)).$ In fact, we choose all parameters for the external potential to be close to previous experimental realizations. However, our approach also allows us to optimize more general situations where the parametrization of the trapping potential is more complicated [45].

Similar saturation functions are commonly used in control theory to realize limits on control parameters. In our particular case χ is implemented using a piecewise cubic hermite interpolating polynomial (PCHIP). Its functional form is shown in figure 4. The interpolating points are chosen such that χ always remains positive. Moreover, $\chi (0)=0$ and $\chi (1)=1.$

Zoom In Zoom Out Reset image size

The following simulations are carried out using ${V}_{0}^{*}=h\times 30\;{\rm{kHz}},$ ${w}_{0}=5\;\mu {\rm{m}},\;$T = 9 ms and N = 5000. The frequencies ${\omega }_{y}=2.5\;{\rm{kHz}}$ and ${\omega }_{z}=5\;{\rm{kHz}}$ are kept constant during the simulation. The initial configuration of the confinement potential is characterized by ${\omega }_{x}^{i}=1\;{\rm{kHz}}$ and ${V}_{0}(t=0)/h=0\;{\rm{kHz}},$ whereas the final configuration is given by ${\omega }_{x}^{f}=2.5\;{\rm{kHz}}$ and ${V}_{0}(T)/h={V}_{0}^{*}.$

As in the previous example we consider first the case where the parameters ω_x and V₀ are changed linearly (see figure 5(a)). Figure 5(b) reveals that the associated infidelity does not drop at all until t = T. For $t\gt T$ we observe a slight decrease of the infidelity. This can be attributed to the fact that, as time evolves, the density of the condensate becomes more evenly distributed in the toroidal trapping potential, bringing its wavefunction closer to ψ_d. However, as can be seen from figures 5(c)–(e), the final wave function still differs strongly from the wave function of the desired state after ${t}^{*}=22\;$ ms.

Zoom In Zoom Out Reset image size

Let us now discuss the result of the optimal control algorithm. An optimal time-evolution of the control parameters is given in figure 5(f). Intuitively this control can be understood as the result of two separate time-scales. During the first half of the control, the trap frequency ω_x is increased, while the limits imposed on λ₂ prohibit any change of V₀. During the second half, on the other hand, V₀ is adjusted to its final value, while ω_x is only subject to small corrections.

Until the end of the control this leads to a drop in the infidelity by approximately three orders of magnitude, as visualized in figure 5(g). Furthermore, the infidelity remains bounded by $3\times {10}^{-3}$ for $t\gt T,$ which is well below the measurement sensitivity in typical experiments. Consequently, only slight deviations from the desired wavefunction at time ${t}^{*}=22$ ms can be observed in figures 5(h)–(j).

In the experiments the splitting is realized by dressing the static magnetic trapping potential with a strong near-field radio-frequency (RF) field. The unscaled static potential is given by ${V}_{\mathrm{static}}={g}_{{\rm{F}}}{\mu }_{B}{m}_{{\rm{F}}}| {\bf{B}}| ,$ with the magnetic field ${\bf{B}}=({B}_{x},{B}_{y},{B}_{z})$ being well approximated by the famous Ioffe–Pritchard form

$$\begin{eqnarray*}{B}_{x} & = & {B}_{1}x-\displaystyle \frac{{B}_{2}}{2}{xy},\\ {B}_{z} & = & -{B}_{1}z-\displaystyle \frac{{B}_{2}}{2}{zy},\\ {B}_{y} & = & {B}_{0}+\displaystyle \frac{{B}_{2}}{2}\left[{y}^{2}-\displaystyle \frac{1}{2}\left({x}^{2}+{z}^{2}\right)\right].\end{eqnarray*}$$

The parameters are given by ${B}_{0}={\hslash }{\omega }_{0}/{m}_{{\rm{F}}}{g}_{{\rm{F}}}{\mu }_{B},{B}_{1}=\sqrt{m{\omega }_{\perp }^{2}{B}_{0}/{m}_{{\rm{F}}}{g}_{{\rm{F}}}{\mu }_{B}}$ and ${B}_{2}=m{\omega }_{\parallel }^{2}/{m}_{{\rm{F}}}{g}_{{\rm{F}}}{\mu }_{B}.$ In the following simulations we consider ${}^{87}$ Rb atoms which are trapped in the $5{{\rm{S}}}_{1/2}\;F=2,{m}_{{\rm{F}}}=2$ state where g_F = 1/2. The trap parameters are

$$\begin{eqnarray*}{\omega }_{0} & = & 2\pi \times 390\;\mathrm{kHz},\\ {\omega }_{x}={\omega }_{z}\equiv {\omega }_{\perp } & = & 2\pi \times 2\;\mathrm{kHz},\\ {\omega }_{y}\equiv {\omega }_{\parallel } & = & 2\pi \times 85\;\mathrm{Hz}.\end{eqnarray*}$$

The resulting dressed-state potential is given by [58]

$$\begin{eqnarray*}{V}_{\lambda } & = & {g}_{{\rm{F}}}{\mu }_{B}{\tilde{m}}_{{\rm{F}}}\sqrt{{\left(\displaystyle \frac{{\hslash }{\omega }_{{RF}}}{| {g}_{{\rm{F}}}| {\mu }_{B}}-| {\bf{B}}| \right)}^{2}+{\left(\displaystyle \frac{{B}_{\mathrm{RF}\perp }(t)}{2}\right)}^{2}},\\ & = & {g}_{{\rm{F}}}{\mu }_{B}{\tilde{m}}_{{\rm{F}}}\sqrt{{{\rm{\Delta }}}_{\mathrm{RF}}{({\bf{r}})}^{2}+{{\rm{\Omega }}}_{\mathrm{rabi}}^{2}(t)}\end{eqnarray*}$$

with ${\tilde{m}}_{F}=2,{\omega }_{\mathrm{RF}}$ the frequency of the RF radiation and B_RF⊥ denoting the component of the linear polarized dressing field ${{\bf{B}}}_{\mathrm{RF}}$ that is aligned perpendicular to the static field. As in [54] we use a detuning of ${{\rm{\Delta }}}_{\mathrm{RF}}(0)=-2\pi \times 30\;\mathrm{kHz}$ from the ${m}_{{\rm{F}}}=2\to {m}_{{\rm{F}}}=1$ transition for the simulation. The Rabi-frequency is parameterized by the control parameter λ

$$\begin{eqnarray*}&&{{\rm{\Omega }}}_{\mathrm{rabi}}(\lambda )={{\rm{\Omega }}}_{\mathrm{rabi}}^{*}\;\chi (\lambda ),\quad {{\rm{\Omega }}}_{\mathrm{rabi}}^{*}=2\pi \times 155\;\mathrm{kHz}\end{eqnarray*}$$

wherein

$$\begin{eqnarray*}&&{\lambda }_{1}(0)=0,\quad {\lambda }_{1}(T)=1.\end{eqnarray*}$$

The control parameter λ mimics the situation in experiments, where the double well potential is controlled by changing the RF field amplitude through an RF current in a wire. For λ = 0 we recover the static harmonic potential, whereas λ = 1 corresponds to a fully separated double well with no wave function overlap between the two halves of the system. Since the Rabi-frequency is strictly positive in experiments we employ the same saturation function χ as in the previous example (see figure 4).

As the trapping potential is significantly changed during the splitting the atoms are radially displaced from their equilibrium position in the harmonic trap. Consequently, strong dipole and breathing oscillations are usually observed in experiments. This poses a strong limitation to the use of such systems as interferometers [56]. The minimization of such excitations is therefore one of the main motivations for our optimization.

We illustrate the splitting procedure for N = 2000 atoms and T = 6 ms.

In a first step we again consider the case where the Rabi-frequency is increased linearly (see figure 7(a)). This procedure is similar to the one that is typically used in experiments [49, 53]. At the final time t = T the infidelity has only decreased slightly as can be seen from figure 7(b). Moreover, the infidelity shows the expected strong oscillations for $t\gt T.$ A snapshot of the density at time ${t}^{*}=22.5$ ms is illustrated in figures 7(c)–(e), revealing that there is a large discrepancy between the computed state ψ and the desired state ψ_d.

Zoom In Zoom Out Reset image size

Next, we consider the result of the optimal control algorithm. We find that, irrespective of the specific choice of the initial guess λ⁰, the algorithm always converges to approximately the same minimizer of the cost functional. The corresponding time-evolution of the Rabi-frequency is shown in figure 7(f). We observe that the Rabi-frequency remains zero for the first few miliseconds. In fact, only about three miliseconds of the optimization time T are used for the transformation of the external potential. This behavior persists even if we increase the optimization time T, with the Rabi-frequency vanishing for an even longer initial period of time. The precise timescale depends on the parameters of the trap, as the optimization algorithm tries to find a compromise between longitudinal and radial directions.

Interestingly, our 3D control qualitatively resembles the result of a previous 1D optimization that included beyond mean-field effects to model the distribution of atoms into the two final gases on the quantum level [57]. In both cases, the initial BEC is first rapidly split into two halves. Subsequently, these two halves are kept close enough to experience a tunnel coupling for a finite time-scale. This qualitative observation is very interesting, as reducing relative number fluctuations can help to significantly enhance the sensitivity of such interferometers. A detailed study of how useful our control can be in this context will be a natural extension of this work.

As a result of the optimal control algorithm the infidelity at the final time T is reduced by more than two orders of magnitude (see figure 7(g)). However, for $t\gt T$ we again observe a strong oscillation. Snapshots of the density distribution at ${t}^{*}=22.5$ ms are given in figures 7(h)–(j).

Interestingly, the 6 ms period of the very regular infidelity oscillation shown in figure 7(g) for the optimized splitting is approximately the same as the period of the infidelity oscillation depicted in figure 7(b) for the simple linear splitting. This suggests that the character of the oscillation is determined by the intrinsic properties of the BEC rather than by the splitting protocol.

Indeed, we demonstrate in the following that the oscillations are caused by collective excitations of the BEC, which are created during, but irrespective of the details of the splitting process. To this end, we show that they are the result of a small deviation δ ψ from the desired state ψ_d, which can be described within the Bogoliubov–de Gennes (BdG) framework.

Let therefore ${\rm{\Phi }}({\bf{r}},t)=\phi ({\bf{r}}){{\rm{e}}}^{-{\rm{i}}\mu t/{\hslash }}$ denote an eigenstate solution of the GPE. Here, μ is the corresponding chemical potential and ϕ is a solution of the stationary GPE, ${H}_{0}\phi +g{| \phi | }^{2}\phi =\mu \phi ,$ with ${H}_{0}=-{{\hslash }}^{2}/2m{\rm{\Delta }}+V.$ We consider a generic state ψ which deviates from the eigenstate solution by a small fluctuation δ ψ, i.e,

$$\begin{eqnarray}&&\psi ({\bf{r}},t)\approx {\rm{\Phi }}({\bf{r}},t)+\delta \psi ({\bf{r}},t).\end{eqnarray} \tag{ 13 }$$

In a linear approximation (with respect to δ ψ) this small deviation is given by

$$\begin{eqnarray}&&\delta \psi ({\bf{r}},t)=\left(u({\bf{r}}){{\rm{e}}}^{-{\rm{i}}\omega t}+{v}^{*}({\bf{r}}){{\rm{e}}}^{{\rm{i}}{\omega }^{*}t}\right){{\rm{e}}}^{-{\rm{i}}\mu t/{\hslash }},\end{eqnarray} \tag{ 14 }$$

where u, v and ω are defined via the solutions of the BdG equations [59, 60]

$$\begin{eqnarray}\left[\begin{array}{lc}{H}_{0}-\mu +2\;g{| \phi | }^{2} & g{\phi }^{2}\\ -{(g{\phi }^{2})}^{*} & -{H}_{0}+\mu -2\;g{| \phi | }^{2}\end{array}\right]\left[\begin{array}{l}u\\ v\end{array}\right]={\hslash }\omega \left[\begin{array}{l}u\\ v\end{array}\right].\end{eqnarray} \tag{ 15 }$$

We want to investigate small fluctuations δ ψ corresponding to some of the lowest energy eigenvalues ${\hslash }\omega$ in equation (15). To this end, we proceed in a conceptually similar way to [61] where numerical methods are used to investigate the stability and decay rates of non-isotropic attractive BECs. Like in [61] we consider the full 3D problem. However, for the discretization of the operators in (15) we employ a high-order finite difference discretization rather than working in a Fourier basis. By gradually increasing the spatial resolution of the finite difference discretization we are able to verify the convergence of the algorithm. A detailed description of our implementation is again given in the appendix.

As an example, we find that the first three eigenvalues converge towards ${\omega }_{1}=\pm 314.54\;{\rm{Hz}},$ ${\omega }_{2}=\pm 523.49\;{\rm{Hz}}$ and ${\omega }_{3}=\pm 734.26\;{\rm{Hz}}.$ Subsequently, the corresponding eigenfunctions $({u}_{i},{v}_{i})$ are normalized according to the norm [60]

$$\begin{eqnarray}&&{\displaystyle \int }_{{{\mathbb{R}}}^{3}}\left({u}_{i}^{2}({\bf{r}})-{v}_{i}^{2}({\bf{r}})\right){\rm{d}}{\bf{r}}=1.\end{eqnarray} \tag{ 16 }$$

Knowing the frequencies ω_i and amplitude functions u_i and v_i, it is possible to investigate the time-evolution of the excitations given by equation (14). It turns out that ${| \delta {\psi }_{i}(t)| }^{2}$ can be well described by a simple periodic oscillation in amplitude, while the shape remains mostly unchanged (see left column in figure 8). As u_i and v_i are purely real-valued functions, which approximately fulfill v_i = −u_i (see right column in figure 8) we find $\delta {\psi }_{i}({\boldsymbol{r}},{{\mathcal{T}}}_{i}/2)\approx \delta {\psi }_{i}({\boldsymbol{r}},{{\mathcal{T}}}_{i})$ and hence the effective oscillation periods are halved with respect to the eigenvalues found above, i.e. ${{\mathcal{T}}}_{\mathrm{eff},i}={{\mathcal{T}}}_{i}/2=\pi /{\omega }_{i}.$ In detail we find ${{\mathcal{T}}}_{\mathrm{eff},1}=9.99$ ms, ${{\mathcal{T}}}_{\mathrm{eff},2}=6.00$ ms and ${{\mathcal{T}}}_{\mathrm{eff},3}=4.28$ ms.

Zoom In Zoom Out Reset image size

Note that the effective period of the second excitation is very close to the period of the oscillation of the infidelity observed above. Indeed, plotting the time-evolution of ${| \langle {\psi }_{{\rm{d}}},\delta {\psi }_{2}(t)\rangle | }^{2}$ along with the time-evolution of the infidelity in figure 7(g) demonstrates clearly that the oscillation of the infidelity is dominated by the second excitation. As further evidence, we extract the deviation of ψ from ψ_d from our simulation. A comparison shows again very good agreement with the time-evolution of $\delta {\psi }_{2}(t)$ (see the appendix).

The fact that only the second but not the first excitation contributes to the observations can be understood from symmetry arguments. The first excitation corresponds to an antisymmetric wave function with respect to the longitudinal direction, whereas the second excitation is symmetric. During the splitting process, the halving of the atom number in each of the two gases, as well as an overall change in the longitudinal trapping potential leads to a symmetric change in the extension of the BEC in this direction. If the control is unable to compensate for this change in extension, the second Bogoliubov–de Gennes mode is automatically excited.

This effect is especially pronounced for the linear splitting. In contrast to that, the optimal control algorithm can still reduce the infidelity at t = T, but even a small deviation of the wave function from the stationary state leads to a strong oscillation in the infidelity for $t\geqslant T.$

Once the wave function differs from the stationary state in the longitudinal direction it is impossible to stop the observed oscillation by a simple variation of the Rabi-frequency. The BEC will thus oscillate for $t\gt T$ after the end of the control.

A central role in this scenario is played by the longitudinal frequency ω_y. The smaller ω_y the longer the extension of the condensate in the longitudinal direction. In analogy to a classical harmonic oscillator this increases the susceptibility to small deviations from the equilibrium position. We have confirmed this intuition with additional simulations, finding an even more pronounced excitation of the second mode for smaller ω_y.

This is particularly noteworthy with respect to experiments studying BECs in the 1D limit, where ${\omega }_{x,z}\gg {\omega }_{y}$ [53]. Intuitively, such experiments should be very well described through a 1D approximation, where only a reduced GPE for the x-direction has to be considered (see the appendix). Our results here show that such an approach will, in general, also lead to a strong breathing oscillation. Even if the 1D control is able to reach the 1D desired state with high precision, it does not necessarily describe the experimental reality and will thus fail in 3D.

In the last part of this article we will show how the oscillations reported above can be eliminated using a more sophisticated control scheme that is made possible by the 3D character of our control and that involves a manipulation of the trapping potential along the longitudinal direction. In experiments on atom chips, this manipulation can, for example, be realized using additional wire structures, which provide longitudinal confinement independent of the main radial trapping structures [62].

In analogy to the previous examples, we consider the following parameterization of ${{\rm{\Omega }}}_{\mathrm{rabi}}$ and ω_y:

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{rabi}}({\lambda }_{1})={{\rm{\Omega }}}_{\mathrm{rabi}}^{*}\;\chi ({\lambda }_{1}),\quad {{\rm{\Omega }}}_{\mathrm{rabi}}^{*}=2\pi \times 155\;\mathrm{kHz},\end{eqnarray} \tag{ 17a }$$

$$\begin{eqnarray}&&{\omega }_{y}({\lambda }_{2})={\omega }_{y}^{*}\;{\lambda }_{2},\quad {\omega }_{y}^{*}=2\pi \times 85\;\mathrm{Hz},\end{eqnarray} \tag{ 17b }$$

with

$$\begin{eqnarray*}{\lambda }_{1}(0) & = & 0,\quad {\lambda }_{1}(T)=1,\\ {\lambda }_{2}(0) & = & 1,\quad {\lambda }_{2}(T)=1.\end{eqnarray*}$$

The only difference to the previous example is thus that the value of the longitudinal trap frequency ω_y is now part of the control. We still fix ${\omega }_{y}(t=0)=2\pi \times 85\;{\rm{Hz}}$ and ${\omega }_{y}(t=T)=2\pi \times 85\;{\rm{Hz}}$ such that the initial and desired final states remain unchanged.

Using ${{\boldsymbol{\lambda }}}^{0}(t)=[t/T,1]$ for $t\in [0,T]$ as an initial guess the optimization algorithm converges to a solution which reduces the cost functional by more than three orders of magnitude. The time-evolution of the corresponding physical parameters is given in figure 9(a). As can be seen from figure 9(b) the infidelity remains very low for $t\geqslant T.$ Snapshots of the density distribution at time ${t}^{*}=22.5$ ms confirm that the deviation from the desired state is extremely small, see figures 9(d)–(f).

Zoom In Zoom Out Reset image size

In the given example we have chosen $T={{\mathcal{T}}}_{\mathrm{eff},2}.$ In contrast to that, for a time $T\lt {{\mathcal{T}}}_{\mathrm{eff},2}$ we find significantly worse results. The minimum time scale T is thus set by the oscillation period of the excitation that the control aims to stop. This oscillation period is in turn set by the geometry of the trap. Each different experimental situation will thus require carefully chosen parameters for the control.

The evaluatation of the reduced cost functional (7) for a given control curve ${\boldsymbol{\lambda }}$ implicitly involves the computation of ${\psi }_{{\boldsymbol{\lambda }}},$ that is, the solution of the GPE. No analytical solutions are available in general, so we use a numerical approximation. For brevity of notation, we write ψ instead of ${\psi }_{{\boldsymbol{\lambda }}}$ in the following.

For the numerical computation of the first term in  (3), that is $1/2\left(1-{| \langle {\psi }_{{\rm{d}}},\psi (T)\rangle | }^{2}\right),$ we have to solve the GPE (5a) with initial data (6a) for $t\in [0,T].$ Our simulations are performed on the spatial domain

$$\begin{eqnarray*}&&{\rm{\Omega }}=[-{L}_{x}/2,{L}_{x}/2]\times [-{L}_{y}/2,{L}_{y}/2]\times [-{L}_{z}/2,{L}_{z}/2]\end{eqnarray*}$$

with ${L}_{x},{L}_{y}$ and L_z chosen sufficiently large to capture the significant part of the rapidly decaying solution ψ.

For numerical discretization in time, we employ the following time-splitting spectral method [32]:

$$\begin{eqnarray}&&\psi ({t}_{n+1})\approx {{\rm{e}}}^{-{\rm{i}}{B}_{n}^{+}\triangle t/2}{{\rm{e}}}^{-{\rm{i}}A\triangle t}{{\rm{e}}}^{-{\rm{i}}{B}_{n}^{-}\triangle t/2}\psi ({t}_{n}),\end{eqnarray} \tag{ 24 }$$

with operators $A=-1/2{\rm{\Delta }},{B}^{\pm }={V}_{{{\boldsymbol{\lambda }}}^{(n+1/2)}}+g{| {\psi }_{n}^{\pm }| }^{2},$ and with ${t}_{n}=n\triangle t,n=0,\ldots ,N-1$ such that $N\triangle t=T.$ Here ${{\boldsymbol{\lambda }}}^{(n+1/2)}=1/2\;({\boldsymbol{\lambda }}({t}_{n})+{\boldsymbol{\lambda }}({t}_{n+1})),$ and the choice of ${\psi }_{n}^{\pm }$ is given below. Thus, the nth time step consists of the following three sub-steps. First, solve ${\rm{i}}{\partial }_{t}\psi =({V}_{{{\boldsymbol{\lambda }}}^{(n+1/2)}}+g{| \psi ({t}_{n})| }^{2})\psi$ for a duration of $\triangle t/2$ with initial value $\psi ({t}_{n});$ thus ${\psi }_{n}^{-}=\psi ({t}_{n}).$ The result is used as initial value for the free Schrödinger equation ${\rm{i}}{\partial }_{t}\psi =-1/2{\rm{\Delta }}\psi ,$ which is then solved for duration of $\triangle t;$ the result is ${\psi }_{n}^{+}.$ Finally, ${\rm{i}}{\partial }_{t}\psi =({V}_{{{\boldsymbol{\lambda }}}^{(n+1/2)}}+g{| {\psi }_{n}^{+}| }^{2})\psi$ is solved with initial value ${\psi }_{n}^{+},$ again for a duration of $\triangle t/2.$ The result of the third sub-step is taken as $\psi ({t}_{n+1}).$

The free Schrödinger equation is solved using the Fourier spectral method. To this end, the wave function ψ is interpolated by a trigonometric polynomial on the grid points of the cartesian grid

$$\begin{eqnarray*}&&({x}_{{j}_{x}},{y}_{{j}_{y}},{z}_{{j}_{z}})=(-{L}_{x}/2+{j}_{x}\triangle x,-{L}_{y}/2+{j}_{y}\triangle y,-{L}_{z}/2+{j}_{z}\triangle z),\end{eqnarray*}$$

where $\triangle x={L}_{x}/{J}_{x}$ with ${j}_{x}=0,\ldots ,{J}_{x}-1$ etc. Thus, at time t_n, the wave function $\psi ({t}_{n})$ is represented by a 3D array of complex numbers ${{\boldsymbol{\psi }}}^{(n)}\in {{\mathbb{C}}}^{{J}_{x},{J}_{y},{J}_{z}}.$

As Matlab code, the nth time step looks as follows:

$$\begin{eqnarray}{\boldsymbol{\psi }} & = & \mathrm{exp}(-1i\;\ast \;({{\bf{V}}}_{{{\boldsymbol{\lambda }}}^{(n+1/2)}}+g\;\ast \;\mathrm{abs}({\boldsymbol{\psi }}).\hat{}2)\;\ast \;\triangle t/2)\;.\;\ast \;\;{\boldsymbol{\psi }};\\ {\boldsymbol{\psi }} & = & \mathrm{fftn}({\boldsymbol{\psi }});\\ {\boldsymbol{\psi }} & = & {\bf{M}}\;.\;\ast \;\;{\boldsymbol{\psi }};\\ {\boldsymbol{\psi }} & = & \mathrm{ifftn}({\boldsymbol{\psi }});\\ {\boldsymbol{\psi }} & = & \mathrm{exp}(-1i\;\ast \;({{\bf{V}}}_{{{\boldsymbol{\lambda }}}^{(n+1/2)}}+g\;\ast \;\mathrm{abs}({\boldsymbol{\psi }}){.}^{}\hat{}2)\;\ast \;\triangle t/2)\;.\;\ast \;\;{\boldsymbol{\psi }};\end{eqnarray} \tag{ 25 }$$

The array ${\bf{M}}$ represents the action of the free Schrödinger operator. Due to the vectorized implementation in Matlab this procedure is highly efficient.

The method is of second order in time and of spectral accuracy in space, provided that ψ₀ and ${V}_{{\boldsymbol{\lambda }}}$ are sufficiently smooth. In comparison to a finite difference Crank–Nicolson scheme (see for example [21, 32]), the solution of a linear evolution equation is avoided, and less grid points J = J_x J_y J_z are needed to achieve the same quality of approximation for ψ.

Typically, the numerical costs for our implementation (25) are dominated by the fast Fourier transforms fftn and ifftn, which are of order ${\mathcal{O}}(J\mathrm{log}J).$ However, in some simulations (splitting), the costs for computing the external potential ${V}_{{{\boldsymbol{\lambda }}}^{(n+1/2)}}$ exceed that of the Fourier transforms.

For the numerical solution of the optimization problem, on the order of 10 to 100 evaluations of the cost functional are needed. The respective solution of the time-dependent GPE is performed on the graphics processing unit (GPU) of a powerful graphics card. Thanks to the vectorized implementation (25), it suffices to initialize the arrays ${\boldsymbol{\psi }}$ and ${\bf{V}}$ once at the beginning, using the Matlab command gpuArray. For handling the intermediate results or for calling the data in the memory of the main processor at the end of the computation we use the command gather. The trap potentials need to be updated in each time step. However, these calculations can be performed in a vectorized way on the GPU as well.

Finally, we compute $1/2\left(1-{| \langle {\psi }_{{\rm{d}}},\psi (T)\rangle | }^{2}\right)$ with $\psi (T)\approx {{\boldsymbol{\psi }}}^{(N)},$ using a quadrature formula. The integral $\displaystyle \frac{\gamma }{2}{\displaystyle \int }_{0}^{T}{| {\partial }_{t}{\boldsymbol{\lambda }}(t)| }^{2}\;{\rm{d}}t$ is computed by a quadrature formula as well, using a finite difference formula of second order for the approximation of the time derivative ${\partial }_{t}{\boldsymbol{\lambda }}.$

According to  (11a), the H¹-gradient $\hat{J}({\boldsymbol{\lambda }})$ is obtained as solution to the second order problem

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}}^{2}}{{\rm{d}}{t}^{2}}\left[{\rm{\nabla }}\hat{J}({\boldsymbol{\lambda }})\right]={\bf{r}}({\boldsymbol{\lambda }}),\quad {\bf{r}}({\boldsymbol{\lambda }})=\gamma \ddot{{\boldsymbol{\lambda }}}+\mathrm{Re}\langle \psi ,({\partial }_{{\boldsymbol{\lambda }}}{V}_{{\boldsymbol{\lambda }}})p\rangle \end{eqnarray} \tag{ 26 }$$

subject to the boundary conditions $\left[{\rm{\nabla }}\hat{J}({\boldsymbol{\lambda }})\right](0)={\bf{0}}$ and $\left[{\rm{\nabla }}\hat{J}({\boldsymbol{\lambda }})\right](T)={\bf{0}}.$ The time derivatives are discretized by second order finite differences.

To evaluate the right-hand side ${\bf{r}}({\boldsymbol{\lambda }}),$ the functions ψ and p need to be determined for $t\in [0,T].$ First, the state equation (5a) is solved as described above. Then, the adjoint equation (5b) is solved backwards in time, for the terminal condition (6b). For solution of the adjoint equation, a time-splitting method is applied as well: we alternately solve the equations ${\rm{i}}{\partial }_{t}p={V}_{{\boldsymbol{\lambda }}}p+2g{| \psi | }^{2}p+g{\psi }^{2}{p}^{*}$ and ${\rm{i}}{\partial }_{t}p=-1/2{\rm{\Delta }}p.$ The free Schrödinger equation is discretized by the Fourier-spectral method, and the value of ${\partial }_{{\boldsymbol{\lambda }}}{V}_{{\boldsymbol{\lambda }}}$ at time $t=(n-1/2)\triangle t$ is computed by means of the complex-step derivative approximation [64].

For integration of the ajoint equation on the time interval $[(n-1)\triangle t,n\triangle t],$ an approximation of the wave function ${\psi }^{(n-1/2)}=1/2\;({\psi }^{(n-1)}+{\psi }^{(n)})$ is needed. Since it is impossible to store the arrays ${{\boldsymbol{\psi }}}^{(n)}$ for every time step $n=0,\ldots ,N$ on the graphics card, the state equation is simultaneously solved backwards in time as well. The procedure is sketched in figure A1 : the calculation of ${{\bf{p}}}^{(n-1)}$ involves only two instances of the wave function and the 'old' adjoint state—that is, ${{\boldsymbol{\psi }}}^{(n)},{{\boldsymbol{\psi }}}^{(n-1)},$ and ${{\bf{p}}}^{(n)}.$ As soon as the approximations of ${{\boldsymbol{\psi }}}^{(n-1)}$ and ${{\bf{p}}}^{(n-1)}$ are available, also ${{\bf{r}}}^{(n-1)}$ can be computed. In this way it is enough to store at each time step four arrays in 3D, and the values of all available ${{\bf{r}}}^{(n)}$ with n = 0, ..., N (the storage space of which is neglegible).

Zoom In Zoom Out Reset image size

A further difficulty in the numerical computation of the adjoint equation arises from the conjugate-complex quantity p* in $g{\psi }^{2}{p}^{*}.$ Without going into details, we refer to the implementation in [27], which can be easily applied to the 3D case. As in the case of the GPE the computation of ${\bf{r}}$ can be significantly accelerated by using the graphics card. Still, the costs for the computation of the gradient are three to four times higher than for the evaluation of the cost functional.

The initial and terminal states ψ₀ and ψ_d are assumed to be ground state solutions of the stationary GPE. We compute them by imaginary time propagation [65, 66] (also known as normalized gradient flow): the time step $\triangle t$ in (25) is replaced by $-i\triangle t,$ and the wave function ϕ is normalized after every time step. By using adaptive time stepping, we reach a sufficiently exact solution with justifiable numerical costs.

For the numerical solution of the considered optimal control problems we use a personal computer $({\rm{i}}7-4770{\rm{K}}$ CPU $@3.50$ Ghz $\times 8$) and Matlab. The parts with the highest numerical costs, thus the solving of the partial differential equations and the computation of the external potentials, are performed on the graphics card (GeForce GTX TITAN), which accelerates the calculations significantly. The evaluation of the Fourier transform, for example, on the finest space discretization can be accelerated by a factor 4−6. In this context, it is important to mention that the CPU-version of fftn in Matlab is parallelized as well and hence uses all cores available on the CPU.

In general it is useful to initially solve each optimal control problem with a small number of Fourier modes J_x, J_y, J_z and with a relatively big time step $\triangle t.$ Subsequently, the same optimal control problem is solved on a finer mesh grid and with smaller time step, whereby as initial data ${{\boldsymbol{\lambda }}}^{0}$ is used, obtained as approximated solution in the computation before. We repeat this procedure until the computed control curve with respect to the old discretization does not differ from the control curve of the finer discretization anymore.

We consider a sequence of discretization parameters

$$\begin{eqnarray*}({J}_{x}^{(1)},{J}_{y}^{(1)},{J}_{z}^{(1)},{(\triangle t)}^{(1)}) & \to & ({J}_{x}^{(2)},{J}_{y}^{(2)},{J}_{z}^{(2)},{(\triangle t)}^{(2)})\\ & \to & \ldots \to ({J}_{x}^{(M)},{J}_{y}^{(M)},{J}_{z}^{(M)},{(\triangle t)}^{(M)})\end{eqnarray*}$$

with ${J}_{x}^{({\ell }+1)}\gt {J}_{x}^{({\ell })},{J}_{y}^{({\ell }+1)}\gt {J}_{y}^{({\ell })},{J}_{z}^{({\ell }+1)}\gt {J}_{z}^{({\ell })}$ and ${(\triangle t)}^{({\ell }+1)}\lt {(\triangle t)}^{({\ell })}$ for ${\ell }=2,\ldots ,M.$ Typically on the order of 10 to 100 iterations of the conjugate gradient method are needed for solving the optimal control problems on the coarse grid. The computational time is of some minutes. The numerical costs for the calculations of each single iteration increase rapidly with each discretization level. In the same time if one gets near to the local minimum, less iterations for finding the local minimum are required. By means of the described strategy each of the presented optimal control problems can be solved in several hours computing time with respect to the finest discretization level.