<figure><img src="image_1.png" /></figure>

HAL
open science

# Optimal control of a Bose-Einstein condensate in an optical lattice: the non-linear and two-dimensional cases

E. Dionis, B Peaudecerf, Stéphane Guérin, David Guéry-Odelin, D. Sugny

► To cite this version:

E. Dionis, B Peaudecerf, Stéphane Guérin, David Guéry-Odelin, D. Sugny. Optimal control of a Bose-Einstein condensate in an optical lattice: the non-linear and two-dimensional cases. Frontiers in Quantum Science and Technology, 2025, 4, pp.1540695. 10.3389/frqst.2025.1540695. hal-04929609

HAL Id: hal-04929609

https://cnrs.hal.science/hal-04929609v1

Submitted on 4 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Optimal control of a Bose-Einstein Condensate in an optical lattice: The non-linear and two-dimensional cases

E. Dionis, B. Peaudecerf,\* S. Guérin,† D. Guéry-Odelin,‡ D. Sugny§

January 31, 2025

## Abstract

We numerically study the optimal control of an atomic Bose-Einstein condensate in an optical lattice. We present two generalizations of the gradient-based algorithm, GRAPE, in the non-linear case and for a two-dimensional lattice. We show how to construct such algorithms from Pontryagin's maximum principle. A wide variety of target states can be achieved with high precision by varying only the laser phases setting the lattice position. We discuss the physical relevance of the different results and the future directions of this work.

## 1 Introduction

Quantum technologies seek to exploit the specific properties of quantum systems for real-world applications in computing, sensing, simulations or communications [1]. In this framework, quantum optimal control (QOC) can be viewed as a set of methods for designing and implementing external electromagnetic fields to realize specific operations on a quantum device in the best possible way [2–4]. QOC is becoming a key tool in many different experimental platforms, ranging from superconducting circuits [5–7] to cold atoms [8–10], molecular physics [11] or NV centers [12]. Despite the maturity and effectiveness of optimal control techniques which rely on a rigorous mathematical framework, namely the Pontryagin Maximum Principle (PMP, see [9, 13–15] for details), developments and adaptations of standard methods are necessary to take into account experimental limitations or additional degrees of freedom in the experiment [16–22]. Bose-Einstein condensates (BEC) constitute a promising system [23, 24] for applications in quantum sensing and simulation [25] in which optimal control may

<footnote>\*Laboratoire Collisions Agrégats Réactivité, UMR 5589, FeRMI, UT3, Université de Toulouse, CNRS, 118 Route de Narbonne, 31062 Toulouse CEDEX 09, France</footnote>

<footnote>†Laboratoire Interdisciplinaire Carnot de Bourgogne, CNRS UMR 6303, Université de Bourgogne, BP 47870, F-21078 Dijon, France</footnote>

<footnote>‡Laboratoire Collisions Agrégats Réactivité, UMR 5589, FeRMI, UT3, Université de Toulouse, CNRS, 118 Route de Narbonne, 31062 Toulouse CEDEX 09, France</footnote>

<footnote>§Laboratoire Interdisciplinaire Carnot de Bourgogne, CNRS UMR 6303, Université de Bourgogne, BP 47870, F-21078 Dijon, France, dominique.sugny@u-bourgogne.fr</footnote>
play a major role to prepare specific states or to improve the estimation of un-
known parameters [4]. This approach has been applied in a variety of works both
theoretically and experimentally, with very good agreement [8–10,26–47]. A spe-
cific example is given by a BEC trapped in an optical lattice. In this case, the
two control parameters are usually the depth and phase of the lattice that can
be precisely adjusted experimentally [26]. However, a majority of studies con-
sider a simplified situation in which the non-linear term of the Gross-Pitaevski
equation describing the dynamics of the BEC is neglected. This approxima-
tion that can be justified by the low density of the BEC is crucial and allows
to study the system in a unitary framework given by the Schrödinger equa-
tion. It also greatly simplifies the implementation of QOC as recently shown by
our group in [26, 46] where a gradient-based algorithm, GRAPE [17], has been
used. Different formulations of optimal control algorithms have been proposed
in the non-linear case. The optimal control of a BEC in a magnetic microtrap
has been investigated numerically in [8, 27, 28]. The control problem has been
solved analytically for a two-level quantum system in [48–52]. The control of
the Gross-Pitaevski equation has also been the subject of a series of mathe-
matical papers (see [53, 54] to mention a few). In this article, we propose to
revisit such works by applying the PMP in the presence of mean-field interac-
tions and deriving the corresponding algorithm from this optimization principle.
Intensive use of the pseudo-spectral approach [55–59] and FBR-DVR bases (Fi-
nite Basis Representation and Discrete Variable Representation) is necessary to
analytically express the different quantities and accelerate the numerical calcu-
lation. We demonstrate the effectiveness of this algorithm and discuss the role
of nonlinearity on the control procedure.

The use of one-dimensional lattices can limit the range of phenomena acces-
sible to quantum simulation: higher dimensionalities e.g. can radically change
the physics of localization [62], and also increase the role of interaction, giving
access to many-body phenomena such as the Mott transition [23]. It is thus of
utmost importance to extend the 1D optimal approach to the 2D or 3D cases.
This problem has been investigated by very few studies such as [33, 47]. In [33],
optimal control is applied to a BEC in a three dimensional magnetic trap. An
optimization algorithm different from GRAPE has been used to control BEC
in 2D and 3D lattices in [47]. We present in this work another numerical im-
plementation of QOC to a 2D lattice. We consider here a triangular geometry,
which is non-separable. This example can be used as a test-bed for other geome-
tries and for 3D lattices. State-to-state transfer can be optimized numerically
with a very good precision. We discuss on this example the controllability of
the target state with respect to the number of independent controls available.

The paper is organized as follows. Section 2 briefly recalls the specifics of
the experimental setup and the application of QOC in a 1D and linear cases.
Section 3 is dedicated to the extension of the GRAPE algorithm to the non-
linear case. Optimal control with two-dimensional lattices is described in Sec. 4.
Conclusion and prospective views are given in Sec. 5.

## 2 The standard one-dimensional case

Cold atom systems are characterized by their large size and the broad range
of controllable parameters they offer, which makes them excellent candidates
for applications in quantum technologies. We consider in this paper the BEC experiment in Toulouse in the group of D. Guéry-Odelin. Recent experimental results have shown the key role of QOC for state-to-state transfer in this setup [9, 26, 42, 46].

The experiment starts by laser cooling followed by evaporation of a Rubidium 87 gas allowing the formation of a BEC. The condensate is composed of $5 \times 10^5$ atoms at a temperature of 90 nK and is held in a hybrid magneto-optical trap to confine the BEC and compensate for gravity. A horizontal one-dimensional optical lattice trap, formed by the interference of two counter propagating laser beams with a wavelength $\lambda = 1064$ nm, is superimposed on the hybrid trap, aligned with the dipole trapping beam. Using acousto-optic modulators, the amplitude and phase of the lattice lasers can be shaped in time. This amounts to modifying the depth of the lattice or translating it, respectively. Following the experimental results of [26, 42, 46], we assume here that the lattice depth is fixed and only the relative phase of the lasers can be controlled.

The wave function $|\psi(t)\rangle$ describing the state of the BEC is governed by the time-dependent Schrödinger equation which can be expressed as

$$i\hbar \frac{d|\psi(t)\rangle}{dt} = \left( \frac{\hat{p}^2}{2m} - \frac{sE_L}{2} \cos(k_L \hat{x} + \varphi(t)) \right) |\psi(t)\rangle, \quad (1)$$

where $\hat{p}$ and $\hat{x}$ are respectively the momentum and the position operators, with $x$ the spatial coordinate along the axis of the optical lattice. We denote by $m$ the atom mass, $k_L = 2\pi/(\lambda/2)$ the wave vector and $E_L = (\hbar^2 k_L^2)/(2m)$ the characteristic energy of the lattice. The parameters $s$ and $\varphi$ correspond respectively to the relative amplitude of the lattice and to its phase. In standard experimental conditions, we stress that the atom interactions and the potential energy due to the hybrid trap, $V_{hyb} = \frac{1}{2}m\omega_{ext}x^2$ with $\omega_{ext} = 2\pi \times f$ where $f$ varies from a few Hz to 25 Hz at most, can be neglected due to the low density of the condensate and the short timescale of the atomic dynamics ($t_{dyn} \ll 2\pi/\omega_{ext}$).

In order to work with dimensionless coordinates, we introduce the following change of variables $t \to \frac{E_L}{\hbar}t$ and $x \to k_Lx$. We obtain

$$i \frac{d|\psi(t)\rangle}{dt} = \left( \hat{p}^2 - \frac{s}{2} \cos(\hat{x} + \varphi(t)) \right) |\psi(t)\rangle, \quad (2)$$

where $\hat{p} = -i \frac{\partial}{\partial x}$ and $\hat{x} = x$ in the position representation. Note that the study of this dimensionless system makes it possible to transfer the results of this paper to other experimental setups with different characteristics. The eigenvectors of the momentum operator are denoted by $|\phi_\alpha\rangle$, of wave function in the $x$ representation $\phi_\alpha(x) = \frac{1}{\sqrt{2\pi}} e^{i\alpha x}$ and of eigenvalue $\alpha$. Since the potential is periodic in $x$, the Bloch theorem states that $\alpha = n + q$, where $n \in \mathbb{Z}$ is a relative integer and $q \in [-0.5, 0.5]$ is the quasi-momentum. The periodicity of the potential also implies that the quasi-momentum is a constant of the motion during the control process. In a sub-Hilbert space characterized by a specific value of $q$, the state of the system can be expanded in the plane wave basis as follows

$$|\psi\rangle = \sum_{n \in \mathbb{Z}} c_{q,n} |\phi_{q+n}\rangle. \quad (3)$$
It is then straightforward to show that the coefficients $c_{q,n}$ are solutions of the following equation

$$i\dot{c}_{q,n} = (n+q)^2 c_{q,n} - \frac{s}{4} \left( e^{i\varphi(t)} c_{q,n-1} + e^{-i\varphi(t)} c_{q,n+1} \right). \quad (4)$$

We deduce that the Schrödinger equation can be written in matrix form as

$$i \frac{d|\psi(t)\rangle}{dt} = \hat{H} |\psi(t)\rangle = (\hat{H}_0 + \cos(\varphi(t)) \hat{H}_1 + \sin(\varphi(t)) \hat{H}_2) |\psi(t)\rangle, \quad (5)$$

with

$$|\psi(t)\rangle = \begin{pmatrix} \vdots \\ c_{q,n-1} \\ c_{q,n} \\ c_{q,n+1} \\ \vdots \end{pmatrix}, \quad (6)$$

$$\hat{H}_0 = \begin{pmatrix} \vdots & & & & \\ \cdots & 0 & ((n-1)+q)^2 & 0 & 0 & \cdots \\ \cdots & 0 & 0 & (n+q)^2 & 0 & \cdots \\ \cdots & 0 & 0 & 0 & ((n+1)+q)^2 & \cdots \\ & & & & & \ddots \end{pmatrix}, \quad (7)$$

and

$$\hat{H}_1 = \begin{pmatrix} \vdots & & & & \\ \cdots & -\frac{s}{4} & 0 & -\frac{s}{4} & 0 & \cdots \\ \cdots & 0 & -\frac{s}{4} & 0 & -\frac{s}{4} & \cdots \\ \cdots & 0 & 0 & -\frac{s}{4} & 0 & \cdots \\ & & & & \ddots & \ddots \\ & & & & & \ddots \end{pmatrix}, \quad \hat{H}_2 = \begin{pmatrix} \vdots & & & & \\ \cdots & -i\frac{s}{4} & 0 & i\frac{s}{4} & 0 & \cdots \\ \cdots & 0 & -i\frac{s}{4} & 0 & i\frac{s}{4} & \cdots \\ \cdots & 0 & 0 & -i\frac{s}{4} & 0 & \cdots \\ & & & & \ddots & \ddots \\ & & & & & \ddots \end{pmatrix}. \quad (8)$$

From a numerical point of view, the infinite-dimensional Hilbert space is truncated to a finite one such that $|n| \le n_{max}$, where $n_{max}$ is chosen as a function of the initial and target states in order to avoid edge effects in the computation. In the numerical simulations, we consider generally $n_{max} = 10$, which leads to a truncated space of dimension $2 \times n_{max} + 1 = 21$. The parameter $n_{max}$ is adjusted in each case depending on the target state.

Different theoretical and experimental implementations of QOC have shown the effectiveness of the procedure and the ability to achieve a variety of target states $|\psi_t\rangle$, with a very good match between theory and experiment [26, 46]. To this aim, the optimal phase is calculated from the application of GRAPE [9, 17] to the dynamical system (5). The goal of the optimal control problem is to maximize the fidelity $F = |\langle\psi_t|\psi(t_f)\rangle|^2$ at time $t_f$. The PMP leads to the Pontryagin Hamiltonian $H_p = \int \langle\chi(t)|\hat{H}|\psi(t)\rangle$, with $|\chi(t)\rangle$ the adjoint state, whose dynamics are also governed by Eq. (5), with the final condition
<figure><img src="image_3.png" /><figcaption>Figure 1: Example of three state-to-state transfers of a BEC from the initial state $|\phi_{0+0}\rangle$. The first, second and third rows (from top to bottom) correspond respectively to the target states $|\phi_{0+2}\rangle$, $\frac{1}{\sqrt{3}}(|\phi_{0-2}\rangle + |\phi_{0+0}\rangle + |\phi_{0+2}\rangle)$ and $|g(0, 0, 1/3)\rangle$. The columns (a), (b) and (c) represent respectively the corresponding controls $\varphi(t)$, the final population distribution in the momentum basis and the probability distribution in position at the final time.</figcaption></figure>

$\chi_0\langle\psi_t|\psi(t_f)\rangle|\psi_t\rangle$, where $\chi_0$ is the dual variable of the cost. $\chi_0$ is set to $1/2$ in the numerical simulation. We refer the interested reader to the recent tutorial [9] for a comprehensive introduction to the PMP for quantum optimal control.

Using the maximization condition of the PMP, the control is iteratively improved as

$$ \varphi(t) \mapsto \varphi(t) + \epsilon \Im \left( \langle \chi(t) | \left( -\sin(\varphi(t)) \hat{H}_1 + \cos(\varphi(t)) \hat{H}_2 \right) | \psi(t) \rangle \right), \quad (9) $$

where $\epsilon$ is a small positive parameter. The states $|\psi(t)\rangle$ and $|\chi(t)\rangle$ are respectively propagated forward and backward in time from their initial and final conditions. Figure 1 shows three numerical examples of state-to-state transfers. For the three examples, the initial state of the BEC is $|\phi_{0+0}\rangle$. The target states are described in the basis $|\phi_{0+n}\rangle$ with $q = 0$. They correspond either to a unique momentum state, a superposition of momentum states or to a Gaussian or a squeezed state (see [46] for the mathematical definition of these states). We denote by $|g(x_c, p, \xi)\rangle$ a squeezed state centered in $(x_c, p_c)$ with the squeezing parameter $\xi$. The following numerical parameters are used in the numerical simulation: $s = 5$, $n_{max} = 10$ and $t_f = 7.6$ (of the order of 150 $\mu$s in real units). More precisely, the target states are then the momentum state $|\phi_{0+2}\rangle$, the superposition $1/\sqrt{3}(|\phi_{0-2}\rangle + |\phi_{0+0}\rangle + |\phi_{0+2}\rangle)$, and the centered squeezed state $|g(x_c = 0, p_c = 0, \xi = 1/3)\rangle$.

### 3 Optimal control in the non-linear case

Although most experiments to date present a very good match with numerical simulations from the linear model system given in Eq. (2), improvement can be
made by taking into account the atom interactions which have been neglected in a first step due to the low density of the condensate. This interaction can be modeled at the mean-field level by the Gross-Pitaevskii equation, i.e. by an effective term in the Hamiltonian (2), proportional to the square modulus of the wave function $|\psi(x,t)|^2$ [60, 61]. This additional term makes the numerical resolution of the problem more complex. Indeed the non-linear interaction term does not have a simple state-independent matrix representation in the basis ($|\phi_{q+n}\rangle$)$_{n \in \mathbb{Z}}$. We present in this section an efficient numerical method which allows rapid propagation of the dynamics, taking into account the interactions between the atoms of the condensate and allows the application of an optimal control algorithm.

## 3.1 The model system

We consider the dimensionless equation (2) which describes the dynamics of a BEC in position representation. The interaction between the atoms can be treated at first order as a mean field effect, leading to an additional non-linear term $|\psi|^2$ and to the Gross-Pitaevskii equation,

$$i \frac{\partial \psi(x,t)}{\partial t} = \left( -\frac{\partial^2}{\partial x^2} - \frac{s}{2} \cos(x + \varphi(t)) + \beta |\psi(x,t)|^2 \right) \psi(x,t). \quad (10)$$

The weight of the non-linear term is given by the dimensionless parameter $\beta$. In the experiments, this parameter is of the order of $0.5 \pm 0.2$ [9]. This value depends on the number of atoms and on the frequency of the hybrid trap.

Expanding the state $|\psi\rangle$ over the eigenvectors of the momentum operator $|\psi\rangle = \sum_{n \in \mathbb{Z}} c_{q,n} |\phi_{q+n}\rangle$, the non-linear term $|\psi|^2$ can be expressed in position representation as,

$$|\psi|^2 = \psi(x,t)^* \psi(x,t) = \left( \frac{1}{\sqrt{2\pi}} \sum_{m \in \mathbb{Z}} c_{q,m}^* e^{-i(q+m)x} \right) \left( \frac{1}{\sqrt{2\pi}} \sum_{l \in \mathbb{Z}} c_{q,l} e^{i(q+l)x} \right)$$

which leads to

$$|\psi|^2 = \frac{1}{2\pi} \sum_{m,l \in \mathbb{Z}} c_{q,m}^* c_{q,l} e^{i(q+l-m)x}.$$

The dynamics of the coefficients $c_{q,n}$ are then given by

$$i \dot{c}_{q,n} = (n+q)^2 c_{q,n} - \frac{s}{4} \left( e^{i\varphi(t)} c_{q,n-1} + e^{-i\varphi(t)} c_{q,n+1} \right) + \frac{\beta}{2\pi} \sum_{m,l \in \mathbb{Z}} c_{q,m}^* c_{q,l} c_{q,n+m-l}. \quad (11)$$

Unlike the Schrödinger equation (5), an analytical matrix representation to this non-linear problem cannot be given here. This means that the simple matrix exponential method used earlier cannot be used anymore to propagate the dynamics. A first option consists in applying the Runge-Kutta of order 4 approach (RK4). However, for large dimensional systems, this method can be very costly in terms of calculation time. Typically a propagation of Eq. (11) takes about 15 seconds, for a constant control $\varphi(t) = 0$ (see Fig. 2). Numerical simulations have been performed on a standard laptop. Optimal control, which requires a hundred propagation, is therefore difficult to apply with this approach, as the design of a single control would require a computational time of the
order of one hour. In order to speed up the optimization process, we propose in this paper to combine optimal control algorithms with the FBR-DVR bases as described in Sec. 3.2.

## 3.2 The FBR-DVR bases

An efficient approach to express the non-linear term of Eq. (10) in matrix form is the FBR-DVR (or pseudo-spectral) approach [55–59]. This method consists in using two different bases for the matrix representation of the operators. In our case, the FBR basis (for Finite Basis Representation) is the basis of the eigenvectors $|\phi_{q+n}\rangle$ of the momentum operator. In this basis, it is straightforward to express the operators $\hat{p}$ and $\cos(\hat{x})$. The DVR basis (for Discrete Variable Representation) corresponds to a basis built from the discretization of a continuous parameter, here the position $x$. In this basis, a matrix representation of the non-linear term $|\psi(x, t)|^2$ can be given. A unitary transformation allows to go from one basis to the other, and therefore to express all the operators in the FBR basis for example. The time propagation of the system can then be done by the matrix exponential approach.

We apply this method to Eq. (10). The FBR basis corresponds to the basis of the eigenvectors $|\phi_{q+n}\rangle$ of the momentum operator. From a numerical point of view, we work in a finite dimensional sub-space $\mathcal{H}_N$ of $\mathcal{H}$ of dimension $N = 2n_{max} + 1$, with $n_{max} \in \mathbb{N}$ and $-n_{max} \le n \le n_{max}$. The potential energy is a periodic spatial function, of period $2\pi$. We choose the interval $[0, 2\pi]$ for the position $x$. The DVR basis is built from the discretization of the variable $x \in [0, 2\pi]$, such that $x_j = \frac{2\pi}{N}j$ with $0 \le j < N$. A scalar product is given by,

$$ \forall \Psi, \Phi \in \mathcal{H}_N, \int_0^{2\pi} \Psi^*(x)\Phi(x)dx = \frac{2\pi}{N} \sum_{j=0}^{N-1} \Psi^*(x_j)\Phi(x_j). $$

From the orthogonality of functions $\phi_{q+n}$ we deduce that we have approximately:

$$ \frac{2\pi}{N} \sum_{j=0}^{N-1} \phi_{q+m}^{*}(x_j)\phi_{q+n}(x_j) = \delta_{m,n}. \qquad (12) $$

From the definition of the plane wave functions, we have also for $n_{max}$ large enough:

$$ \frac{2\pi}{N} \sum_{n=-n_{max}}^{n_{max}} \phi_{q+n}(x_i)\phi_{q+n}^{*}(x_j) = \delta_{i,j}. \qquad (13) $$

We introduce the matrix $\hat{R}$ defined as

$$ \hat{R} = \begin{pmatrix} R_{0,-n_{\max}} & R_{0,-n_{\max}+1} & \cdots \\ \vdots & \vdots & \cdots \\ R_{N-1,-n_{\max}} & R_{N-1,n_{\max}} & \cdots \end{pmatrix}, \qquad (14) $$

with matrix elements $R_{j,n}$ given by

$$ R_{j,n} = \frac{1}{\sqrt{N}} e^{i \frac{2\pi}{N} (q+n) j}. \qquad (15) $$
Equations (12) and (13) are equivalent to $\hat{R}^\dagger \hat{R} = \hat{I}$ and $\hat{R} \hat{R}^\dagger = \hat{I}$. The matrix $\hat{R}$ is thus a unitary matrix. The DVR basis ($|u_j\rangle$)$_{j \in [0,N-1]}$ is then built from the following kets as

$$|u_j\rangle = \sum_{n=n_{\max}}^{-n_{\max}} \hat{R}_{jn}^* |\phi_n\rangle. \qquad (16)$$

This basis connects a basis of $\mathcal{H}_N$ to the spatial discretization basis. We deduce that an operator $\hat{W}$ depending on $x$ has a diagonal representation in the DVR basis given by

$$\langle u_i | \hat{W} | u_j \rangle \approx W(x_j) \delta_{i,j}. \qquad (17)$$

The matrix $\hat{R}$ allows to express a known operator in the DVR basis to the FBR basis, and vice versa

$$\hat{W}^{\text{FBR}} = \hat{R}^{\dagger} \hat{W}^{\text{DVR}} \hat{R}. \qquad (18)$$

Since the square modulus of the wave function, $|\psi(x)|^2$, depends on the position $x$ as shown in Eq (11), it is straightforward to obtain its matrix representation in the DVR basis, and then in the FBR basis. The exponential matrix approach can then be used to propagate the dynamics. We denote by $\hat{G}^{\text{DVR}}$ the matrix of $|\psi(x)|^2$ in the DVR basis,

$$\hat{G}^{\text{DVR}} = \begin{pmatrix} G_0 & & & \\ & G_1 & & \\ & & \ddots & \\ & & & G_{N-1} \end{pmatrix}, \qquad (19)$$

where

$$G_j = \left| \frac{1}{\sqrt{2\pi}} \sum_{n=-n_{\max}}^{n_{\max}} c_{q,n} e^{i \frac{2\pi}{N} (q+n)j} \right|^2. \qquad (20)$$

The parameters $c_{q,n}$ are the coefficients of the state in the FBR basis. Using the notations introduced in Sec. 2 and the relations (14) and (19), Eq. (10) can be written in matrix form in the FBR basis as follows

$$i \frac{d|\psi\rangle}{dt} = \hat{H}_{GP} |\psi\rangle = \left( \hat{H}_0 + \cos(\varphi(t)) \hat{H}_1 + \sin(\varphi(t)) \hat{H}_2 + \beta \hat{R}^\dagger \hat{G}^{\text{DVR}} \hat{R} \right) |\psi\rangle, \qquad (21)$$

with $\hat{H}_{GP}$ the Gross-Pitaevskii Hamiltonian.

For comparison, the propagation for a constant control is here of the order of 0.6 seconds, i.e. a total of around 2 minutes for the calculation time of the optimal control. Figure 2 shows the time evolution of the projection of the state $|\psi\rangle$ on $|\phi_{0+0}\rangle$, calculated from the matrix exponential approach and the RK4 method. The initial state is the squeezed state $|g(0, 0, 1/2)\rangle$ and the phase $\varphi(t) = 0$. For $\beta = 1$, we observe that the two methods give the same result.

## 3.3 The non-linear GRAPE

We show in this section how to extend the standard GRAPE algorithm to this non-linear case. The corresponding algorithm is called non-linear GRAPE.
<figure><img src="image_4.png" /><figcaption>Figure 2: Time evolution of $|\langle\psi(t)|\phi_{0+0}\rangle|^2$ for a zero control, $\varphi = 0$. The initial state is $|\psi(0)\rangle = |g(0, 0, 1/2)\rangle$. The blue and orange solid lines represent respectively the propagation of the state with the matrix exponential approach for $\beta = 0$ and $\beta = 1$. The dashed black line corresponds to a dynamic computed with the RK4 method for $\beta = 1$.</figcaption></figure>

We consider the state-to-state transfer defined with the fidelity $F = |\langle\psi_t|\psi(t_f)\rangle|^2$, with the dynamics of Eq. (10). In this case, the Pontryagin Hamiltonian $H_p$ can be written as

$$H_p = \mathfrak{F}[\langle\chi|\hat{H}_{GP}|\psi\rangle], \quad (22)$$

with $\hat{H}_{GP}$, the Gross-Pitaevskii Hamiltonian in position representation,

$$\hat{H}_{GP} = -\frac{\partial^2}{\partial x^2} - \frac{s}{2} \cos(x + \varphi(t)) + \beta |\psi(x, t)|^2. \quad (23)$$

The time evolution of the adjoint state is given by the relation, $\dot{\chi}(x, t) = -\partial_\psi H_p$, which leads to

$$i \frac{\partial \chi}{\partial t} = \left( -\frac{\partial^2}{\partial x^2} - \frac{s}{2} \cos(x + \varphi(t)) + 2\beta |\psi|^2 \right) \chi - \beta \psi^2 \chi^*. \quad (24)$$

Note that the nonlinearity of the Gross-Pitaeski equation breaks the symmetry between the state and the adjoint state which are no longer solutions of the same differential equation. The transversality condition on the adjoint state $\chi(x, t_f) = \chi_0 \partial_\psi F$ yields

$$\chi(x, t_f) = 2\chi_0 \langle\psi_t|\psi(t_f)\rangle \psi_t, \quad (25)$$

with $\chi_0 = 1/2$ in the numerical simulations. Finally, the maximization condition of the Pontryagin Hamiltonian leads to the iterative procedure of GRAPE

$$\varphi(t) \mapsto \varphi(t) + \epsilon \mathfrak{F}[\langle\chi|\frac{\partial \hat{H}_{GP}}{\partial \varphi}|\psi\rangle]$$

where $\epsilon$ is a small positive parameter. The state and the adjoint state are respectively propagated forward and backward in time. The differential equation (24)
associated with the adjoint state can be written as

$$i \frac{\partial \chi}{\partial t} = \left( -\frac{\partial^2}{\partial x^2} - \frac{s}{2} \cos(x + \varphi(t)) + \beta |\psi|^2 \right) \chi - 2i\beta \Im[\chi^* \psi] \psi. \quad (26)$$

The backward propagation can be done in matrix form with the final conditions $|\psi(t_f)\rangle = |\psi_f\rangle$ and $|\chi(t_f)\rangle = \langle\psi_t|\psi(\tau_f)\rangle|\psi_t\rangle$, with the following extended system

$$i \begin{pmatrix} |\dot{\psi}\rangle \\ |\dot{\chi}\rangle \end{pmatrix} = \begin{pmatrix} \hat{H}_{\text{GP}} & 0 \\ -2i\beta \hat{R}^\dagger \hat{I}^\text{DVR} \hat{R} & \hat{H}_{\text{GP}} \end{pmatrix} \begin{pmatrix} |\psi\rangle \\ |\chi\rangle \end{pmatrix}, \quad (27)$$

where $\hat{H}_{\text{GP}}$ is the Gross-Pitaevskii Hamiltonian in matrix form given by Eq. (21), and $\hat{I}^\text{DVR}$ is the diagonal matrix associated with the term $\Im[\chi^* \psi]$ in the DVR basis. The matrix elements of $\hat{I}^\text{DVR}$ are $I_j = \Im\left[\frac{1}{2\pi} \sum_{m,n} d_{q,m}^* c_{q,n} e^{i \frac{2\pi}{N} (q+n-m)j}\right]$, $d_{q,n}$ being the coefficients of the adjoint state $|\chi\rangle$ in the FBR basis, $-n_{\text{max}} \le n, m \le n_{\text{max}}$.

As an illustrative example, we consider the transfer from the state $|\psi(0)\rangle = |\phi_{0+0}\rangle$ to the squeezed state $|g(0, 0, 3/2)\rangle$. The control time is set to $t_f = 150 \ \mu s$, the lattice depth to $s = 5$, and the quasi-momentum to $q = 0$. A first control is optimized for $\beta = 0$ as represented in Fig. 3. The dynamics are then propagated for this control with a non-linear coefficient $\beta$ going from 0 to 1, $\beta \in [0, 1]$. Figure 4 shows the evolution of the fidelity with respect to the nonlinear parameter $\beta$. The grey area indicates the experimental uncertainty on this parameter. We observe that for $\beta = 0.7$, the upper bound of the uncertainty interval, the fidelity is only of the order of 0.93, which means that the target state is not reached by the system. A second control is optimized by using the non-linear version of GRAPE for $\beta = 0.5$. In this case, the fidelity is equal to 0.994 for the same control time. Figure 3 depicts the different populations in the momentum eigenbasis. The role of $\beta$ on the dynamics depends on the initial and target states considered for the optimal control.

# 4 Optimal control of a BEC in a two-dimensional optical lattice

In this section we show how to apply the GRAPE algorithm to a two-dimensional optical lattice. We assume here that the non-linearity of the Gross-Pitaevski equation can be neglected.

Before presenting the modeling of an experiment with a 2D lattice, it may be useful to recall the concept of direct and reciprocal lattices, as well as Bloch's theorem in the multidimensional case. The direct lattice is a set of vectors describing the nodes of the lattice, i.e. here the positions where the potential is minimum. We denote by $\mathbf{a}_i$ the vectors of the basis of the direct lattice and $N$ the corresponding dimension. A node is an element of this set,

$$\mathcal{D} = \left\{ \sum_{i=1}^{N} j_i \mathbf{a}_i, j_i \in \mathbb{Z} \right\}. \quad (28)$$

The reciprocal lattice is described by the set,

$$\mathcal{R} = \left\{ \sum_{i=1}^{N} j_i \mathbf{b}_i, j_i \in \mathbb{Z} \right\} \quad (29)$$
<figure><img src="image_5.png" /><figcaption>Figure 3: (a)) Time evolution of the optimal controls for $\beta = 0$ (red line) and $\beta = 0.7$ (green line). (b) Populations $|c_{0,n}|^2$ at the final time. The blue, green and red columns correspond respectively to the populations of the target state $|g(0, 0, 3/2)\rangle$, of the final state for a control optimized with $\beta = 0.7$ and of the final state for a control with $\beta = 0$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4: Evolution of the fidelity $F$ for the target state $|g(0, 0, 3/2)\rangle$ as a function of $\beta$. The blue and orange lines correspond respectively to $\beta = 0$ and $\beta = 0.5$.</figcaption></figure>

where the vectors $\mathbf{b}_i$ are defined from the relation,

$$ \mathbf{a}_i \cdot \mathbf{b}_j = 2\pi\delta_{ij}, \qquad (30) $$

with $\mathbf{a}_i \in \mathcal{D}$ and $\mathbf{b}_j \in \mathcal{R}$. Using the Bloch's theorem, it can be shown that the wave function $\psi(\mathbf{r})$ solution of the time-independent Schrödinger equation,

$$ \left(-\frac{\hbar^2}{2m} \frac{\partial^2}{\partial \mathbf{r}^2} + V(\mathbf{r})\right) \psi(\mathbf{r}) = E\psi(\mathbf{r}), \qquad (31) $$

where the potential $V$ is periodic, $V(\mathbf{r} + \mathbf{R}) = V(\mathbf{r})$ for any $\mathbf{R} \in \mathcal{D}$ and $\mathbf{r} = \sum_{i=1}^{N} r_i \mathbf{a}_i$ is the position of the system, can be expressed as

$$ \psi_{\mathbf{q}}(\mathbf{r}) = e^{i\mathbf{q}\cdot\mathbf{r}}u_{\mathbf{q}}(\mathbf{r}), \qquad (32) $$

with $\mathbf{q} = \sum_{i=1}^{N} q_i \mathbf{b}_i$ the two-dimensional quasi-momentum and $u_{\mathbf{q}}(\mathbf{r})$ a periodic function in $\mathcal{D}$.

## 4.1 The model system

A one-dimensional optical lattice can be generated by the superposition of two lasers. A two-dimensional optical lattice requires three or more lasers. A large variety of two-dimensional configurations can be realized experimentally [62]. As an illustrative example, we consider the following configuration, three lasers with the same angular frequency $\omega$ which propagate in the plane $(\mathbf{e}_x, \mathbf{e}_y)$ of the laboratory frame with a linear polarization along the $z$-axis. The wave vectors are defined in this plane by

$$ \mathbf{k}_1 = k(1, 0) $$

$$ \mathbf{k}_2 = \frac{k}{2}(-1, \sqrt{3}) $$

$$ \mathbf{k}_3 = -\frac{k}{2}(1, \sqrt{3}), $$
where $k = \frac{\omega}{c}$ is the wave number. The angle between the vectors is $2\pi/3$. The total electric field is then given by

$$ \mathbf{E} = E_0 \sum_i \cos (\mathbf{k}_i \cdot \mathbf{r} - \omega t + \varphi_i) \mathbf{e}_z, \quad (33) $$

where $E_0$ is the amplitude of the electric field, $\mathbf{r}$ the position vector and $\varphi_i$ is the phase of the laser $i$, with $i = 1, 2, 3$. The interaction between the electric field and the atom is described by the Hamiltonian $\hat{H}_I = -\mu \cdot \mathbf{E}$, where $\mu$ is the induced dipole moment. At first order, the dipole moment is given by $\mu = -\frac{\alpha}{2}\mathbf{E}$, where $\alpha$ is the polarizability of the atom. This leads to $\hat{H}_I = -\frac{\alpha}{2}\mathbf{E} \cdot \mathbf{E}$. The electric field is not resonant with the transition frequencies of the atom, and the period associated with the $\omega$ frequency is much shorter than the characteristic time of the experiment, $T_{dyn}$. The interaction Hamiltonian can be averaged to only consider the long-time dynamics. The BEC is then subjected to the following dipole potential

$$ V = \langle \hat{H}_I \rangle = -\frac{3}{4} \alpha E_0^2 - \frac{\alpha}{2} E_0^2 \left[ \begin{aligned} & \cos ((\mathbf{k}_1 - \mathbf{k}_2) \cdot \mathbf{r} + \varphi_{12}) \\ & + \cos ((\mathbf{k}_2 - \mathbf{k}_3) \cdot \mathbf{r} + \varphi_{23}) \\ & + \cos ((\mathbf{k}_3 - \mathbf{k}_1) \cdot \mathbf{r} + \varphi_{31}) \end{aligned} \right], \quad (34) $$

where $\varphi_{jk} = \varphi_j - \varphi_k$, and $\langle \cdot \rangle$ denotes the time average over one period $\tau$, such that $\frac{2\pi}{\omega} \gg \tau \gg T_{dyn}$. Finally, we deduce that this average potential governs the 2D-dynamics (31), with $\mathbf{r} = (x, y)$.

The direct lattice nodes correspond to the minima of the potential energy for $\varphi_{jk} = 0$. We have

$$ \begin{aligned} (\mathbf{k}_1 - \mathbf{k}_2) \cdot \mathbf{r} &= 0 \pmod{2\pi} \\ (\mathbf{k}_2 - \mathbf{k}_3) \cdot \mathbf{r} &= 0 \pmod{2\pi} \\ (\mathbf{k}_3 - \mathbf{k}_1) \cdot \mathbf{r} &= 0 \pmod{2\pi} \end{aligned} \quad (35) $$

A node of the direct lattice belongs to the set $\mathcal{D}$ (28) with,

$$ \begin{aligned} \mathbf{a}_1 &= -\frac{4\pi}{3k^2}\mathbf{k}_3, \\ \mathbf{a}_2 &= \frac{4\pi}{3k^2}\mathbf{k}_2. \end{aligned} $$

A node of the reciprocal lattice is an element of the set $\mathcal{R}$ (29) with

$$ \begin{aligned} \mathbf{b}_1 &= \mathbf{k}_1 - \mathbf{k}_3, \\ \mathbf{b}_2 &= \mathbf{k}_2 - \mathbf{k}_1. \end{aligned} $$

Since the potential is periodic, Bloch's theorem applies, and eigenvectors are as in Eq. (32). Within a given sub-Hilbert space of quasi-momentum $\mathbf{q}$, a generic wave function can be written as:

$$ \psi(\mathbf{r}) = \sum_{m,n \in \mathbb{Z}} c_{\mathbf{q},m,n} e^{i\mathbf{q} \cdot \mathbf{r}} e^{i(m\mathbf{b}_1 + n\mathbf{b}_2) \cdot \mathbf{r}}. \quad (36) $$

We denote by $\phi_{m,n} = e^{i(m\mathbf{b}_1 + n\mathbf{b}_2) \cdot \mathbf{r}}$ the functions of the Fourier basis. The coefficients $c_{\mathbf{q},m,n}$ are the coefficients of the state in this basis. In order to
simplify the study, we consider below that the quasi-momentum is zero, which is a reasonable experimental assumption [26]. Plugging the expression of $\psi(\mathbf{r}, t)$ of Eq. (36) into Eq. (31), we obtain,

$$i\hbar \sum_{m,n \in \mathbb{Z}} \dot{c}_{m,n} e^{i(m\mathbf{b}_1 + n\mathbf{b}_2)\mathbf{r}} = \sum_{m,n \in \mathbb{Z}} \frac{\hbar^2}{2m} c_{m,n} (m\mathbf{b}_1 + n\mathbf{b}_2)^2 e^{i(m\mathbf{b}_1 + n\mathbf{b}_2)\mathbf{r}} \\ - \sum_{m,n \in \mathbb{Z}} \frac{V_0}{4} c_{m,n} \left[ e^{i((\mathbf{k}_1 - \mathbf{k}_2)\cdot\mathbf{r} + \varphi_{12})} e^{i(m\mathbf{b}_1 + n\mathbf{b}_2)\mathbf{r}} + e^{-i((\mathbf{k}_1 - \mathbf{k}_2)\cdot\mathbf{r} + \varphi_{12})} e^{i(m\mathbf{b}_1 + n\mathbf{b}_2)\mathbf{r}} \right] \\ - \sum_{m,n \in \mathbb{Z}} \frac{V_0}{4} c_{m,n} \left[ e^{i((\mathbf{k}_2 - \mathbf{k}_3)\cdot\mathbf{r} + \varphi_{23})} e^{i(m\mathbf{b}_1 + n\mathbf{b}_2)\mathbf{r}} + e^{-i((\mathbf{k}_2 - \mathbf{k}_3)\cdot\mathbf{r} + \varphi_{23})} e^{i(m\mathbf{b}_1 + n\mathbf{b}_2)\mathbf{r}} \right] \\ - \sum_{m,n \in \mathbb{Z}} \frac{V_0}{4} c_{m,n} \left[ e^{i((\mathbf{k}_3 - \mathbf{k}_1)\cdot\mathbf{r} + \varphi_{31})} e^{i(m\mathbf{b}_1 + n\mathbf{b}_2)\mathbf{r}} + e^{-i((\mathbf{k}_3 - \mathbf{k}_1)\cdot\mathbf{r} + \varphi_{31})} e^{i(m\mathbf{b}_1 + n\mathbf{b}_2)\mathbf{r}} \right],$$

where the constant term $-\frac{3}{4}\alpha E_0^2$ of the potential $V$ is removed, and $V_0 = \alpha E_0^2$ is the lattice depth. From the relations, $\mathbf{k}_1 - \mathbf{k}_2 = -\mathbf{b}_2$, $\mathbf{k}_2 - \mathbf{k}_3 = \mathbf{b}_1 + \mathbf{b}_2$, $\mathbf{k}_3 - \mathbf{k}_1 = -\mathbf{b}_1$, $\mathbf{b}_1 \cdot \mathbf{b}_1 = \mathbf{b}_2 \cdot \mathbf{b}_2 = 3k^2$ and $\mathbf{b}_1 \cdot \mathbf{b}_2 = -\frac{3}{2}k^2$, we arrive at

$$i\hbar \sum_{m,n \in \mathbb{Z}} \dot{c}_{m,n} e^{i(m\mathbf{b}_1 + n\mathbf{b}_2)\mathbf{r}} = \sum_{m,n \in \mathbb{Z}} \frac{3k^2\hbar^2}{2m} c_{m,n} (m^2 + n^2 - mn) e^{i(m\mathbf{b}_1 + n\mathbf{b}_2)\mathbf{r}} \\ - \sum_{m,n \in \mathbb{Z}} \frac{V_0}{4} c_{m,n} \left[ e^{i\varphi_{12}} e^{i(m\mathbf{b}_1 + (n-1)\mathbf{b}_2)\mathbf{r}} + e^{-i\varphi_{12}} e^{i(m\mathbf{b}_1 + (n+1)\mathbf{b}_2)\mathbf{r}} \right] \\ - \sum_{m,n \in \mathbb{Z}} \frac{V_0}{4} c_{m,n} \left[ e^{i\varphi_{23}} e^{i((m+1)\mathbf{b}_1 + (n+1)\mathbf{b}_2)\mathbf{r}} + e^{-i\varphi_{23}} e^{i((m-1)\mathbf{b}_1 + (n-1)\mathbf{b}_2)\mathbf{r}} \right] \\ - \sum_{m,n \in \mathbb{Z}} \frac{V_0}{4} c_{m,n} \left[ e^{i\varphi_{31}} e^{i((m-1)\mathbf{b}_1 + n\mathbf{b}_2)\mathbf{r}} + e^{-i\varphi_{31}} e^{i((m+1)\mathbf{b}_1 + n\mathbf{b}_2)\mathbf{r}} \right].$$

Projecting this equation onto the state $\phi_{m,n}$ and introducing the dimensionless time $\tau = \frac{3\hbar k^2}{2m}t$ and lattice depth $s = 2mV_0/(3\hbar^2 k^2)$, we obtain the differential equation governing the dynamics of the coefficients $c_{m,n}$,

$$i\dot{c}_{m,n} = (m^2 + n^2 - mn)c_{m,n} - \frac{s}{4} [e^{i\varphi_{12}}c_{m,n+1} + e^{-i\varphi_{12}}c_{m,n-1} + e^{i\varphi_{23}}c_{m-1,n-1} + e^{-i\varphi_{23}}c_{m+1,n+1} + e^{i\varphi_{31}}c_{m+1,n} + e^{-i\varphi_{31}}c_{m-1,n}] \quad (37)$$

## 4.2 Optimal control problem

At this point, Eq. (37) can be put in matrix from. We have

$$i \frac{d|\psi(t)\rangle}{dt} = (\hat{H}_0 + e^{i\varphi_{12}}\hat{H}_{12}^+ + e^{-i\varphi_{12}}\hat{H}_{12}^- + e^{i\varphi_{23}}\hat{H}_{23}^+ + e^{-i\varphi_{23}}\hat{H}_{23}^- + e^{i\varphi_{31}}\hat{H}_{31}^+ + e^{-i\varphi_{31}}\hat{H}_{31}^-) |\psi(t)\rangle,$$

with

$$\left\{
\begin{aligned}
\hat{H}_0 &= \sum_{m,n \in \mathbb{Z}} (m^2 + n^2 - mn) |\phi_{m,n}\rangle\langle\phi_{m,n}| \\
\hat{H}_{12}^+ &= -\frac{s}{4} \sum_{m,n \in \mathbb{Z}} |\phi_{m,n}\rangle\langle\phi_{m,n+1}| \\
\hat{H}_{12}^- &= -\frac{s}{4} \sum_{m,n \in \mathbb{Z}} |\phi_{m,n}\rangle\langle\phi_{m,n-1}| \\
\hat{H}_{23}^+ &= -\frac{s}{4} \sum_{m,n \in \mathbb{Z}} |\phi_{m,n}\rangle\langle\phi_{m-1,n-1}| \\
\hat{H}_{23}^- &= -\frac{s}{4} \sum_{m,n \in \mathbb{Z}} |\phi_{m,n}\rangle\langle\phi_{m+1,n+1}| \\
\hat{H}_{31}^+ &= -\frac{s}{4} \sum_{m,n \in \mathbb{Z}} |\phi_{m,n}\rangle\langle\phi_{m+1,n}| \\
\hat{H}_{31}^- &= -\frac{s}{4} \sum_{m,n \in \mathbb{Z}} |\phi_{m,n}\rangle\langle\phi_{m-1,n}|
\end{aligned}
\right.$$
From a numerical point of view, we use a Hilbert space of finite dimension such that $-M \le m \le M$ and $-N \le n \le N$. The dimension of this space is $d_H = (2M+1) \times (2N+1)$. In the numerical simulations, $M$ and $N$ have been set to 5. The GRAPE algorithm can then be applied to this case [9]. For each control $\varphi_{12}$, $\varphi_{23}$ and $\varphi_{31}$, the maximization condition given by the PMP leads to the iterative procedure of GRAPE

$$ \varphi_{\ell} \mapsto \varphi_{\ell} + \epsilon \Im \left[ \langle \chi(t) | t (\hat{H}_{\ell}^{+} - \hat{H}_{\ell}^{-}) | \psi(t) \rangle \right], $$

where $|\chi\rangle$ is the adjoint state and $\ell$ corresponds to the index 12, 23 or 31.

As an illustrative example, we consider the state-to-state transfer from the initial state $|\psi(0)\rangle = |0,0\rangle$ to a target state $|\psi_t\rangle = \frac{1}{\sqrt{2}}(|3,3\rangle + |-3,-3\rangle)$ with a control time $t_f = 250 \ \mu s$. The values of the other parameters are the same as in Sec. 2. Figure 5 represents the optimization result for this control problem. Note that the controls are not independent since $\varphi_{12}$ and $\varphi_{31}$ are equal. This observation also illustrates the role of the number of controls for the controllability of this system. Figure 6 shows this point for another target state. In this case, only two controls $\varphi_{2,3}$ and $\varphi_{3,1}$ are optimized, the third one is kept constant in time. Despite this constraint, the target state is reached with a very good precision. Finally, in Fig. 7, we consider a non-symmetric state in $m$ and $n$ with $|\varphi_t\rangle = \frac{1}{\sqrt{2}}(|\phi_{1,2}\rangle + |\phi_{-3,-1}\rangle)$. Again, our algorithm designs a very efficient control process with a fidelity above 0.99. We note that a broad range of states can be reached with just two variable phases, which raises the question of the controllability of this system, that is the relationship between the set of accessible states, the number of available controls and the lattice configuration.

# 5 Conclusion and prospective views

We have proposed two extensions of the gradient-based algorithm, GRAPE, for controlling a BEC in an optical lattice. We have shown how to adapt this algorithm to the non-linear case where the atomic interaction is not neglected and for a two-dimensional optical lattice. The two generalizations are supported by a mathematical analysis of the optimal control problem based on the PMP. Numerical examples show the effectiveness of the different procedures in experimentally realistic cases. We emphasize that such algorithms have the advantage of simplicity and general applicability, regardless of the structure of the optical lattice at two or three dimensions or the number of available controls. Thus, based on the material presented in this paper, GRAPE can be adapted to other geometric configurations.

This work describes in detail how to numerically solve the optimal control of a BEC in an optical lattice, but it also raises a number of issues that go beyond the scope of this paper. A first problem concerns the controllability of this system under experimental constraints and limitations on pulse amplitude and duration. A general problem is to describe the set of reachable sets for a given lattice configuration and physical constraints on the control parameters. Some mathematical results have shown the small-time controllability of this system in an ideal situation [63, 64]. However, such results do not use a directly implementable control strategy and consider pulses with a very large
<figure><img src="image_7.png" /><figcaption>Figure 5: Transfer from the state $|\phi_{0,0}\rangle$ to the target $\frac{1}{\sqrt{2}}(|\phi_{-3,-3}\rangle + |\phi_{3,3}\rangle)$; (a) Time evolution of the controls $\phi_{12}, \phi_{23}$ and $\phi_{31}$, (b) Populations of the initial (orange) and target (blue) states.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: Same as Fig. 5 but for the target state $\frac{1}{\sqrt{2}}(|\phi_{-1,-1}\rangle + |\phi_{3,3}\rangle)$. Only the controls $\phi_{3,1}$ and $\phi_{2,3}$ are optimized.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 7: Same as Fig. 5 but for the target state $\frac{1}{\sqrt{2}}(|\phi_{1,2}\rangle + |\phi_{-3,-1}\rangle)$. Only the controls $\phi_{3,1}$ and $\phi_{2,3}$ are optimized.</figcaption></figure>

amplitude. Another interesting question is to evaluate the role of the nonlinearity on the optimal solution. It has been shown in [65] that the quantum speed limit [66] grows with the nonlinearity strength of the BEC. An intriguing question is to test this conjecture in our system using time-optimal control protocols. If this statement is verified, it could be of utmost importance to minimize the preparation time of specific states by playing with the nonlinearity of the BEC dynamics. This work assumes that the Hamiltonian parameters of the system are exactly known. It is not the case in practice where the experimental parameters are estimated within a given range. This limitation could be partly avoided by generalizing to this system the robust control pulses known for two-level quantum systems [67–70]. An optimal control can be designed in this case by considering the simultaneous control of a set of quantum systems characterized by a different value of the parameter. In this study, we only consider state-to-state transfer, but it will be interesting to extend this approach to quantum gates [71, 72]. Recently, a rigorous foundation of quantum computing in the nonlinear case has been developed [73]. In this direction, a first goal will be to show that the non-linear GRAPE algorithm described in this work can be used to implement quantum gates. Work is underway in our group to answer such questions both theoretically and experimentally.

## Funding

We thank the support from the Erasmus Mundus Master QuanTeem (Project Number: 101050730), the project QuanTEdu-France (ANR-22-CMAS-0001) on
quantum technologies and the ANR project QuCoBEC (ANR-22-CE47-0008-02). This project has received financial support from the CNRS through the MITI interdisciplinary programs (QUSPIDE project).

## References

[1] Acín, A., Bloch, I., Buhrman, H., Calarco, T., Eichler, C., Eisert, J., et al. (2018). The quantum technologies roadmap: a european community view. *New J. Phys.* 20, 080201

[2] Glaser, S. J., Boscain, U., Calarco, T., Koch, C. P., Köckenberger, W., Kosloff, R., et al. (2015). Training Schrödinger's cat: quantum optimal control. Strategic report on current status, visions and goals for research in Europe. *European Physical Journal D* 69, 279

[3] Brif, C., Chakrabarti, R., and Rabitz, H. (2010). Control of quantum phenomena: Past, present, and future. *New J. Phys.* 12, 075008

[4] Koch, C. P., Boscain, U., Calarco, T., Dirr, G., Filipp, S., Glaser, S. J., et al. (2022). Quantum optimal control in quantum technologies. strategic report on current status, visions and goals for research in europe. *EPJ Quantum Technology* 9, 19. doi:10.1140/epjqt/s40507-022-00138-x

[5] Werninghaus, M., Egger, D. J., Roy, F., Machnes, S., Wilhelm, F. K., and Filipp, S. (2021). Leakage reduction in fast superconducting qubit gates via optimal control. *NPJ Quantum Information* 7, 14

[6] Abdelhafez, M., Baker, B., Gyenis, A., Mundada, P., Houck, A. A., Schuster, D., et al. (2020). Universal gates for protected superconducting qubits using optimal control. *Phys. Rev. A* 101, 022321

[7] Wilhelm, F. K., Kirchhoff, S., Machnes, S., Wittler, N., and Sugny, D. (2020). An introduction into optimal control for quantum technologies. *Lecture notes for the 51st IFF Spring School*

[8] Hohenester, U., Rekdal, P. K., Borzì, A., and Schmiedmayer, J. (2007). Optimal quantum control of bose-einstein condensates in magnetic microtraps. *Phys. Rev. A* 75, 023602

[9] Ansel, Q., Dionis, E., Arrouas, F., Peaudecerf, B., Guérin, S., Guéry-Odelin, D., et al. (2024). Introduction to theoretical and experimental aspects of quantum optimal control. *Journal of Physics B: Atomic, Molecular and Optical Physics* 57, 133001

[10] van Frank, S., Bonneau, M., Schmiedmayer, J., Hild, S., Gross, C., Cheneau, M., et al. (2016). Optimal control of complex atomic quantum systems. *Sci. Rep.* 6, 34187

[11] Koch, C. P., Lemeshko, M., and Sugny, D. (2019). Quantum control of molecular rotation. *Rev. Mod. Phys.* 91, 035005
[12] Rembold, P., Oshnik, N., Müller, M. M., Montangero, S., Calarco, T., and Neu, E. (2020). Introduction to quantum optimal control for quantum sensing with nitrogen-vacancy centers in diamond. *AVS Quantum Science* 2

[13] Pontryagin, L. S., Boltianski, V., Gamkrelidze, R., and Mitchtenko, E. (1962). *The Mathematical Theory of Optimal Processes* (John Wiley and Sons, New York)

[14] Boscain, U., Sigalotti, M., and Sugny, D. (2021). Introduction to thepontryagin maximum principle for quantum optimal control. *PRX Quantum* 2, 030203

[15] Liberzon, D. (2012). *Calculus of variations and optimal control theory* (Princeton University Press, Princeton, NJ)

[16] Bryson, A. E. (1975). *Applied optimal control: optimization, estimation and control* (CRC Press, New York). doi:https://doi.org/10.1201/9781315137667

[17] Khaneja, N., Reiss, T., Kehlet, C., Schulte-Herbrüggen, T., and Glaser, S. J. (2005). Optimal control of coupled spin dynamics: design of NMR pulse sequences by gradient ascent algorithms. *J. Magn. Res.* 172, 296–305. doi:10.1016/j.jmr.2004.11.004

[18] Reich, D. M., Ndong, M., and Koch, C. P. (2012). Monotonically convergent optimization in quantum control using Krotov’s method. *The Journal of Chemical Physics* 136, 104103

[19] Werschnik, J. and Gross, E. K. U. (2007). Quantum optimal control theory. *J. Phys. B* 40, R175

[20] Goerz, M. H., Carrasco, S. C., and Malinovsky, V. S. (2022). Quantum Optimal Control via Semi-Automatic Differentiation. *Quantum* 6, 871

[21] Dionis, E. and Sugny, D. (2023). Time-optimal control of two-level quantum systems by piecewise constant pulses. *Phys. Rev. A* 107, 032613

[22] Harutyunyan, M., Holweck, F., Sugny, D., and Guérin, S. (2023). Digital optimal robust control. *Phys. Rev. Lett.* 131, 200801

[23] Bloch, I., Dalibard, J., and Zwerger, W. (2008). Many-body physics with ultracold gases. *Rev. Mod. Phys.* 80, 885–964

[24] Eckardt, A. (2017). Colloquium: Atomic quantum gases in periodically driven optical lattices. *Rev. Mod. Phys.* 89, 011004

[25] Gross, C. and Bloch, I. (2017). Quantum simulations with ultracold atoms in optical lattices. *Science* 357, 995

[26] Dupont, N., Chatelain, G., Gabardos, L., Arnal, M., Billy, J., Peaudecerf, B., et al. (2021). Quantum state control of a bose-einstein condensate in an optical lattice. *PRX Quantum* 2, 040303
[27] Jäger, G. and Hohenester, U. (2013). Optimal quantum control of bose-einstein condensates in magnetic microtraps: Consideration of filter effects. *Phys. Rev. A* 88, 035601

[28] Jäger, G., Reich, D. M., Goerz, M. H., Koch, C. P., and Hohenester, U. (2014). Optimal quantum control of bose-einstein condensates in magnetic microtraps: Comparison of gradient-ascent-pulse-engineering and krotov optimization schemes. *Phys. Rev. A* 90, 033628

[29] Rodzinka, T., Dionis, E., Calmels, L., Beldjoudi, S., Béguin, A., Guéry-Odelin, D., et al. (2024). Optimal floquet engineering for large scale atom interferometers. *Nature Comm.* 15, 10281

[30] Saywell, J., Carey, M., Belal, M., Kuprov, I., and Freegarde, T. (2020). Optimal control of raman pulse sequences for atom interferometry. *Journal of Physics B: Atomic, Molecular and Optical Physics* 53, 085006

[31] Sørensen, J. J. W. H., Aranburu, M. O., Heinzel, T., and Sherson, J. F. (2018). Quantum optimal control in a chopped basis: Applications in control of bose-einstein condensates. *Phys. Rev. A* 98, 022119

[32] Hocker, D., Yan, J., and Rabitz, H. (2016). Optimal nonlinear coherent mode transitions in bose-einstein condensates utilizing spatiotemporal controls. *Phys. Rev. A* 93, 053612

[33] Mennemann, J.-F., Matthes, D., Weishäupl, R.-M., and Langen, T. (2015). Optimal control of bose-einstein condensates in three dimensions. *New Journal of Physics* 17, 113027

[34] Chen, X., Torrontegui, E., Stefanatos, D., Li, J.-S., and Muga, J. G. (2011). Optimal trajectories for efficient atomic transport without final excitation. *Phys. Rev. A* 84, 043415

[35] Zhang, Q., Muga, J. G., Guéry-Odelin, D., and Chen, X. (2016). Optimal shortcuts for atomic transport in anharmonic traps. *Journal of Physics B: Atomic, Molecular and Optical Physics* 49, 125503

[36] Bucker, R., Berrada, T., van Frank, S., Schaff, J.-F., Schumm, T., Schmiedmayer, J., et al. (2013). Vibrational state inversion of a bose-einstein condensate: optimal control and state tomography. *Journal of Physics B: Atomic, Molecular and Optical Physics* 46, 104012

[37] van Frank, S., Negretti, A., Berrada, T., Bucker, R., Montangero, S., Schaff, J.-F., et al. (2014). Interferometry with non-classical motional states of a bose-einstein condensate. *Nature Comm.* 5, 4009

[38] Pötting, S., Cramer, M., and Meystre, P. (2001). Momentum-state engineering and control in bose-einstein condensates. *Phys. Rev. A* 64, 063613

[39] Weidner, C. A., Yu, H., Kosloff, R., and Anderson, D. Z. (2017). Atom interferometry using a shaken optical lattice. *Phys. Rev. A* 95, 043624

[40] Weidner, C. A. and Anderson, D. Z. (2018). Simplified landscapes for optimization of shaken lattice interferometry. *New Journal of Physics* 20, 075007
[41] Bason, M., Viteau, M., Malossi, N., Huillery, P., Arimondo, E., Ciampini, D., et al. (2012). High-fidelity quantum driving. *Nat. Phys.* 8, 147

[42] Arrouas, F., Ombredane, N., Gabardos, L., Dionis, E., Dupont, N., Billy, J., et al. (2023). Floquet operator engineering for quantum state stroboscopic stabilization. *Comptes Rendus. Physique* 24, 173-185

[43] Amri, S., Corgier, R., Sugny, D., Rasel, E. M., Gaaloul, N., and Charron, E. (2019). Optimal control for the transport of bose einstein condensates with an atom chip. *Sci. Rep.* 9, 5346

[44] Adriazola, J. and Goodman, R. H. (2022). Reduction-based strategy for optimal control of bose-einstein condensates. *Phys. Rev. E* 105, 025311

[45] Hohenester, U. (2014). Octbec—a matlab toolbox for optimal quantum control of bose-einstein condensates. *Computer Physics Communications* 185, 194-216

[46] Dupont, N., Arrouas, F., Gabardos, L., Ombredane, N., Billy, J., Peaudecerf, B., et al. (2023). Phase-space distributions of bose-einstein condensates in an optical lattice: optimal shaping and reconstruction. *New Journal of Physics* 25, 013012

[47] Zhou, X., Jin, S., and Schmiedmayer, J. (2018). Shortcut loading a bose-einstein condensate into an optical lattice. *New J. Phys.* 20, 055005

[48] Zhang, Y., Lapert, M., Sugny, D., Braun, M., and Glaser, S. J. (2011). Time-optimal control of spin 1/2 particles in the presence of radiation damping and relaxation. *The Journal of Chemical Physics* 134, 054103

[49] Chen, X., Ban, Y., and Hegerfeldt, G. C. (2016). Time-optimal quantum control of nonlinear two-level systems. *Phys. Rev. A* 94, 023624

[50] Dorier, V., Gevorgyan, M., Ishkhanyan, A., Leroy, C., Jauslin, H. R., and Guérin, S. (2017). Nonlinear stimulated raman exact passage by resonance-locked inverse engineering. *Phys. Rev. Lett.* 119, 243902

[51] Zhu, J.-J., Chen, X., Jauslin, H.-R., and Guérin, S. (2020). Robust control of unstable nonlinear quantum systems. *Phys. Rev. A* 102, 052203

[52] Zhu, J.-J. and Guérin, S. (2024). Control of nonlinear quantum systems: From nonlinear rabi oscillations to robust two-stage optimal inverse engineering. *Phys. Rev. A* 110, 063503

[53] Feng, B. and Zhao, D. (2016). Optimal bilinear control of gross-pitaevskii equations with coulombian potentials. *Journal of Differential Equations* 260, 2973-2993

[54] Hintermüller, M., Marahrens, D., Markowich, P. A., and Sparber, C. (2013). Optimal bilinear control of gross-pitaevskii equations. *SIAM Journal on Control and Optimization* 51, 2509-2543

[55] Guérin, S. and Jauslin, H. (1999). Grid methods and hilbert space basis for simulations of quantum dynamics. *Computer Physics Communications* 121-122, 496-498. Proceedings of the Europhysics Conference on Computational Physics CCP 1998
[56] Light, J. C. (1992). *Discrete Variable Representations in Quantum Dynamics* (Boston, MA: Springer US). 185–199

[57] Littlejohn, R. G., Cargo, M., Carrington, J., Tucker, Mitchell, K. A., and Poirier, B. (2002). A general framework for discrete variable representation basis sets. *The Journal of Chemical Physics* 116, 8691–8703

[58] Kosloff, D. and Kosloff, R. (1983). A fourier method solution for the time dependent schrödinger equation as a tool in molecular dynamics. *Journal of Computational Physics* 52, 35–53

[59] Leforestier, C., Roncero, O., Bisseling, R., Cerjan, C., Feit, M., Friesner, R., et al. (1990). A comparison of different propagation schemes for the time-dependent schrödinger equation. *Journal of Computational Physics* 89, 490–491

[60] Dalfovo, F., Giorgini, S., Pitaevskii, L. P., and Stringari, S. (1999). Theory of bose-einstein condensation in trapped gases. *Rev. Mod. Phys.* 71, 463–512

[61] Leggett, A. J. (2001). Bose-einstein condensation in the alkali gases: Some fundamental concepts. *Rev. Mod. Phys.* 73, 307–356

[62] Morsch, O. and Oberthaler, M. (2006). Dynamics of bose-einstein condensates in optical lattices. *Rev. Mod. Phys.* 78, 179–215

[63] Pozzoli, E. (2024). Small-time global approximate controllability of bilinear wave equations. *Journal of Differential Equations* 388, 421–438

[64] Chambrion, T. and Pozzoli, E. (2023). Small-time bilinear control of schrödinger equations with application to rotating linear molecules. *Automatica* 153, 111028

[65] Deffner, S. (2022). Nonlinear speed-ups in ultracold quantum gases. *Europhysics Letters* 140, 48001

[66] Deffner, S. and Campbell, S. (2017). Quantum speed limits: from heisenberg's uncertainty principle to optimal quantum control. *Journal of Physics A: Mathematical and Theoretical* 50, 453001

[67] Kobzar, K., Ehni, S., Skinner, T. E., Glaser, S. J., and Luy, B. (2012). Exploring the limits of broadband 90 and 180 universal rotation pulses. *Journal of Magnetic Resonance* 225, 142–160

[68] Lapert, M., Ferrini, G., and Sugny, D. (2012). Optimal control of quantum superpositions in a bosonic Josephson junction. *Phys. Rev. A* 85, 023611

[69] Li, J.-S. and Khaneja, N. (2009). Ensemble control of bloch equations. *IEEE Transactions on Automatic Control* 54, 528–536. doi:10.1109/TAC.2009.2012983

[70] Van Damme, L., Ansel, Q., Glaser, S. J., and Sugny, D. (2017). Robust optimal control of two-level quantum systems. *Phys. Rev. A* 95, 063403
[71] Palao, J. P. and Kosloff, R. (2002). Quantum computing by an optimal control algorithm for unitary transformations. *Phys. Rev. Lett.* 89, 188301

[72] Palao, J. P. and Kosloff, R. (2003). Optimal control theory for unitary transformations. *Phys. Rev. A* 68, 062308

[73] Xu, S., Schmiedmayer, J., and Sanders, B. C. (2022). Nonlinear quantum gates for a bose-einstein condensate. *Phys. Rev. Res.* 4, 023071