Optimal control of Bose-Einstein condensates in three dimensions

Ultracold gases promise many applications in quantum metrology, simulation and computation. In this context, optimal control theory (OCT) provides a versatile framework for the efficient preparation of complex quantum states. However, due to the high computational cost, OCT of ultracold gases has so far mostly been applied to one-dimensional (1D) problems. Here, we realize computationally efficient OCT of the Gross-Pitaevskii equation (GPE) to manipulate Bose-Einstein condensates in all three spatial dimensions. We study various realistic experimental applications where 1D simulations can only be applied approximately or not at all. Moreover, we provide a stringent mathematical footing for our scheme and carefully study the creation of elementary excitations and their minimization using multiple control parameters. The results are directly applicable to recent experiments and might thus be of immediate use in the ongoing effort to employ the properties of the quantum world for technological applications.


I. INTRODUCTION
Over the last decade, the ever increasing experimental toolbox of atomic, optical and molecular physics has lead to an exciting improvement in the control and understanding of complex quantum systems [1]. Recently, this has resulted in an important shift of paradigm. While quantum systems were previously mostly studied to check the validity of theoretical models, interest has now increased in their manipulation for specific technological applications. Prototypical examples for this shift of paradigm are atomic interferometers for quantum enhanced metrology [2][3][4], atomic field probes [5] and microscopes [6,7], inertial sensors [8], atomic clocks [9], or applications in quantum computing [10,11] and quantum simulation [12].
In many cases, these applications rely on the controlled preparation of a well-defined quantum many-body state with particular properties. One of the key experimental challenges is thus the efficient transfer of a system to such a state. Optimal control theory (OCT) is a mathematical tool to devise control strategies for this transfer [13]. It is well studied in many physical systems, ranging from atoms and molecules to solid-state systems [14][15][16][17][18].
In this work, we apply it to the control of a dilute atomic Bose-Einstein condensate (BEC), a system which is well described by the three-dimensional (3D) Gross-Pitaevskii equation (GPE) [19,20]. Such BECs form a versatile experimental platform for the storage, manipulation and probing of interacting quantum fields with * Electronic address: mennemann@acin.tuwien.ac.at † Electronic address: tim.langen@colorado.edu high precision [1]. In a seminal work Hohenester et al. [21] demonstrated that OCT [22] provides a highly efficient way to realize the transfer of a BEC to a target state, vastly outperforming more simple schemes. In this context, it has also been shown that OCT is robust against fluctuations and decoherence, and can also specifically take into account experimental constraints [23]. This has recently lead to first experimental demonstrations [24,25].
OCT of BECs has so far been mostly used in onedimensional (1D) settings, as the computational cost scales exponentially with the number of dimensions [26]. However, many experimental situations can only approximately be described by a 1D model, potentially limiting the applicability of OCT to real-life situations. In the following we demonstrate the first OCT of a BEC in all three spatial dimensions. We go beyond situations where a 1D approximation is feasible, thus significantly expanding the range of applicability of OCT for BECs. Moreover, we perform an analysis of the collective excitations that are created as a result of the control. These excitations are directly connected to the non-linear nature of the GPE and can only be fully captured and minimized in a 3D treatment including multiple control parameters. They are thus highly relevant in realistic experimental situations.

II. THE CONTROL PROBLEM
We start with a brief review of OCT, as well as of the description of BECs in terms of the GPE.
Here, V λ ≡ V (r, λ(t)) is an external potential that is characterized by a single or several control parameters denoted by the vector λ. Assuming that the wave function is normalized to unity, the coupling constant g = N 4π 2 a s /m is defined by the mass m, the s-wave scattering length a s and the number N of atoms in the BEC. For example, for ultracold gases of 87 Rb atoms, the atomic mass is given by m = 1.44 × 10 −25 kg and the s-wave scattering length by a s = 5.24 nm. Measuring length in units of l 0 = 1×µm, mass in units of the atomic mass and time in units of t 0 = ml 2 0 / , equation (1) can be written as which is the starting point for our considerations below.

B. Optimal control problem
We seek to find an optimal time-evolution of the mcomponent control parameter λ : (0, T ) → R m , λ(0) = λ 0 , λ(T ) = λ T , which steers the system from the initial state ψ 0 at time zero to a desired state ψ d at final time T . Without loss of generality we assume that ψ 0 and ψ d are ground state solutions of the stationary GPE corresponding to the smooth external potentials V λ0 and V λ T at times t = 0 and t = T , respectively, with fixed parameters λ 0 , λ T . To find the time evolution we apply well-known techniques from optimal control theory [21]. As cost functional, we use where u, v = R 3 u(r) * v(r) dr denotes the standard scalar product of u, v ∈ L 2 (R 3 ; C). The definition (3) is the generalization of the functional used in Refs. [21,23,[26][27][28][29][30] to a multi-component control parameter λ. The first term in J measures the proximity of ψ to the desired state ψ d at the end of the steering process. The expression F(ψ) = 1−| ψ d , ψ | 2 is known as the infidelity and provides a measure for the difference of ψ and ψ d . In detail, it quantifies the L 2 -norm of ψ's component that is orthogonal to ψ d . The second term regularizes the control trajectory to account for the fact that parameters can never be changed infinitely fast in a real experiment.
Here, γ > 0 sets the penalty for fast variations of λ(t).
For our examples below we find that already a very small value γ = 1 × 10 −6 yields a satisfactory regularization. Our goal is to minimize J(λ, ψ) subject to the constraint that ψ solves the GPE (Eq. (2)) with the initial condition given by the respective ψ 0 . To this end, one introduces the Lagrange function where p(r, t) acts as a generalized Lagrange multiplier [22]. At a local minimum (λ, ψ, p) of J, all three variational derivatives D p L(λ, ψ, p)[δp], D ψ L(λ, ψ, p)[δψ] and D λ L(λ, ψ, p)[δλ] vanish for all admissible variations δp, δψ and δλ, respectively. The corresponding three conditions constitute the optimality system together with the initial and terminal conditions In general, no analytical solutions are available for (5) with (6). Here we use an iterative method to find a numerical approximation of the solution. For this purpose it is useful to introduce the reduced cost functional where ψ λ denotes the unique solution of the Gross-Pitaevskii equation for a given control parameter curve λ. The goal is to find a local (or, preferably, even global) minimizer λ * ofĴ. The most straight-forward iterative procedure that can be employed is the method of steepest descent, To determine an appropriate step size α k , we perform a line search in each iteration: Here the upper index denotes the iteration step. A comment is due on the use of the gradient ∇Ĵ(λ k ) in (8).
Recall that the gradient ofĴ at λ with respect to a specific inner product (·, ·) X on the space X of admissible variations δλ is the uniquely determined element ∇Ĵ ∈ X such that (∇Ĵ, δλ) X = D λĴ (λ)[δλ] for all admissible variations δλ ∈ X. The gradient thus depends sensitively on the choice of the inner product (·, ·) X on X. It has been pointed out already in Ref. [27] that any admissible variation δλ must have a finite value in the penalty term, i.e., its weak time derivative ∂ t δλ must be square-integrable on (0, T ), and must respect the boundary conditions in (6c), i.e., δλ(0) = δλ(T ) = 0. A natural choice for (·, ·) X is thus the H 1 A calculation, which we present in the appendix, shows that this choice of (·, ·) X yields wherein ψ and p are solutions of (5a) and (6a) or (5b) and (6b), respectively. By definition, ∇Ĵ vanishes at the boundaries t = 0 and t = T , and so the iteration (8) preserves the boundary conditions (6c). We emphasize that the seemingly canonical choice of (·, ·) X as the standard L 2 -scalar product would not allow to specify boundary data for ∇Ĵ, which would result in a severe loss of stability of the optimization algorithm.

C. Implementation
In the situations considered below we found that the method of steepest descent (see Eq. (8)) works reliably. However, using more advanced methods the number of iterations needed to ensure convergence of the algorithm can be reduced significantly. In fact, our solver is based on the non-linear conjugate gradient scheme of Hager and Zhang [31], which has also been employed in Ref. [27] for optimal control of the one-dimensional GPE. We stress that all inner products and norms related to the nonlinear conjugate gradient scheme need to be expressed in terms of the inner product given in Eq. (10).
The reduced cost functional (7) needs to be evaluated several times per iteration. Moreover, at the beginning of each iteration a gradient vector needs to be determined using Eq. (11). Solutions to the time-dependent GPE (5a) and the adjoint equation (5b) are obtained via the time-splitting spectral method [32]. Initial and desired final states for a given potential are found by imaginary time propagation.
In order to accelerate the solving of the optimal control problem we perform all computations on the graphics processing unit (GPU) of a powerful graphics card. To speed up the calculations and ensure convergence of the algorithm we start each optimization with a coarse spatial grid and a relatively big time step ∆t. The result for λ is used as an input for another round of optimization on a finer grid. This procedure is repeated until the algorithm converges to a final time-evolution for λ. A detailed description of our implementation is given in Appendix B.

III. EXAMPLES
In the following we demonstrate the results of our scheme by considering three applications of increasing complexity, which are directly connected to recent experiments.

A. Harmonic oscillator potential
In the first application we study a Bose-Einstein condensate in an elongated harmonic potential. Initially, the trap frequencies are chosen such that the condensate is aligned along the y direction. Using a suitable timeevolution of the trap frequencies, we aim to rotate the condensate by π/2, while keeping it in the ground state of the external potential.
An example of the transition is visualized in Fig. 1. It can be understood as a toy example of a broad class of experimental protocols in which the trapping geometry is changed, e.g. to mode match different traps [34], to (de)compress a trap [35] or to transfer condensates into dynamical potentials for atomtronics [36,37]. Conceptually similar pulsed manipulations are also performed to focus BECs in time-of-flight expansion [38].

Trapping potential
The harmonic potential in this example is given by wherein the frequencies ω x and ω y can be set independently via the control parameters λ 1 and λ 2 . More precisely, we transform the external potential from an initial configuration with ω x = ω i x and ω y = ω i y at time t = 0 to a final configuration with ω x = ω f x and ω y = ω f y at the final time t = T . To this end, we parametrize ω x and ω y as with λ 1 (0) = 0, λ 1 (T ) = 1, λ 2 (0) = 0, λ 2 (T ) = 1. t = 0 ms t = 3 ms t = 6 ms t = 9 ms FIG. 1: Two-parameter optimal control of an elongated harmonic potential. The timescale of the control is T = 9 ms. Initially aligned along the y-direction, the condensate is dynamically transformed to be aligned along the x-direction. The black isosurface corresponds to the external trapping potential that is controlled using OCT, the blue isosurfaces visualize the atomic density. Note that for clarity only the lower half of the potential is shown. Also, here and throughout this work any trivial potential offset has been removed for simplicity and easier visualization. Its only effect is an overall phase shift of the wave function which is of no relevance to the optimization procedure. Animations of the full dynamics are available online [33].
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.

Numerical simulations
In the following simulations the number of atoms is N = 5000, the final time is set to T = 9 ms and ω z = 5 kHz. The initial configuration of the trapping potential is given by ω i x = 5 kHz and ω i y = 0.75 kHz, the final configuration by ω f x = 0.75 kHz and ω f y = 5 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 λ 1 and λ 2 are varied linearly. The corresponding time-evolution of the trap frequencies ω x and ω y is depicted in Fig. 2a.
In order to investigate the overlap of ψ with ψ d beyond the end of the control we continue the time-evolution with λ(t) = λ(T ) for t > T . We proceed analogously in the other examples. As can be seen from Fig. 2b the infidelity decreases only slightly until t = T and shows a strong oscillation for t > 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 Figs. 2c-e.
Next, we consider the result of the optimal control algorithm.
Using λ 0 (t) = [0.25 sin(πt/T ) + t/T, −0.25 sin(πt/T ) + t/T ] for t ∈ [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 Fig. 2f, the time-evolution of the corresponding infidelity in Fig. 2g. 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 > 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 Figs. 2h-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].

B. Loading of a toroidal trap
In the second application we consider the loading of a toroidal trap as shown in Fig. 3. Such toroidal traps have recently been employed to realize atomic analogues of electrical circuits to study superflow and hysteresis [39][40][41][42][43].

Trapping potential
The trapping potential is given by a slightly elongated harmonic potential and a Gaussian function centered at the origin of our coordinate system [44] . In an experiment this Gaussian function could for example correspond to a red-detuned laser beam realizing a repulsive dipole potential.
As illustrated in Fig. 3 we consider the transformation of the potential from an initial harmonic configuration with ω x = ω i x and V 0 = 0 at time t = 0 to a toroidal configuration with ω x = ω f x and V 0 = V * 0 at the final time t = T . Hence, a suitable parameterization of ω x and V 0 is given by where λ 1 (0) = 0, λ 1 (T ) = 1, λ 2 (0) = 0, λ 2 (T ) = 1.
In Eq. (12b), χ plays the role of a saturation function. The use of the saturation function ensures that V 0 remains positive -and thus experimentally realizable -for any possible choice of λ 2 . This does not restrict the original control problem, as every experimentally realizable trajectory V 0 (t) ≥ 0 can be parametrized through a suitable λ 2 (t) in V * 0 χ(λ 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 Fig. 4. The interpolating points are chosen such that χ always remains positive. Moreover, χ(0) = 0 and χ(1) = 1.

Numerical simulations
The following simulations are carried out using V * 0 = h × 30 kHz, w 0 = 5 µm , T = 9 ms and N = 5000. The frequencies ω y = 2.5 Hz and ω z = 5 kHz are kept constant during the simulation. The initial configuration of the confinement potential is characterized by ω i x = 1 kHz and V 0 (t = 0)/h = 0 kHz, whereas the final configuration is given by ω f x = 2.5 kHz and V 0 (T )/h = V * 0 . As in the previous example we consider first the case where the parameters ω x and V 0 are changed linearly (see Fig. 5a). Fig. 5b reveals that the associated infidelity does not drop at all until t = T . For t > 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 Figs. 5c-e, the final wave function still differs strongly from the wave function of the desired state after t * = 22 ms.
Let us now discuss the result of the optimal control algorithm. An optimal time-evolution of the control parameters is given in Fig. 5f. Intuitively this control can be understood as the result of two separate time-scales. During the first halve of the control, the trap frequency ω x is increased, while the limits imposed on λ 2 prohibit any change of V 0 . During the second halve, on the other hand, V 0 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 Fig. 5g. Furthermore, the infidelity remains bounded by 3 × 10 −3 for t > 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 Figs. 5h-j.

C. Splitting
In terms of technological applications, a particular noteworthy realization of BECs is achieved using atom chips [47,48]. On these chips micro-fabricated wires allow the precise manipulation of BECs using static, radio and microwave fields. As a third application we thus consider the splitting of a single condensate into two identical halves using such an atom chip [46]. A visualization is presented in Fig. 6. This splitting protocol has recently been used to study the non-equilibrium dynamics of 1D Bose gases, revealing subtle effects, such as prethermalization [49][50][51][52], generalized statistical ensembles [53] and the light-cone-like emergence of thermal correlations [54,55]. Moreover, it forms the basic building block for integrated matter-wave interferometers [56,57].

Trapping potential
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 The parameters are given by In the following simulations we consider 87 Rb atoms which are trapped in the 5S 1/2 F = 2, m F = 2 state where g F = 1/2. The trap parameters are The resulting dressed-state potential is given by [58] , ω RF the frequency of the RF radiation and B RF ⊥ denoting the component of the linear polarized dressing field B RF that is aligned perpendicular to the static field. As in [54] we use a detuning of ∆ RF (0) = −2π × 30 kHz from the m F = 2 → m F = 1 transition for the simulation. The Rabi-frequency is parameterized by the control parameter λ   The splitting of a Bose-Einstein condensate, as realized by a radial deformation of an initially harmonic potential into a double well [46]. The two gases in the final picture are completely decoupled, with no more overlap between the respective wave functions. Animations of the full dynamics are available online [33].
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 (cf. Fig. 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.

Numerical simulations: single-parameter control
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 Fig. 7a). This procedure is identical 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 Fig. 7b. Moreover, the infidelity shows the expected strong oscillations for t > T . A snapshot of the density at time t * = 22.5 ms is illustrated in Figs. 7c-e, revealing that there is large discrepancy between the computed state ψ and the desired state ψ d .
Next, we consider the result of the optimal control algorithm. We find that, irrespective of the specific choice of λ 0 , the algorithm always converges to approximately the same minimizer of the cost functional. The corresponding time-evolution of the Rabi-frequency is shown in Fig. 7f. We observe that the Rabi-frequency remains zero for the first few milliseconds. In fact, only about three milliseconds 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 timescale. 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 Fig. 7g). However, for t > T we again observe a strong oscillation. Snapshots of the density distribution at t * = 22.5 ms are given in Figs. 7h-j.

Bogoliubov-de Gennes analysis
Interestingly, the 6 ms period of the very regular infidelity oscillation shown in Fig. 7g for the optimized splitting is approximately the same as the period of the infidelity oscillation depicted in Fig. 7b 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 Bogoliubovde Gennes (BdG) framework.
In a linear approximation (with respect to δψ) this small deviation is given by where u, v and ω are defined via the solutions of the BdG equations [59,60] (15) We want to investigate small fluctuations δψ corresponding to some of the lowest energy eigenvalues ω 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 Bose-Einstein condensates. Like in Ref. [61] we consider the full three-dimensional 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 ω 1 = ±314.54 Hz, ω 2 = ±523.49 Hz and ω 3 = ±734. 26 Hz. Subsequently, the corresponding eigenfunctions (u i , v i ) are normalized according to the norm [60] R 3 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 Eq. (14). It turns out that |δψ i (t)| 2 can be well described by a simple periodic oscillation in amplitude, while the shape remains mostly unchanged (see left column in Fig. 8). As u i and v i are purely real-valued functions, which approximately fulfill v i = −u i (see right column in Fig. 8) we find δψ i (r, T i /2) ≈ δψ i (r, T i ) and hence the effective oscillation periods are halved with respect to the eigenvalues found above, i.e. T eff,i = T i /2 = π/ω i . In detail we find T eff,1 = 9.99 ms, T eff,2 = 6.00 ms and T eff,3 = 4.28 ms.
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 timeevolution of | ψ d , δψ 2 (t) | 2 along with the time-evolution of the infidelity in Fig. 7g 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 δψ 2 (t) (see 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 ≥ T . Left: Density of the first three (scaled) excitations δψ1(r, t), δψ2(r, t) and δψ3(r, t) at t = T eff,1 /2, t = T eff,2 /2 and t = T eff,3 /2. Right: Normalized (with respect to the inner product (16)) amplitude functions u and v evaluated along the longitudinal direction at x = xs and z = 0. All functions are purely real-valued.
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 Rabifrequency. The BEC will thus oscillate for t > 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 one-dimensional limit, where ω x,z ω 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 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.

Numerical simulations: two-parameter control
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].
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 ω y (t = 0) = 2π × 85 Hz and ω y (t = T ) = 2π × 85 Hz such that the initial and desired final states remain unchanged. Using λ 0 (t) = [t/T, 1] for t ∈ [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 Fig. 9a. As can be seen from Fig. 9b the infidelity remains very low for t ≥ T . Snapshots of the density distribution at time t = 22.5 ms confirm that the deviation from the desired state is extremely small, see Figs. 9d-f.
In the given example we have chosen T = T eff,2 . In contrast to that, for a time T < T 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.

IV. CONCLUSION AND OUTLOOK
In this work we have presented the first optimal control of the GPE in 3D. As we have shown, this situation is inherently more difficult than the optimal control of the 1D GPE because of the non-linear coupling of different coordinate directions. We have performed a detailed analysis of the resulting small excitations, which we were able to minimize by extending previous control schemes from a single to a multi-parameter control.
In contrast to 1D approximations our 3D approach allows the study of realistic trapping potentials, which will have direct impact on the quality of experiments and therefore provide an important step in the ongoing effort to use the properties of the quantum world for real life applications. Importantly, our scheme is not limited to the examples discussed in this work but rather very flexible, with many more applications conceivable.
A straight-forward extension of our numerical solver could include the treatment of excited states. This would allow the three-dimensional study of a recent experiment, where the BEC was transferred to the first excited state of the trapping potential via a 1D optimal control sequence [24]. Based on our observations we expect an even stronger excitation of BdG-Modes in such an experiment. In that context, another interesting application would be to replace the cigar-shaped confinement potentials used in the splitting and vibrational state inversion experiments by torus-shaped trapping potentials. Due to the different topology the issues related to the excitation of small perturbations are expected to be strongly reduced.
Another obvious extension of this work could be to consider different cost functionals. More precisely, it would be interesting to investigate whether it is possible to reduce the optimization times T by using other cost functionals which are not based on the infidelity but rather on a conserved quantity like the total energy. Finally, interesting further directions include the study of beyond mean-field effects using the multiconfigurational time-dependent Hartree framework for bosons [63] or the optimization of finite temperature states.   In the main text we have introduced the cost functional J(λ, ψ) in (3) and the reduced cost functional J(λ) = J(λ, ψ λ ) in (7). Recall that, for a given control λ : (0, T ) → R m satisfying the boundary conditions (6c), ψ λ is the solution to the initial value problem (5a) and (6a) for the GPE with the corresponding potential V λ . Below, we argue why the H 1 -gradient ofĴ is given by the component Λ : (0, T ) → R m of the solution (ψ, p, Λ) to the system consisting of (5a), (5b) and subject to the initial and terminal conditions (6a), (6b) and Before discussing the gradient, we first calculate the variational derivative ofĴ. As it is customary in the context of optimization problems, we express the validity of the GPE (5a) in the form of a constraint Z = 0, with the contraint functional By definition, ψ λ satisfies Z(λ, ψ λ ) = 0, hencê where L denotes the Lagrangian which was defined in Eq. (4) of the main text.
Eq. (20) holds for arbitrary smooth functions p : (0, T ) → L 2 (R 3 ; C). For fixed p, differentiation ofĴ in the direction δλ yields where δψ is the variation in ψ λ induced by the variation δλ of λ, i.e., it satisfies D λ Z(λ, ψ λ )[δλ] + D ψ Z(λ, ψ λ )[δψ] = 0 and δψ(0) = 0. For simplification of D λĴ , we choose p, which has been arbitrary up to this point, such that the last two terms in (21) cancel. Indeed, taking p as a solution to the terminal value problem (5b) and (6b), it follows that To arrive at this result, we have performed an integration by parts with respect to time, using the terminal condition (6b) and the fact that δψ(0) = 0 thanks to the initial condition (6a). In view of these cancellations, equation (21) simplifies to We are now in the position to calculate the H 1 -gradient ofĴ. Recall that the Sobolev space H 1 (0, T ; R m ) consists of all square integrable functions λ ∈ L 2 (0, T ; R m ) that possess a weak derivative ∂ t λ ∈ L 2 (0, T ; R m ). Functions λ ∈ H 1 (0, T ; R m ) are actually Hölder continuous, and therefore, they have well-defined boundary values at t = 0 and t = T . It is natural to consider the reduced cost functionalĴ as defined on H 1 * (0, T ; R m ), which is the affine subspace of functions λ ∈ H 1 (0, T ; R m ) that satisfy the boundary conditions (6c). Indeed, any admissible control λ : (0, T ) → R m must produce a finite value in the penalty term in J, which implies that ∂ t λ ∈ L 2 (0, T ; R m ). The tangent space to H 1 * (0, T ; R m ), i.e., the space of possible variations δλ, is the linear subspace H 1 with vanishing boundary values, Λ(0) = Λ(T ) = 0. This is a Hilbert space with respect to the inner product By definition, the gradient ofĴ with respect to the inner product (·, ·) is the uniquely determined element Λ ∈ H 1 0 (0, T ; R m ) such that (Λ, δλ) = D λĴ (λ)[δλ] for all variations δλ ∈ H 1 0 (0, T ; R m ). In view of (22), Λ satisfies and Λ ∈ H 1 0 (0, T ; R m ) induces the boundary conditions (19). To verify that the solution Λ to the boundary value problem (18) and (19) satisfies (23), it sufficies to integrate by parts in the time integral on the left-hand side, using that δλ(0) = δλ(T ) = 0.
B. Algorithms and implementation

Numerical evaluation of the cost functional
The evaluatation of the reduced cost functional (7) for a given control curve λ implicitly involves the computation of ψ λ , 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 ψ λ in the following.
For the numerical computation of the first term in (3), that is 1/2 1 − | ψ d , ψ(T ) | 2 , we have to solve the GPE (5a) with initial data (6a) for t ∈ [0, T ]. Our simulations are performed on the spatial domain with L x , L y and L z chosen sufficiently large to capture the significant part of the rapidly decaying solution ψ.
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 (x jx , y jy , z jz ) = where x = L x /J x with j x = 0, ..., J x − 1 etc. Thus, at time t n , the wave function ψ(t n ) is represented by a threedimensional array of complex numbers ψ (n) ∈ C Jx,Jy,Jz .
The method is of second order in time and of spectral accuracy in space, provided that ψ 0 and V λ 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 O(J log J). However, in some simulations (Splitting), the costs for computing the external potential V λ (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 ψ and 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 1 − | ψ d , ψ(T ) | 2 with ψ(T ) ≈ ψ (N ) , using a quadrature formula. The integral γ 2 T 0 |∂ t λ(t)| 2 dt is computed by a quadrature formula as well, using a finite difference formula of second order for the approximation of the time derivative ∂ t λ.

Numerical computation of the gradient
According to (11a), the H 1 -gradientĴ(λ) is obtained as solution to the second order problem ψ (0) M −1 b. We found that the algorithm converges reasonably fast when the factors M 1 and M 2 are given by the matrices L and U obtained from a sparse incomplete LU -factorization. Such an approximate factorization of A − σI can be computed using another Matlab function called ilu. For further information about the Matlab functions mentioned above we refer to the Matlab documentation and the literature cited therein. The time needed to compute a few eigenvalue-eigenvector solutions of the Bogoliubov-de Gennes equations depends strongly on the number of grid points J x , J y and J z . For the number of grid points reported above the whole computation takes on the order of five hours computing time utilizing the above mentioned CPU. The small perturbation which causes the oscillation of the infidelity in the splitting example can be extracted directly from the time-evolution of the wave function ψ. To this end, we assume that ψ(r, t) and ψ d (r) are almost identical for t = T , i.e., ψ(t = T ) ≈ e iθ ψ d , θ = arg min θ ψ(T ) − e iθ ψ d .
This assumption is in good agreement with our observations, where the minimum value of the infidelity is reached at this point in time. In analogy to equation (13) we define the difference ψ := ψ(r, t) − Φ(r, t), which leads to the result ψ(r, t) := ψ(r, t) − e iθ ψ d (r)e −iµ(t−T )/ , t ≥ T.
Here, we have introduced an additional phase factor e iµT / in order to take into account that we consider the time-evolution of ψ starting at t = T . A snapshot of the density | ψ(r, t)| 2 for t = 9 ms is shown in Fig. 11. It is quite obvious that the distribution of the density is very similar to the distribution of the density of the second excitation depicted in Fig. 8.

E. One-dimensional approximation for the splitting of a BEC
We briefly discuss the 1D approximation for the splitting example. In this case, the reduced GPE for the x-direction is given by where the effective 1D interaction strength g 1d is found by integrating out the two transversal dimensions [68,69] g 1d ≈ g Here,φ(y, z) := φ(0, y, z) corresponds to the normalized ground-state solution of the 3D model in the (x ≡ 0)plane. With this approximation we find g 1d ≈ h × 1300.44 Hz µm for N = 2000 atoms. This value describes the situation along the whole x-axis and also leads to reasonable results away from the center of the cloud, as can be seen from Figs. 12a-b. We then follow the same procedures as in the 3D case to find an optimal control trajectory for the Rabi frequency. Here, φ3D,1(x) = φ3D(x, 7.5 µm, 1 µm), φ3D,2(x) = φ3D(x, 7.5 µm, 0 µm), φ3D,3(x) = φ3D(x, 0 µm, 1 µm) and φ3D,4(x) = φ3D(x, 0 µm, 0 µm). Each wave function has been normalized to unity. c) Optimal control of the Rabi-frequency corresponding to the 1D approximation. d) Infidelity (1D-1D) corresponding to the one-dimensional model and infidelity (1D-3D) when the same trajectory of the Rabi-frequency is applied in a simulation using the 3D model.