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

1 Introduction

Quantum technologies seek to exploit the specific properties of quantum systems for real-world applications in computing, sensing, simulations or communications (Acín et al., 2018). 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 (Glaser et al., 2015; Brif et al., 2010; Koch et al., 2022). QOC is becoming a key tool in many different experimental platforms, ranging from superconducting circuits (Werninghaus et al., 2021; Abdelhafez et al., 2020; Wilhelm et al., 2020) to cold atoms (Hohenester et al., 2007; Ansel et al., 2024; van Frank et al., 2016), molecular physics (Koch et al., 2019) or NV centers (Rembold et al., 2020). Despite the maturity and effectiveness of optimal control techniques which rely on a rigorous mathematical framework, namely, the Pontryagin Maximum Principle [PMP, see (Ansel et al., 2024; Pontryagin et al., 1962; Boscain et al., 2021; Liberzon, 2012) for details], developments and adaptations of standard methods are necessary to take into account experimental limitations or additional degrees of freedom in the experiment (Bryson, 1975; Khaneja et al., 2005; Reich et al., 2012; Werschnik and Gross, 2007; Goerz et al., 2022; Harutyunyan et al., 2023; Dionis and Sugny, 2023). Bose-Einstein condensates (BEC) constitute a promising system (Eckardt, 2017; Bloch et al., 2008) for applications in quantum sensing and simulation (Gross and Bloch, 2017) in which optimal control may play a major role to prepare specific states or to improve the estimation of unknown parameters (Koch et al., 2022). This approach has been applied in a variety of works both theoretically and experimentally, with very good agreement (Hohenester et al., 2007; Ansel et al., 2024; van Frank et al., 2016; Dupont et al., 2021; Jäger and Hohenester, 2013; Jäger et al., 2014; Rodzinka et al., 2024; Saywell et al., 2020; Sørensen et al., 2018; Hocker et al., 2016; Mennemann et al., 2015; Chen et al., 2011; Zhang et al., 2016; Bucker et al., 2013; van Frank et al., 2014; Pötting et al., 2001; Weidner et al., 2017; Weidner and Anderson, 2018; Bason et al., 2012; Zhou et al., 2018; Arrouas et al., 2023; Dupont et al., 2023; Amri et al., 2019; Adriazola and Goodman, 2022; Hohenester, 2014). A specific 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 (Dupont et al., 2021). However, a majority of studies consider a simplified situation in which the non-linear term of the Gross-Pitaevski equation describing the dynamics of the BEC is neglected. This approximation 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 equation. It also greatly simplifies the implementation of QOC as recently shown by our group in Dupont et al. (2021); Dupont et al. (2023) where a gradient-based algorithm, GRAPE (Khaneja et al., 2005), 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 Hohenester et al. (2007), Jäger and Hohenester (2013), Jäger et al. (2014). The control problem has been solved analytically for a two-level quantum system in Zhang et al. (2011), Chen et al. (2016), Dorier et al. (2017), Zhu and Guérin (2024), Zhu et al. (2020). The control of the Gross-Pitaevski equation has also been the subject of a series of mathematical papers [see (Feng and Zhao, 2016; Hintermüller et al., 2013) to mention a few]. In this article, we propose to revisit such works by applying the PMP in the presence of mean-field interactions and deriving the corresponding algorithm from this optimization principle. Intensive use of the pseudo-spectral approach (Littlejohn et al., 2002; Light, 1992; Leforestier et al., 1990; Kosloff and Kosloff, 1983; Guérin and Jauslin, 1999) and FBR-DVR bases (Finite Basis Representation and Discrete Variable Representation) is necessary to analytically express the different quantities and accelerate the numerical calculation. 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 accessible to quantum simulation: higher dimensionalities e.g., can radically change the physics of localization (Morsch and Oberthaler, 2006), and also increase the role of interaction, giving access to many-body phenomena such as the Mott transition (Bloch et al., 2008). 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 Mennemann et al. (2015), and Zhou et al. (2018). In Mennemann et al. (2015), 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 Zhou et al. (2018). We present in this work another numerical implementation 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 geometries 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 Section 4. Conclusion and prospective views are given in Section 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 (Ansel et al., 2024; Dupont et al., 2021; Arrouas et al., 2023; Dupont et al., 2023).

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 × 1 0 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 λ = 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 Dupont et al. (2021), Arrouas et al. (2023), Dupont et al. (2023), we assume here that the lattice depth is fixed and only the relative phase of the lasers can be controlled.

The wave function | ψ ( t ) 〉 describing the state of the BEC is governed by the time-dependent Schrödinger equation which can be expressed aswhere p ̂ and 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 π / ( λ / 2 ) the wave vector and E L = ℏ 2 k L 2 / 2 m the characteristic energy of the lattice. The parameters s and φ 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 = 1 2 m ω ext x 2 with ω ext = 2 π × 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 ≪ 2 π / ω ext ).

ı ℏ d | ψ t 〉 d t = p ̂ 2 2 m − s E L 2 cos k L x ̂ + φ t | ψ t 〉 ,

In order to work with dimensionless coordinates, we introduce the following change of variables t → E L ℏ t and x → k L x. We obtainwhere p ̂ = − ı ∂ ∂ x and 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 | ϕ α 〉, of wave function in the x representation ϕ α ( x ) = 1 2 π e ı α x and of eigenvalue α. Since the potential is periodic in x, the Bloch theorem states that α = n + q, where n ∈ Z is a relative integer and q ∈ [ − 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

ı d | ψ t 〉 d t = p ̂ 2 − s 2 cos x ̂ + φ t | ψ t 〉 , (1)

It is then straightforward to show that the coefficients c q , n are solutions of the following equation

ı c ̇ q , n = n + q 2 c q , n − s 4 e ı φ t c q , n − 1 + e − ı φ t c q , n + 1 .

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

ı d | ψ t 〉 d t = H ̂ | ψ t 〉 = H ̂ 0 + cos φ t H ̂ 1 + sin φ t H ̂ 2 | ψ t 〉 , (2)

with

H ̂ 0 = ⋱ … 0 n − 1 + q 2 0 0 0 … … 0 0 n + q 2 0 0 … … 0 0 0 n + 1 + q 2 0 … ⋱ ,

and

H ̂ 1 = ⋱ ⋱ … − s 4 0 − s 4 0 0 … … 0 − s 4 0 − s 4 0 … … 0 0 − s 4 0 − s 4 … ⋱ ⋱ , H ̂ 2 = ⋱ ⋱ … − ı s 4 0 ı s 4 0 0 … … 0 − ı s 4 0 ı s 4 0 … … 0 0 − ı s 4 0 ı s 4 … ⋱ ⋱ .

From a numerical point of view, the infinite-dimensional Hilbert space is truncated to a finite one such that | n | ≤ 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 × 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 | ψ t 〉, with a very good match between theory and experiment (Dupont et al., 2021; Dupont et al., 2023). To this aim, the optimal phase is calculated from the application of GRAPE (Ansel et al., 2024; Khaneja et al., 2005) to the dynamical system Equation 2. The goal of the optimal control problem is to maximize the fidelity F = | 〈 ψ t | ψ ( t f ) 〉 | 2 at time t f. The PMP leads to the Pontryagin Hamiltonian H p = I ( 〈 χ ( t ) | H ̂ | ψ ( t ) 〉 ), with | χ ( t ) 〉 the adjoint state, whose dynamics are also governed by Equation 2, with the final condition χ 0 〈 ψ t | ψ ( t f ) 〉 | ψ t 〉, where χ 0 is the dual variable of the cost. χ 0 is set to 1/2 in the numerical simulation. We refer the interested reader to the recent tutorial (Ansel et al., 2024) for a comprehensive introduction to the PMP for quantum optimal control.

Using the maximization condition of the PMP, the control is iteratively improved aswhere ϵ is a small positive parameter. The states | ψ ( t ) 〉 and | χ ( t ) 〉 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 | ϕ 0 + 0 〉. The target states are described in the basis | ϕ 0 + n 〉 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 Dupont et al. (2023) for the mathematical definition of these states]. We denote by | g ( x c , p , ξ ) 〉 a squeezed state centered in ( x c , p c ) with the squeezing parameter ξ. 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 μ s in real units). More precisely, the target states are then the momentum state | ϕ 0 + 2 〉, the superposition 1 / 3 ( | ϕ 0 − 2 〉 + | ϕ 0 + 0 〉 + | ϕ 0 + 2 〉 ), and the centered squeezed state | g ( x c = 0 , p c = 0 , ξ = 1 / 3 ) 〉.

φ t ↦ φ t + ϵ I 〈 χ t | − sin φ t H ̂ 1 + cos φ t H ̂ 2 | ψ t 〉 ,

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 Equation 1, 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 Equation 1, proportional to the square modulus of the wave function | ψ ( x , t ) | 2 (Dalfovo et al., 1999; Leggett, 2001). 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 ( | ϕ q + n 〉 ) n ∈ 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 1) 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 | ψ | 2 and to the Gross–Pitaevskii equation,

ı ∂ ψ x , t ∂ t = − ∂ 2 ∂ x 2 − s 2 cos x + φ t + β | ψ x , t | 2 ψ x , t . (3)

The weight of the non-linear term is given by the dimensionless parameter β. In the experiments, this parameter is of the order of 0.5 ± 0.2 Ansel et al. (2024). This value depends on the number of atoms and on the frequency of the hybrid trap.

Expanding the state | ψ 〉 over the eigenvectors of the momentum operator | ψ 〉 = ∑ n ∈ Z c q , n | ϕ q + n 〉, the non-linear term | ψ | 2 can be expressed in position representation as,

ψ 2 = ψ x , t * ψ x , t = 1 2 π ∑ m ∈ Z c q , m * e − ı q + m x 1 2 π ∑ l ∈ Z c q , l e ı q + l x

which leads to

ψ 2 = 1 2 π ∑ m , l ∈ Z c q , m * c q , l e ı q + l − m x .

The dynamics of the coefficients c q , n are then given by

ı c ̇ q , n = n + q 2 c q , n − s 4 e ı φ t c q , n − 1 + e − ı φ t c q , n + 1 + β 2 π ∑ m , l ∈ Z c q , m * c q , l c q , n + m − l . (4)

Unlike the Schrödinger Equation 2, 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 Equation 4 takes about 15 s, for a constant control φ ( t ) = 0 (see Figure 2). Numerical simulations have been performed on a standard laptop. Optimal control, which requires a hundred propagations, is therefore difficult to apply with this approach, as the design of a single control would require a computational time of the order of 1 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 Section 3.2.

3.2 The FBR-DVR bases

An efficient approach to express the non-linear term of Equation 3 in matrix form is the FBR-DVR (or pseudo-spectral) approach (Littlejohn et al., 2002; Light, 1992; Leforestier et al., 1990; Kosloff and Kosloff, 1983; Guérin and Jauslin, 1999). 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 | ϕ q + n 〉 of the momentum operator. In this basis, it is straightforward to express the operators p ̂ and cos ( 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 | ψ ( 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 Equation 3. The FBR basis corresponds to the basis of the eigenvectors | ϕ q + n 〉 of the momentum operator. From a numerical point of view, we work in a finite dimensional sub-space H N of H of dimension N = 2 n max + 1, with n max ∈ N and − n max ≤ n ≤ n max. The potential energy is a periodic spatial function, of period 2 π. We choose the interval [ 0,2 π ] for the position x. The DVR basis is built from the discretization of the variable x ∈ 0,2 π, such that x j = 2 π N j with 0 ≤ j < N. A scalar product is given by,

∀ Ψ , Φ ∈ H N , ∫ 0 2 π Ψ * x Φ x d x = 2 π N ∑ j = 0 N − 1 Ψ * x j Φ x j .

From the orthogonality of functions ϕ q + n we deduce that we have approximately:

2 π N ∑ j = 0 N − 1 ϕ q + m * x j ϕ q + n x j = δ m , n . (5)

From the definition of the plane wave functions, we have also for n max large enough:

2 π N ∑ n = − n max n max ϕ q + n x i ϕ q + n * x j = δ i , j . (6)

We introduce the matrix R ̂ defined as

R ̂ = R 0 , − n max R 0 , − n max + 1 … ⋮ ⋱ R N − 1 , − n max R N − 1 , n max , (7)

with matrix elements R j , n given by

R j , n = 1 N e ı 2 π N q + n j .

Equations 5, 6 are equivalent to R ̂ † R ̂ = I ̂ and R ̂ R ̂ † = I ̂. The matrix R ̂ is thus a unitary matrix. The DVR basis | u j 〉 j ∈ [ 0 , N − 1 ] is then built from the following kets as

This basis connects a basis of H N to the spatial discretization basis. We deduce that an operator W ̂ depending on x has a diagonal representation in the DVR basis given by

〈 u i | W ̂ | u j 〉 ≃ W x j δ i , j .

The matrix R ̂ allows to express a known operator in the DVR basis to the FBR basis, and vice versa

W ̂ FBR = R ̂ † W ̂ DVR R ̂ .

Since the square modulus of the wave function, | ψ ( x ) | 2, depends on the position x as shown in Equation 4, 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 G ̂ DVR the matrix of | ψ ( x ) | 2 in the DVR basis,

G ̂ DVR = G 0 G 1 ⋱ G N − 1 , (8)

where

G j = 1 2 π ∑ n = − n max n max c q , n e ı 2 π N q + n j 2 .

The parameters c q , n are the coefficients of the state in the FBR basis. Using the notations introduced in Section 2 and the relations Equations 7, 8, Equation 3 can be written in matrix form in the FBR basis as follows

ı d | ψ 〉 d t = H ̂ G P | ψ 〉 = H ̂ 0 + cos φ t H ̂ 1 + sin φ t H ̂ 2 + β R ̂ † G ̂ DV R R ̂ | ψ 〉 , (9)

with H ̂ GP the Gross–Pitaevskii Hamiltonian.

For comparison, the propagation for a constant control is here of the order of 0.6 s, i.e., a total of around 2 min for the calculation time of the optimal control. Figure 2 shows the time evolution of the projection of the state | ψ 〉 on | ϕ 0 + 0 〉, calculated from the matrix exponential approach and the RK4 method. The initial state is the squeezed state | g ( 0,0,1 / 2 ) 〉 and the phase φ ( t ) = 0. For β = 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.

We consider the state-to-state transfer defined with the fidelity F = | 〈 ψ t | ψ ( t f ) 〉 | 2, with the dynamics of Equation 3. In this case, the Pontryagin Hamiltonian H p can be written as

H p = I 〈 χ | H ̂ G P | ψ 〉 ,

with H ̂ G P, the Gross–Pitaevskii Hamiltonian in position representation,

H ̂ G P = − ∂ 2 ∂ x 2 − s 2 cos x + φ t + β | ψ x , t | 2 .

The time evolution of the adjoint state is given by the relation, χ ̇ ( x , t ) = − ∂ ψ H p, which leads to

ı ∂ χ ∂ t = − ∂ 2 ∂ x 2 − s 2 cos x + φ t + 2 β ψ 2 χ − β ψ 2 χ * . (10)

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 χ ( x , t f ) = χ 0 ∂ ψ F yields

χ x , t f = 2 χ 0 〈 ψ t | ψ t f 〉 ψ t ,

with χ 0 = 1 / 2 in the numerical simulations. Finally, the maximization condition of the Pontryagin Hamiltonian leads to the iterative procedure of GRAPEwhere ϵ is a small positive parameter. The state and the adjoint state are respectively propagated forward and backward in time. The differential Equation 10 associated with the adjoint state can be written as

φ t ↦ φ t + ϵ I 〈 χ | ∂ H ̂ GP ∂ φ | ψ 〉

ı ∂ χ ∂ t = − ∂ 2 ∂ x 2 − s 2 cos x + φ t + β ψ 2 χ − 2 ı β I χ * ψ ψ .

The backward propagation can be done in matrix form with the final conditions | ψ ( t f ) 〉 = | ψ f 〉 and | χ ( t f ) 〉 = 〈 ψ t | ψ ( τ f ) 〉 | ψ t 〉, with the following extended systemwhere H ̂ GP is the Gross–Pitaevskii Hamiltonian in matrix form given by Equation 9, and I ̂ DVR is the diagonal matrix associated with the term I χ * ψ in the DVR basis. The matrix elements of I ̂ DVR are I j = I 1 2 π ∑ m , n d q , m * c q , n e ı 2 π N ( q + n − m ) j, d q , n being the coefficients of the adjoint state | χ 〉 in the FBR basis, − n max ≤ n , m ≤ n max.

ı | ψ ̇ 〉 | χ ̇ 〉 = H ̂ GP 0 − 2 ı β R ̂ † I ̂ DVR R ̂ H ̂ GP | ψ 〉 | χ 〉 ,

As an illustrative example, we consider the transfer from the state | ψ ( 0 ) 〉 = | ϕ 0 + 0 〉 to the squeezed state | g 0,0,3 / 2 〉. The control time is set to t f = 150 μ s, the lattice depth to s = 5, and the quasi-momentum to q = 0. A first control is optimized for β = 0 as represented in Figure 3. The dynamics are then propagated for this control with a non-linear coefficient β going from 0 to 1, β ∈ [ 0,1 ]. Figure 4 shows the evolution of the fidelity with respect to the nonlinear parameter β. The grey area indicates the experimental uncertainty on this parameter. We observe that for β = 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 system. A second control is optimized by using the non-linear version of GRAPE for β = 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 β 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 a i the vectors of the basis of the direct lattice and N the corresponding dimension. A node is an element of this set,

D = ∑ i = 1 N j i a i , j i ∈ Z . (11)

The reciprocal lattice is described by the set,where the vectors b i are defined from the relation,

R = ∑ i = 1 N j i b i , j i ∈ Z (12)

a i ⋅ b j = 2 π δ i j ,

with a i ∈ D and b j ∈ R. Using the Bloch’s theorem, it can be shown that the wave function ψ ( r ) solution of the time-independent Schrödinger equation,where the potential V is periodic, V ( r + R ) = V ( r ) for any R ∈ D and r = ∑ i = 1 N r i a i is the position of the system, can be expressed aswith q = ∑ i = 1 N q i b i the two-dimensional quasi-momentum and u q ( r ) a periodic function in D.

− ℏ 2 2 m ∂ 2 ∂ r 2 + V r ψ r = E ψ r , (13)

ψ q r = e ı q ⋅ r u q r , (14)

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 (Morsch and Oberthaler, 2006). As an illustrative example, we consider the following configuration, three lasers with the same angular frequency ω which propagate in the plane e x , e y of the laboratory frame with a linear polarization along the z- axis. The wave vectors are defined in this plane bywhere k = ω c is the wave number. The angle between the vectors is 2 π / 3. The total electric field is then given bywhere E 0 is the amplitude of the electric field, r the position vector and φ 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 H ̂ I = − μ ⋅ E, where μ is the induced dipole moment. At first order, the dipole moment is given by μ = − α 2 E, where α is the polarizability of the atom. This leads to H ̂ I = − α 2 E ⋅ E. The electric field is not resonant with the transition frequencies of the atom, and the period associated with the ω 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 potentialwhere φ j k = φ j − φ k, and 〈 ⋅ 〉 denotes the time average over one period τ, such that 2 π ω ≫ τ ≫ T dyn. Finally, we deduce that this average potential governs the 2D-dynamics (31), with r = x , y.

k 1 = k 1,0 k 2 = k 2 − 1 , 3 k 3 = − k 2 1 , 3 ,

E = E 0 ∑ i cos k i ⋅ r − ω t + φ i e z ,

V = 〈 H ̂ I 〉 = − 3 4 α E 0 2 − α 2 E 0 2 cos k 1 − k 2 ⋅ r + φ 12 + cos k 2 − k 3 ⋅ r + φ 23 + cos k 3 − k 1 ⋅ r + φ 31 ,

The direct lattice nodes correspond to the minima of the potential energy for φ j k = 0. We have

k 1 − k 2 ⋅ r = 0 m o d 2 π k 2 − k 3 ⋅ r = 0 m o d 2 π k 3 − k 1 ⋅ r = 0 m o d 2 π

A node of the direct lattice belongs to the set D Equation 11 with,

a 1 = − 4 π 3 k 2 k 3 , a 2 = 4 π 3 k 2 k 2 .

A node of the reciprocal lattice is an element of the set R Equation 12 with

b 1 = k 1 − k 3 , b 2 = k 2 − k 1 .

Since the potential is periodic, Bloch’s theorem applies, and eigenvectors are as in Equation 14. Within a given sub-Hilbert space of quasi-momentum q, a generic wave function can be written as:

ψ r = ∑ m , n ∈ Z c q , m , n e ı q ⋅ r e ı m b 1 + n b 2 ⋅ r . (15)

We denote by ϕ m , n = e ı m b 1 + n b 2 ⋅ r the functions of the Fourier basis. The coefficients c 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 (Dupont et al., 2021). Plugging the expression of ψ ( r , t ) of Equation 15 into Equation 13, we obtain,where the constant term − 3 4 α E 0 2 of the potential V is removed, and V 0 = α E 0 2 is the lattice depth. From the relations, k 1 − k 2 = − b 2, k 2 − k 3 = b 1 + b 2, k 3 − k 1 = − b 1, b 1 ⋅ b 1 = b 2 ⋅ b 2 = 3 k 2 and b 1 ⋅ b 2 = − 3 2 k 2, we arrive at

ı ℏ ∑ m , n ∈ Z c ̇ m , n e ı m b 1 + n b 2 r = ∑ m , n ∈ Z ℏ 2 2 m c m , n m b 1 + n b 2 2 e ı m b 1 + n b 2 r − ∑ m , n ∈ Z V 0 4 c m , n e ı k 1 − k 2 ⋅ r + φ 12 e ı m b 1 + n b 2 r + e − ı k 1 − k 2 ⋅ r + φ 12 e ı m b 1 + n b 2 r − ∑ m , n ∈ Z V 0 4 c m , n e ı k 2 − k 3 ⋅ r + φ 23 e ı m b 1 + n b 2 r + e − ı k 2 − k 3 ⋅ r + φ 23 e ı m b 1 + n b 2 r − ∑ m , n ∈ Z V 0 4 c m , n e ı k 3 − k 1 ⋅ r + φ 31 e ı m b 1 + n b 2 r + e − ı k 3 − k 1 ⋅ r + φ 31 e ı m b 1 + n b 2 r ,

ı ℏ ∑ m , n ∈ Z c ̇ m , n e ı m b 1 + n b 2 r = ∑ m , n ∈ Z 3 k 2 ℏ 2 2 m c m , n m 2 + n 2 − m n e ı m b 1 + n b 2 r − ∑ m , n ∈ Z V 0 4 c m , n e ı φ 12 e ı m b 1 + n − 1 b 2 r + e − ı φ 12 e ı m b 1 + n + 1 b 2 r − ∑ m , n ∈ Z V 0 4 c m , n e ı φ 23 e ı m + 1 b 1 + n + 1 b 2 r + e − ı φ 23 e ı m − 1 b 1 + n − 1 b 2 r − ∑ m , n ∈ Z V 0 4 c m , n e ı φ 31 e ı m − 1 b 1 + n b 2 r + e − ı φ 31 e ı m + 1 b 1 + n b 2 r .

Projecting this equation onto the state ϕ m , n and introducing the dimensionless time τ = 3 ℏ k 2 2 m t and lattice depth s = 2 m V 0 / ( 3 ℏ 2 k 2 ), we obtain the differential equation governing the dynamics of the coefficients c m , n,

ı c ̇ m , n = m 2 + n 2 − m n c m , n − s 4 e ı φ 12 c m , n + 1 + e − ı φ 12 c m , n − 1 + e ı φ 23 c m − 1 , n − 1 + e − ı φ 23 c m + 1 , n + 1 + e ı φ 31 c m + 1 , n + e − ı φ 31 c m − 1 , n . (16)

4.2 Optimal control problem

At this point, Equation 16 can be put in matrix from. We have

ı d | ψ t 〉 d t = H ̂ 0 + e ı φ 12 H ̂ 12 + + e − ı φ 12 H ̂ 12 − + e ı φ 23 H ̂ 23 + + e − ı φ 23 H ̂ 23 − + e ı φ 31 H ̂ 31 + + e − ı φ 31 H ̂ 31 − | ψ t 〉 ,

with

H ̂ 0 = ∑ m , n ∈ Z m 2 + n 2 − m n | ϕ m , n 〉 〈 ϕ m , n | H ̂ 12 + = − s 4 ∑ m , n ∈ Z | ϕ m , n 〉 〈 ϕ m , n + 1 | H ̂ 12 − = − s 4 ∑ m , n ∈ Z | ϕ m , n 〉 〈 ϕ m , n − 1 | H ̂ 23 + = − s 4 ∑ m , n ∈ Z | ϕ m , n 〉 〈 ϕ m − 1 , n − 1 | H ̂ 23 − = − s 4 ∑ m , n ∈ Z | ϕ m , n 〉 〈 ϕ m + 1 , n + 1 | H ̂ 31 + = − s 4 ∑ m , n ∈ Z | ϕ m , n 〉 〈 ϕ m + 1 , n | H ̂ 31 − = − s 4 ∑ m , n ∈ Z | ϕ m , n 〉 〈 ϕ m − 1 , n |

From a numerical point of view, we use a Hilbert space of finite dimension such that − M ≤ m ≤ M and − N ≤ n ≤ N. The dimension of this space is d H = ( 2 M + 1 ) × ( 2 N + 1 ). In the numerical simulations, M and N have been set to 5. The GRAPE algorithm can then be applied to this case (Ansel et al., 2024). For each control φ 12, φ 23 and φ 31, the maximization condition given by the PMP leads to the iterative procedure of GRAPEwhere | χ 〉 is the adjoint state and ℓ corresponds to the index 12, 23 or 31.

φ ℓ ↦ φ ℓ + ϵ I 〈 χ t | ı H ̂ ℓ + − H ̂ ℓ − | ψ t 〉 ,

As an illustrative example, we consider the state-to-state transfer from the initial state | ψ ( 0 ) 〉 = | 0,0 〉 to a target state | ψ t 〉 = 1 2 | 3,3 〉 + | − 3 , − 3 〉 with a control time t f = 250 μ s. The values of the other parameters are the same as in Section 2. Figure 5 represents the optimization result for this control problem. Note that the controls are not independent since φ 12 and φ 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 φ 2,3 and φ 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 Figure 7, we consider a non-symmetric state in m and n with | φ t 〉 = 1 2 ( | ϕ 1,2 〉 + | ϕ − 3 , − 1 〉 ). 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 (Chambrion and Pozzoli, 2023; Pozzoli, 2024). However, such results do not use a directly implementable control strategy and consider pulses with a very large amplitude. Another interesting question is to evaluate the role of the nonlinearity on the optimal solution. It has been shown in Deffner (2022) that the quantum speed limit Deffner and Campbell (2017) 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 (Kobzar et al., 2012; Li and Khaneja, 2009; Van Damme et al., 2017; Lapert et al., 2012). 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 (Palao and Kosloff, 2003; Palao and Kosloff, 2002). Recently, a rigorous foundation of quantum computing in the nonlinear case has been developed (Xu et al., 2022). 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.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.