# Efficient Manipulation of Bose-Einstein Condensates in a Double-Well
Potential

∗J. Adriazola(a) ,  R. H. Goodman(b) Department of Mathematical Sciences, New
Jersey Institute of Technology, University Heights, Newark, NJ, USA, 07102
and  P. G. Kevrekidis(c) Department of Mathematics and Statistics, University
of Massachusetts Amherst, Amherst, MA, USA, 01003-4515

###### Abstract.

We pose the problem of transferring a Bose-Einstein Condensate (BEC) from one
side of a double-well potential to the other as an optimal control problem for
determining the time-dependent form of the potential. We derive a reduced
dynamical system using a Galerkin truncation onto a finite set of
eigenfunctions and find that including three modes suffices to effectively
control the full dynamics, described by the Gross-Pitaevskii model of BEC. The
functional form of the control is reduced to finite dimensions by using
another Galerkin-type method called the chopped random basis (CRAB) method,
which is then optimized by a genetic algorithm called differential evolution
(DE). Finally, we discuss the extent to which the reduction-based optimal
control strategy can be refined by means of including more modes in the
Galerkin reduction.

Keywords. Quantum Control, Hamiltonian Dynamical Systems, Galerkin Methods

∗Corresponding author (Currently at the Universtiy of California, Santa
Barbara).

##  1\. Introduction

The dynamics of solitary waves in dispersive media with external potentials is
a topic of widespread scientific interest, as it arises in many areas of
application. For instance, in Bose-Einstein Condensates (BEC) of, e.g., alkali
gases, external potentials may be created using a variety of physical
mechanisms including optical and magnetic fields, and may consist of one or a
few wells or a periodic array, and may effectively confine the BEC to one,
two, or three space dimensions [23, 35, 37, 38]. Another appealing
experimental setting is the nonlinear propagation of light through photonic
crystals, and in the quasi-discrete realm of optical waveguides [20, 28].
Here, a spatially-dependent index of refraction induces an effective potential
[29].

In both these applications, the simplest potential enabling bifurcation
phenomena and nontrivial dynamics is arguably the double well potential. It
has been studied intensely in the atomic realm, following the hallmark
theoretical work of [39]. These predicted Josephson oscillations between the
wells and quantum self-trapping were subsequently realized experimentally in
[2], as well as a dynamical symmetry-breaking bifurcation later observed in
[50]. More recent experiments have added damping and driving, which may
present novel phenomena including stochastic resonance [43]. Relevant double-
well experiments have been conducted in the optical setting as well. The work
of [7] considered the double-well potential in the context of twin-core self-
guided laser beams in Kerr media, while [24] probed two-well dynamics using
photorefractive crystals. In this latter setting, additional phenomena were
demonstrated in a _triple-well_ potential [21].

Naturally, this large volume of experimental developments and control has
motivated a wide range of theoretical explorations in numerous further
directions. The relevant list is too long to do it proper justice, but we
mention some related studies. Some more mathematical examples include the
analysis of the double-well bifurcation structure [3, 18], the low-dimensional
representation of the associated dynamical problem (and its fidelity) [27,
30], and the effect of changing the nonlinear exponent on the bifurcation [40,
26]. Among the many more physical examples are the interactions of multiple
dispersive (e.g., atomic) species [48, 13, 46], incorporating beyond-mean-
field (i.e., many-body) effects [31, 36], and the effect of larger spatial
dimensions (and possibly four wells) [49], among others.

In this work, we aim to apply the deep understanding of the existence,
stability and nonlinear dynamics to the context of optimal control of BEC [25,
5]. The latter methodology has long been recognized as a versatile tool for
engineering on-demand quantum states of interest. In the early stages, the
framework of magnetic microtraps controllable by external parameters including
radio-frequency fields or/and wire currents was used to enable the preparation
of desired states [15] (see also the detailed analysis of the relevant
methodologies and their numerical Matlab-based implementation in [19, 14]).
Subsequently, such ideas have been applied to fully three-dimensional
settings, e.g., in the work of [32] and have been used recently by a subset of
the present authors in order to re-orient the density distribution of an
atomic BEC and alter the topology of its support [1].

Here, more concretely, we intend to show how the low-dimensional
representation available in the context of double (and more generally few [10,
49]) well potentials can be used as a basis for performing optimal control
analysis and for achieving desired end quantum states, both in the low-
dimensional setting, but also in the full mean-field model of the Gross-
Pitaevskii (nonlinear Schrödinger-type) partial differential equation (PDE)
[37, 38]. Along the way, we learn the following lesson, which we find
interesting and important. The two-mode expansion is prevalent in the study of
the double-well system, and widely acknowledged to describe the dynamics, both
qualitatively and even quantitatively. By contrast, we find that in the
context of optimal control, we must include (at least) the third mode in the
expansion to achieve useful agreement. We believe that such lessons may prove
useful for other practitioners in related contexts.

Our presentation is structured as follows. In Section 2, we present the
physical and mathematical setup and its reduced (two- and three-mode)
representation. Subsequently, in Section 3, we present the proposed optimal
control strategy. In Section 4, we display numerical results. Finally, Section
5 provides a summary of our conclusions, as well as a number of directions for
future study.

##  2\. Derivation of Reduced Model Systems

The approach to optimization we propose here is to apply optimal control to a
finite-dimensional model system, whose derivation we outline in this section.
In particular, we use Galerkin truncation to derive a low-dimensional
Hamiltonian system whose dynamics capture the essence of the full dynamics.
The latter, in turn, is described by a Gross-Pitaevskii equation (GPE) in one
spatial dimension. This system may be derived from a three-dimensional model
in the presence of an anisotropic potential that squeezes the condensate into
an effectively one-dimensional arrangement; see details and
nondimensionalization, e.g., in [23]. The initial and (final) desired
conditions used throughout are set to the stable fixed points of the finite-
dimensional model Hamiltonian we have derived. These fixed points correspond
to the two asymmetric states which exist in the presence of the barrier.
Indeed, our aim is to drive the atomic mass from a state predominantly
localized in one well to a state localized in the other. The efficiency will
be determined on the basis of how successful such a transfer is with an
appropriate definition of “fidelity” to the intended end state; more detail
will follow the analysis.

The GPE model is a nonlinear Schrödinger equation with a spatial potential

(2.1)
$i\partial_{t}\psi=\mathcal{L}\left(x,\mathbf{w}(t)\right)\psi+\mathcal{N}(\psi)=-\frac{1}{2}\partial_{x}^{2}\psi+V\left(x,\mathbf{w}(t)\right)\psi+\absolutevalue*{\psi}^{2}\psi,$

where $\mathbf{w}(t)$ is a time-dependent vector of $C^{0}$ control functions.
In addition to conserving a Hamiltonian energy, this system conserves the mass

(2.2) $M=\int_{-\infty}^{\infty}\absolutevalue*{\psi}^{2}\differential x,$

which in the BEC context is interpreted as the total number of atoms in the
condensate. The potential $V$ is chosen as the superposition of a quadratic
confining potential, usually implemented via magnetic fields [37, 38], and a
thin, tall barrier at the center, typically induced by an optical beam [50].
Added together, these form a prototypical double-well potential; see, also,
[44]. In particular, the potential takes the form

(2.3) $V=\frac{1}{2}u(t)x^{2}+v(t)\delta(x),$

where the first term models the magnetic confinement and the second term
models the localized repulsive barrier at the center, as we $v(t)>0$. The
time-dependent parameter vector is thus given by
$\mathbf{w}(t)=(u(t),v(t))^{T}.$

Let $\Phi_{\mathbf{w}}=\left\\{\varphi_{n}(x;\mathbf{w})\in
L^{2}({\mathbb{R}})\ |\ n=0,1,\ldots\right\\}$ be the set of normalized
eigenfunctions of the linear Schrödinger eigenvalue problem

(2.4) $\mathcal{L}(x,\mathbf{w})\varphi_{n}=E_{n}\varphi_{n},$

for a fixed parameter vector $\mathbf{w}$. Because $\Phi_{\mathbf{w}}$ is
complete in $L^{2}({\mathbb{R}})$, we may represent the solution to Eq. (2.1)
at time $t$ by the infinite series

$\psi(x,t)=\sum_{n=0}^{\infty}c_{n}(t)\varphi_{n}(x;\mathbf{w}(t)).$

Plugging this representation into Eq. (2.1), and projecting both sides of the
equation onto $\varphi_{n}$—in the $L^{2}({\mathbb{C}})$ sense—yields an
evolution equation for $c_{n}(t)$. Together, the evolution of the infinite
vector $\mathbf{c}=(c_{0}(t),c_{1}(t),\ldots)$ of complex amplitudes is then
equivalent to the evolution of $\psi$ under GPE.

To derive an approximate reduced system, we consider the truncated series
superposition of instantaneous eigenfunctions

(2.5) $\psi_{N+1}^{\rm
Gal}:=\sum_{n=0}^{N}c_{n}(t)\varphi_{n}(x;\mathbf{w}(t))$

for some fixed value of $N<\infty$. Ignoring any truncation error due to the
terms in the omitted tail of the series yields the system of interest. Such a
truncation has been rigorously justified in certain very simple cases, e.g.
[30, 11], but the method is commonly applied without rigorous justification.

###  2.1. Two-Mode Expansion

By setting $N=1$ in Expansion (2.5), we find the following two-mode
Hamiltonian system [44]:

(2.6)
$\begin{split}i\dot{c}_{0}&=\frac{\partial\mathcal{H}}{\partial\bar{c}_{0}}=\alpha
c_{0}+\gamma_{0}\absolutevalue{c_{0}}^{2}c_{0}+\gamma_{2}\left(c_{1}^{2}\bar{c}_{0}+2c_{0}\absolutevalue{c_{1}}^{2}\right),\\\
i\dot{c}_{1}&=\frac{\partial\mathcal{H}}{\partial\bar{c}_{1}}=\beta
c_{1}+\gamma_{1}\absolutevalue{c_{1}}^{2}c_{1}+\gamma_{2}\left(c_{0}^{2}\bar{c}_{1}+2c_{1}\absolutevalue{c_{0}}^{2}\right),\end{split}$

with instantaneous projection coefficients given by

(2.7)
$\displaystyle\alpha=\left\langle\mathcal{L}\varphi_{0},\varphi_{0}\right\rangle,\
\beta=\left\langle\mathcal{L}\varphi_{1},\varphi_{1}\right\rangle,\
\gamma_{0}=\norm*{\varphi_{0}}^{4},\ \gamma_{1}=\norm*{\varphi_{1}}^{4},\
\gamma_{2}=\left\langle\varphi_{0}^{2},\varphi_{1}^{2}\right\rangle.\ $

Its Hamiltonian reads:

(2.8)
$\mathcal{H}=\alpha\absolutevalue{c_{0}}^{2}+\beta\absolutevalue{c_{1}}^{2}+\frac{\gamma_{0}}{2}\absolutevalue{c_{0}}^{4}+\frac{\gamma_{1}}{2}\absolutevalue{c_{1}}^{4}+\gamma_{2}\left(\real\left\\{c_{0}^{2}\bar{c}_{1}^{2}\right\\}+2\absolutevalue{c_{0}}^{2}\absolutevalue{c_{1}}^{2}\right).$

This expansion holds under the general assumption that
$V(x,\mathbf{w}(t))=V(-x,\mathbf{w}(t))$, i.e., that the potential is even.
This system conserves a discrete form of the mass defined in Eq. (2.2),

(2.9) $M_{\rm
d}(t)=\absolutevalue*{c_{0}(t)}^{2}+\absolutevalue*{c_{1}(t)}^{2}.$

This system has stationary solutions of the form
$(c_{0}(t),c_{1}(t))=(\rho_{0},\rho_{1})e^{-i\Omega t}$. In particular, it has
a solution corresponding to the nonlinear continuation of the ground state
with $\rho_{1}=0$ and $\Omega=\alpha+\gamma_{0}\rho_{0}^{2}$, and a second
solution corresponding to the nonlinear continuation of the excited state,
with $\rho_{0}=0$ and $\Omega=\beta+\gamma_{1}\rho_{1}^{2}$. In the absence of
a barrier, i.e., for $v=0$, and total mass $M_{\rm d}=1$, these are the only
such states, and both are linearly stable. In the presence of a barrier,
however, the excited state can become unstable in a symmetry-breaking
pitchfork bifurcation.

We may take advantage of the conservation law of Eq. (2.9) to reduce the
system from two degrees of freedom to one as follows, which provides a
convenient visualization of the dynamics and bifurcation. Consider the
canonical transformation

(2.10) $c_{0}=Ae^{i\theta},\quad c_{1}=(q+ip)e^{i\theta}.$

We reduce the number of degrees of freedom from two to one using the conserved
mass (2.9) which now reads $A^{2}=M_{\rm d}-q^{2}-p^{2}$. The Hamiltonian in
these coordinates is given by

$\displaystyle\mathcal{H}=\frac{\gamma_{0}M_{\rm d}^{2}}{2}+\alpha M_{\rm
d}+q^{2}\left(\beta-\alpha+M_{\rm
d}\left(3\gamma_{2}-\gamma_{0}\right)\right)+p^{2}\left(\beta-\alpha+M_{\rm
d}\left(\gamma_{2}-\gamma_{0}\right)\right)$ (2.11)
$\displaystyle+\left(\frac{\gamma_{0}}{2}+\frac{\gamma_{1}}{2}-\gamma_{2}\right)p^{4}+\left(\gamma_{0}+\gamma_{1}-4\gamma_{2}\right)p^{2}q^{2}+\left(\frac{\gamma_{0}}{2}+\frac{\gamma_{1}}{2}-3\gamma_{2}\right)q^{4}.$

We show the phase portraits associated with Hamiltonian (2.1) for values of
$v=0$ and $v=10$ with fixed $M_{\rm d}=1$ with $u=1$ in the first column of
Figure 2.1. In the reduced system, the ground state standing wave becomes a
fixed point at the origin, and the excited state standing wave becomes the
boundary circle $p^{2}+q^{2}=1$. For $v=10$, two new asymmetric states have
emerged from the odd solution and appear as fixed points on the $q$-axis. The
right column shows the standing waves constructed from the Galerkin ansatz,
including the initial and desired states of the control problem $\psi_{0}$ and
$\psi_{\rm d}$.

Figure 2.1. Left panels: the phase space of system (2.1) for $v=0$ (top) and
$v=10$ (bottom), with the ground state marked with a black star, showing that
two fixed points have bifurcated from the bounding circle as $v$ was raised,
now sitting at approximately $(q,p)=(\pm 0.7550,0)$. Right panels, same
parameter values, show in black the ground state standing wave and in red/blue
the odd-symmetric standing waves (top) and the symmetry broken standing waves
$\psi_{0}$ and $\psi_{\rm d}$ (bottom).

For $v=0$, the double-well structure is absent (i.e., the the setting is one
of a parabolic trap with equidistant linear eigenvalues), hence the two-mode
reduction is not expected to provide an adequate representation of the
dynamics (except for very low masses $M_{\rm d}$). Of course, to find the
coefficients $c_{0}$ and $c_{1}$, we back substitute using the above canonical
transformations through

(2.12) $c_{0}=\sqrt{M_{\rm d}-q^{2}-p^{2}}e^{i\theta},\quad
c_{1}=(q+ip)e^{i\theta},$

where $\theta\in[0,2\pi]$ is arbitrary by the phase invariance of the
Hamiltonian (2.1). Without loss of generality, we choose $\theta=0$ so that
$\psi_{2}^{\rm Gal}$, given by Equation (2.5), is real.

Hereafter, we fix the barrier height to $v=10$, and choose the value of
$M_{\rm d}$ such that the stable fixed points with $q\neq 0$ are asymmetric
states, as shown in Figure 2.1. That is, we operate within the symmetry-broken
regime of the double-well potential. We note that for any finite-strength
barrier, the stable state is “partially fragmented” in that a nonzero fraction
of the mass resides in each well; see [41] for a theoretical account within
the many-body formalism of BEC. As the value of $v$ is further increased, the
fragmentation is lessened, which provides an obvious strategy for mitigating
fragmentation.

Moving forward, the phase space shown in (the bottom panels of) Figure 2.1
makes clear the goal of our optimal control problem: to find functions $u(t)$
and $v(t)$, with fixed and identical initial and terminal conditions, that
drive the system from one asymmetric steady state $\psi_{0}(x)$ to the other
one $\psi_{\rm d}(x)$. From a physical perspective, our aim is to drive atoms
from a state in which most reside in one well, into one in which most reside
in the other, using the experimentally-developed ability to temporally drive
double-well potentials [43] and more specifically magnetic and optical
confining beams [15, 37, 38]. This control problem is mathematically
formulated in Section 3.

#### Numerical Validation

Before describing the optimal control problem in detail, we numerically test
the ability of the two-mode system (2.6) to approximate the dynamics of the
GPE (2.1) with appropriate initial conditions. We consider the evolution of
the initial condition $\psi(x,0)=\varphi_{0}(x)$ as shown in Figure 2.1
subject to the GPE with imposed controls

(2.13) $u_{\rm trial}(t)=1+2\sin\left(\frac{\pi t}{T}\right)\mbox{\quad
and\quad}v_{\rm trial}(t)=10\cos^{4}\left(\frac{\pi t}{T}\right)$

over the time interval $t\in[0,T]$, with $T=2.$ We have proposed these trial
controls based on ad hoc reasoning, using on the following partial intuition:
before performing an optimization, we suspect that an optimal potential would
allow the barrier at the origin to lower so that the initial mass on the one
well can be transferred significantly to the other well, followed by a raising
of the barrier anew in order to localize the wavefunction into the desired
well of the potential. The choices in the modulation of the parabolic trap,
the strength of the localized barrier, and the length of time in the
simulation were all chosen arbitrarily in this dynamical example.

It does not escape us that in physical settings involving quasi-1d double
wells in atomic BECs, the width parameter $u$ is constrained to be $u(t)\ll 1$
for the quasi-1d reduction to hold. We have considered such scenarios as well,
finding qualitatively similar results, as regards the optimal control
framework discussed later on in Section 3 within the Galerkin truncation,
although over considerably longer time scales.

We numerically integrate the two-mode system (2.6) using Matlab’s ode45. Its
time-dependent coefficients from Eq. (2.7) depend on $\mathbf{w}(t)$ through
the instantaneous eigenfunctions $\varphi_{0}$ and $\varphi_{1}$. While closed
form expressions for these eigenfunctions are determined, for each value of
$\mathbf{w}$, in terms of hypergeometric functions [47], we find it simpler to
solve the associated eigenproblem numerically at each time step using Matlab’s
eig command. Indeed, while the former possibility is particular to the
potential considered herein, the latter can be extended to arbitrary time-
dependent potentials.

We solve the GPE (2.1) using a second-order Fourier split-step method and
approximate the delta function by a narrow Gaussian

(2.14) $\delta(x)=\lim_{a\to\infty}\frac{a}{\sqrt{\pi}}e^{-a^{2}{x^{2}}},$

with $a=12$ here and in all subsequent computations. Although the Fourier
split-step method we use to solve the GPE is quite standard, we provide
details about the method in A for completeness. Throughout this work, we find
that a uniform spatial discretization on the truncated interval
$x\in[-5\pi,5\pi]$ of $2^{11}$ points and a temporal discretization of $T/h$
points, choosing $h=2^{-8},$ yields an accurate and stable computation of the
GPE dynamics.

In what follows, we need a way to measure the agreement between the solutions
to GPE and the finite-dimensional approximation defined by expansion (2.5). We
define the projected wavefunction as

(2.15) $\psi^{\rm proj}_{N+1}=\sum_{n=0}^{N}\left<\psi^{\rm
GPE}(\cdot,t),\varphi_{n}(\cdot,u(t))\right>\varphi_{n}(x,u(t)),$

where $\psi^{\rm GPE}$ solves Equation (2.1). We quantify an expected upper
bound on the Galerkin approximation (2.5) through the relative error

(2.16) $\mathcal{E}_{N+1}(t)=\frac{\norm*{\psi^{\rm GPE}-\psi_{N+1}^{\rm
proj}}^{2}}{M}=\frac{\left\lVert\sum_{n=N+1}^{\infty}\left\langle\psi^{\rm
GPE}(\cdot,t),\varphi_{n}(\cdot,\mathbf{w}(t))\right\rangle\varphi_{n}(x,\mathbf{w}(t))\right\rVert^{2}}{M},$

where $M$ is defined by Eq. (2.2) and the norm is taken in
$L^{2}({\mathbb{C}})$. The second equality, which is interpreted as the
relative mass content which has been excited beyond the few-mode
representation at order $N$, relies on the fact that the GPE (2.1) is well-
approximated by Formula (2.15) in $L^{2}\left({\mathbb{C}}\right)$, i.e.,
$\mathcal{E}_{N+1}\to 0$ as $N\to\infty.$ Therefore, we rely on various
comparisons among Equations (2.5), (2.15), and (2.16) when discussing the
extent to which few-mode representations effectively shadow the full-dynamical
picture given by the GPE (2.1). Also, note that although wavefunctions are
normalized to have unit mass throughout this work, we include $M$ in the above
definition of $\mathcal{E}_{N+1}$ to maintain clarity of how the analysis
should be performed given a different scaling of Equation (2.1).

With those preliminaries, we are ready to show the results of the simulations.
Figure 2.2 shows simulations the GPE (2.1) and the two-mode model system
(2.6). It presents three colormaps: the wave function $\psi^{\rm GPE}$, the
projection of the wave function onto the first two instantaneous
eigenfunctions, i.e., $\psi^{\rm Gal}_{2}$, and the wave function
$\psi_{N+1}^{\rm Gal}$ constructed from the solution to the two-mode model
(2.6) using Formula (2.5). A fourth plot shows the solutions of Equation (2.6)
as well as the coefficients defining the projected wavefunction (2.15) from
the GPE solution.

Figure 2.2. A comparison between numerical solutions to the GPE (2.1) and the
two-mode system given by Equations (2.6). Top left: a full numerical
simulation of the GPE (2.1) with the colormap showing the squared amplitude.
Top right: the projected wavefunction $\psi_{2}^{\rm proj},$ in absolute value
squared, as in Equation (2.15). Bottom left: a comparison of numerically
computed and projected coefficients $c_{n}(t)$ for $n=0$ and $n=1$, as well as
the loss of mass to higher modes given by the error formula (2.16). Bottom
right: the wavefunction $\psi_{2}^{\rm Gal}$, as defined by Equation (2.5) and
also in absolute value squared, with coefficients shown in the bottom left
panel.

We make two observations based on these plots. First, our naively chosen
control functions (2.13) crudely transfer the bulk of the solution from the
left potential well to the right, in both the GPE system and the two-mode
model system. Second, and despite the first observation, the agreement between
the two dynamics is poor. This is seen in the poor agreement to the computed
values of the coefficients and through the large relative error
$\mathcal{E}_{2}$(t).

This is a central observation of this study: while the two-mode reduction
effectively describes the bifurcation structure and the dynamical evolution in
the vicinity of the symmetric and asymmetric equilibrium. However, for the
more highly non-equilibrium transfer proposed herein, a representation
requiring more modes becomes necessary. In that light, we now pursue a three-
mode reduction of the system.

###  2.2. Three-Mode Expansion

In order to increase the fidelity of the reduced model to the full GPE, we now
compute a model equation using $N=2$ in the expansion (2.5). The form of the
Hamiltonian system is long and fairly unenlightening, so we display the
associated Hamiltonian (B.1) in Appendix B. To simplify the search for
stationary solutions, we perform a canonical transformation similar to Eq.
(2.12), yielding Hamiltonian (B). Not surprisingly, we find asymmetric
standing waves similar to $\psi_{0}$ and $\psi_{\rm d}$ from Figure 2.1. In
fact, the contribution to these modes from $c_{2}$ is no greater than one part
in $10^{-10}$ in absolute-value squared. This numerically justifies the
simplification of using the fixed point, corresponding to the left-asymmetric
state from Figure 2.1, with no contribution from the third mode. Additionally,
this also allows for consistency in testing the accuracy of the three-mode
model.

Using the above-mentioned initial condition, we run the same test shown in
Figure 2.2, but for the three-mode model, as shown in Fig. 2.3. In this case,
the approximation is considerably more accurate, especially up to intermediate
times. For example at $t=1$, the relative error $\mathcal{E}_{3}$ is about
0.05. More importantly, we see that the even modes in the reduction
effectively capture the even-projected dynamics of the GPE (2.1) for up to
intermediate times.

Figure 2.3. A comparison between numerical solutions to the GPE (2.1) and the
three-mode system given by Hamiltonian (B.1) as was similarly presented in
Figure 2.2. Here, we see a substantial reduction in the error
$\mathcal{E}_{3}(t)$ in contrast with the error $\mathcal{E}_{2}(t)$ shown in
Figure 2.2, allowing the three-mode model to more effectively capture the full
dynamics given by Equation (2.1).

However, we lose a great deal of accuracy for the remainder of the simulation,
where, by $t=T$ the relative error has grown to about 0.40. The plot of
$c_{1}$ and $c_{1}^{\rm proj}$ shows that this is due to a significant
excitation of higher-order odd modes which have been neglected. Despite this,
we find that the three-mode model performs well enough in our pursuit of
optimal controls. In Section 4, we quantify the contribution of the higher
modes to the relative error (2.16) while the condensate is being controlled in
the full dynamical setting. We find that the inclusion of just one more odd
mode substantially reduces the relative error. There, we discuss this in
greater detail, and, for now, leave the pursuit of higher-dimensional models
as a subject for future work.

We make a final comment on the role that the mass $M$ has on the relative
error $\mathcal{E}_{3}$. Indeed, for smaller values of the mass $M$, we see a
reduction in $\mathcal{E}_{3}$, as expected since nonlinear interactions in
Equation (2.1) are substantially smaller. Using $M=0.2$ which is only slightly
above the symmetry-breaking bifurcation value, we find $\mathcal{E}(T)$ to be
near 0.25. Although this error has been reduced, enhancing the efficacy of the
reduced dynamical system, the mode $\psi_{0}$ corresponding to the asymmetric
stationary state of Hamiltonian (2.1) is only slightly asymmetric, failing to
meet our stated goal of concentrating the mass in a single well. For this
reason, we continue to use $M=1$ in what follows, since this strikes a good
balance between a small enough relative error and an adequate initial
condition $\psi_{0}$.

##  3\. Optimal Control Strategy

The optimal control problem we pursue, motivated by the previous section, is
to construct a function $\mathbf{w}(t)$ that drives an initial state
$\mathbf{c}_{0}$ to the desired state $\mathbf{c}_{\rm d}$ under the dynamics
of the two- or three-mode model. More precisely, the control problem is to
find a local minimizer of the problem

(3.1)
$\min_{\mathbf{w}(t)\in\mathcal{W}}\mathcal{J}=\min_{\mathbf{w}\in\mathcal{W}}\left\\{M^{2}_{\rm
d}-\absolutevalue*{\left\langle\mathbf{c}_{\rm
d},\mathbf{c}_{T}\right\rangle}^{2}\right\\},$

where the admissible space of controls is given by
$\mathcal{W}=\\{\mathbf{w}(t)\in
C^{0}([0,T]):\mathbf{w}(0)=\mathbf{w}(T)=\mathbf{w}_{b}\\}$, for a prescribed
value of $\mathbf{w}_{b}$. We call this objective functional the discrete
infidelity, as opposed to the full infidelity

(3.2) $\mathcal{J}_{\rm full}=M^{2}-\absolutevalue*{\left\langle\psi_{\rm
d},\psi_{T}\right\rangle}^{2}$

first used by Hohenester, et al. [15]. In both cases, the infidelity penalizes
misalignments of the final computed state with the desired state. Note that
since the discrete mass $M_{\rm d}$ is conserved by the dynamics governing
$\mathbf{c}$, the optimal, perhaps unachievable, final state is the desired,
which yields an optimal infidelity of zero.

An advantage of the infidelity (3.1) is that it is insensitive to the global
phase of the dynamics, which is physically unimportant. A more traditional
least-squares approach can be introduced via the objective

(3.3) $\mathcal{J}_{\rm d}^{\rm LSQ}=\frac{1}{M_{\rm
d}}\absolutevalue*{\mathbf{c}_{\rm d}-e^{iS}\mathbf{c}_{T}}^{2},$

where the phase $S\in[0,2\pi].$ Although we do not pursue such an optimization
we report the values of both $\mathcal{J}_{\rm d}^{\rm LSQ}$ and the value of

(3.4) $\mathcal{J}_{\rm full}^{\rm LSQ}=\frac{1}{M}\norm*{\psi_{\rm
d}-e^{iS}\psi_{T}}^{2},$

after a minimization over $S$, as a relative, and more familiar, measure of
optimality at the level of the full-dynamical picture.

The optimal control problem (3.1) is posed over the admissible space
$\mathcal{W}$, which is infinite-dimensional. We approximate this by a finite-
dimensional admissible space, constructed using a Galerkin-type method called
the chopped random basis method (CRAB), first used by [8, 6] and explained in
great detail in the work of [1]. We use the following basis and trial
functions

(3.5) $\mathbf{w}_{\rm CRAB}=\mathbf{w}_{\rm
trial}+\mathbf{w}_{b}\sum_{j=1}^{N_{D}}\frac{\epsilon_{j}^{w}}{j^{2}}\sin{\left(\frac{j\pi
t}{T}\right)},$

where the value of $\mathbf{w}_{b}$ is consistent with the boundary conditions
implied by the trial controls $u_{\rm trial}$ and $v_{\rm trial}$ in Equation
(2.13).

The amplitudes $\varepsilon_{j}^{w}$ are random variables drawn uniformly from
$[-1,1]$. We choose the coefficients $A_{j}=j^{-2}\epsilon_{j}^{w}$ to decay
quadratically because the Fourier series of absolutely continuous functions
exhibit the same type of decay [45]. In this way, the search space for optimal
control ${\bf w}$ is not severely restricted, yet candidate controls remain
technically feasible. To find the coefficients $\varepsilon_{j}^{w}$, we use
the differential evolution (DE) method [42] outlined in C.

#### Remark

The numerical optimization problem associated with the objective functional in
Equation (3.1) is often stated in the variational form of Euler-Lagrange
equations and solved using a form of gradient descent [4, 1]. This requires
functional derivatives of the objective $\mathcal{J}$ with respect to the
control vector $\mathbf{w}$ and thus involves derivatives of the basis
functions $\varphi_{n}$ in Equation (2.5). Since these derivatives cannot be
written in closed form, this renders gradient-based methods cumbersome. Thus
we choose not to pursue such a strategy here. In previous work [1], we
numerically solve a similar optimization problem using both the CRAB method,
and when possible, a combination of the CRAB method and gradient descent and
find that the CRAB method alone is fairly successful in finding efficient
control policies.

##  4\. Results of Numerical Optimization

In this section, we present the results of our numerical optimization. We
briefly summarize the strategy outlined over the past sections. We perform the
optimization on the three-mode reduced system described by the Hamiltonian
(B.1). We use the initial and desired profiles for the two-mode system shown
in Fig. 2.1, as we found that $c_{2}$ was negligibly small in the
corresponding stationary solutions of the three-mode system. We minimize the
phase-insensitive infidelity $\mathcal{J}$ defined in Eq. (3.1). The control
function $\mathbf{w}(t)$ is constructed using the CRAB method (3.5), and the
optimization is performed using the DE method described in C. We find, through
trial and error, that $N_{D}=15$ basis functions are sufficient for the CRAB
approximation. That is, a CRAB basis (3.5) with this many modes reduces the
optimal control problem (3.1) to a nonlinear programming problem in $2N_{D}$
dimensions whose solution we find acceptable.

Figure 4.1 shows the result of this numerical optimization. This computation
yielded an objective function $\mathcal{J}=0.0501$, as defined by Equation
(3.1), and a least squares infidelity $\mathcal{J}_{\rm LSQ}=0.0422$, as
defined by Equation (3.3). The figure also shows the full infidelity
$\mathcal{J}_{\rm full}$ and the full least-squares objective
$\mathcal{J}^{\rm LSQ}_{\rm full}$ resulting from using optimal controls in
simulating the GPE (2.1).

It is important to recognize that the dynamics selected by the optimizer
somewhat conforms to the anticipated physical intuition about the optimal
strategy. It can indeed be observed in the figure that the barrier height
$v(t)$ decreases to about $0$, while the parabolic confinement becomes
tighter, which enables the mass to be transferred from one side of the lowered
barrier to the other. This transfer is visible at time $t\approx 1$ in the the
top right panel of the figure. Subsequently, $v(t)$ rises sharply again, as
the parabolic confinement returns to its original value, so that the
combination of the two now ensures confinement of the transmitted mass to the
right well. It is through this procedure that the infidelity is substantially
decreased in the lower left panel, and indeed subsequently remains small,
during the return of the confinement conditions to their original settings.

Figure 4.1. The result of using numerical optimal control theory. Top left:
optimal controls identified via our numerical methodology. Top right: The
solution of the GPE (2.1) in absolute-value squared. The dotted white line, at
$T=2$, represents the moment the controls are held at their constant terminal
values $\mathbf{w}_{b}$. Bottom left: the full infidelity $\mathcal{J}_{\rm
full}$ and full modified least-squares objective $\mathcal{J}^{\rm LSQ}_{\rm
full}$. Bottom right: wavefunction profiles, in absolute-value squared, of the
initial state $\psi_{0}$, the desired state $\psi_{\rm d}$ the state
$\psi^{\rm Gal}$ computed via the three-mode model at $T=2$, and the state
$\psi^{\rm GPE}$ computed via the GPE (2.1) at $T=2.$

We also show, in Figure 4.2, an error analysis. The left panel in Figure 4.2
is equivalent to the lower-left panel of Figure 2.3, showing the coefficients
$c_{j}(t)$ and the relative error $\mathcal{E}_{3}(t)$ defined by Equation
(2.16). While the $\max_{t}\mathcal{E}_{3}(t)$ is lower in the simulation
using the optimal control, than in the trial control, it is interesting to
note that using the optimal control yields significantly better agreement
between the dynamics of the full GPE system and its projection onto the first
three modes at the final time $T$, decreasing from a maximum of about 18% to
just about 5% at time $T$. The right panel of the figure computes the error
$\mathcal{E}_{N}(t)$ for values of $N\leq 6$. We find $\mathcal{E}_{N}(t)$
decreases monotonically, pointwise in $t,$ as expected. Although we have shown
that three modes suffice to control the GPE, this result indicates the
expected accuracy for $N>3$.

Figure 4.2. Left: Numerical solution of the three-mode model, projected
coefficients, and the error consistent with Figure 4.1. Right: the error of
modes two through six, with the natural number $N$ used consistently with the
Galerkin reduction defined by Equation (2.5).

##  5\. Conclusions

In the present work, we have applied the methodology of optimal control to
dynamics in double-well potentials, one of the most prototypical (and highly
controlled) experimental settings, both in atomic Bose-Einstein condensates
and in nonlinear optics. We have adapted the relevant methodology to the well-
established description of Galerkin truncations within this setting.
Typically, two-mode truncations are used in order to explore the steady
states, stability, bifurcations, and dynamics of such double-well systems. A
key finding of the present work is that such a two-mode approximation is not
sufficient in order to characterize the dynamics associated with optimal
transport.

Indeed, it was found that the involvement of a three-mode approximation was
crucial in order to adequately describe the relevant dynamical process,
identified as optimal by the differential evolution algorithm deployed herein.
Despite quantitative inaccuracies in using a low-order Galerkin reduction, we
see that controlling the three-mode model effectively controls the GPE
wavefunction in numerical simulations. Furthermore, our numerical results
showcase the extent to which the model can be refined by taking the Galerkin
reduction out to higher orders (recall that the basis we use is complete in
$L^{2}$, hence the Galerkin representation converges to the wavefunction
$\psi$ satisfying the GPE in $L^{2}$ as the number of modes N increases). In
particular, our computations, shown in Figure 4.2, indicate that substantial
gain can be made by including a fourth mode into the model. We leave the
potential inclusion of these higher modes in a reduction-based optimal control
strategy as a subject for future work.

Additionally, the adaptation of such optimal control methodologies to low-
dimensional truncated Galerkin dynamics is a technique that could find
significant potential for further applications. Some possibilities include
multi-component and spinor condensates [22] where few-mode approximations have
proven useful [48, 13, 46]. Moreover, extending such control strategies beyond
the mean-field framework and into the realm of many-body effects [31, 36], is
of particular interest in its own right; for a review of the latter, see,
e.g., the recent preprint of [33]. In regards to complementary or alternative
low-dimensional control strategies, it should be possible to couple the
optimal control methods used here with methods based on the so-called shortcut
to adiabaticity [12]. We additionally note that other methods more tailored to
the driving of specific potentials [34], as well as ones that are model
agnostic [16] also exist in the literature and are of interest in their own
right.

It is also natural to extend these ideas to higher-dimensional systems where
few-well arrangements have also been explored [49] or to systems with
disorder, such as recent work which makes use of deep learning [17]. The
Galerkin approach here may be generalized to systems with disorder so that
optimal control policies can be computed via a low-dimensional dynamic
program. Such studies are currently under consideration and will be presented
in future publications.

## Acknowledgments

This material is based upon work supported by the US National Science
Foundation under Grants No. DMS-1809074 and PHY-2110030 (P.G.K.).

## Disclosures

The authors declare no conflicts of interest.

## Data Availability

Data underlying the results presented in this paper are not publicly available
at this time but may be obtained from the authors upon reasonable request.

## References

  * [1] J. Adriazola and R. H. Goodman, Reduction-based strategy for optimal control of Bose-Einstein Condensates, Phys. Rev. E, 105 (2022), p. 025311.
  * [2] M. Albiez, R. Gati, J. Fölling, M. C. S. Hunsmann, and M. K. Oberthaler, Direct observation of tunneling and nonlinear self-trapping in a single bosonic Josephson junction, PRL, 95 (2005), p. 010402.
  * [3] W. H. Aschbacher, J. Fröhlich, G. M. Graf, K. Schnee, and M. Troyer, Symmetry breaking regime in the nonlinear Hartree equation, J. Math. Phys., 43 (2002), pp. 3879–3891.
  * [4] A. Borzi, G. Ciaramella, and M. Sprengel, Formulation and Numerical Solution of Quantum Control Problems, SIAM, Philadeplphia, 2017.
  * [5] C. Brif, R. Chakrabarti, and H. Rabitz, Control of quantum phenomena: past, present and future, New J. Phys., 12 (2010), p. 075008.
  * [6] T. Calarco, E. A. Hinds, D. Jaksch, J. Schmiedmayer, J. I. Cirac, and P. Zoller, Quantum gates with neutral atoms: Controlling collisional interactions in time-dependent traps, Phys. Rev. A, 61 (2000), p. 022304.
  * [7] C. Cambournac, T. Sylvestre, H. Maillotte, B. Vanderlinden, P. Kockaert, P. Emplit, and M. Haelterman, Symmetry-breaking instability of multimode vector solitons, PRL, 89 (2002), p. 083901.
  * [8] T. Caneva, T. Calarco, and S. Montangero, Chopped random-basis quantum optimization, Phys. Rev. A, 84 (2011), p. 022326.
  * [9] R. Glowinski, S. J. Osher, and W. Yin, Splitting Methods in Communication, Imaging, Science, and Engineering, Springer International Publishing AG, Cham, 2017.
  * [10] R. H. Goodman, Bifurcations of relative periodic orbits in NLS/GP with a triple-well potential, Phys. D, 359 (2017), pp. 39–59.
  * [11] R. H. Goodman, J. L. Marzuola, and M. I. Weinstein, Self-trapping and Josephson tunneling solutions to the nonlinear Schrödinger / Gross-Pitaevskii equation, Disc. Cont. Dyn. Sys. A, 35 (2015), pp. 225–246.
  * [12] D. Guéry-Odelin, A. Ruschhaupt, A. Kiely, E. Torrontegui, S. Martínez-Garaot, and J. G. Muga, Shortcuts to adiabaticity: Concepts, methods, and applications, Rev. Mod. Phys., 91 (2019), p. 045001.
  * [13] M. Günay, Binary mixture of Bose-Einstein condensates in a double-well potential: Berry phase and two-mode entanglement, Phys. Rev. A, 101 (2020), p. 043608.
  * [14] U. Hohenester, OCTBEC—A Matlab toolbox for optimal quantum control of Bose–Einstein condensates, Computer Physics Communications, 185 (2014), pp. 194–216.
  * [15] U. Hohenester, P. K. Rekdal, A. Borzi, and J. Schmiedmayer, Optimal quantum control of Bose-Einstein condensates in magnetic microtraps, Phys. Rev. A, 75 (2007), p. 023602.
  * [16] T. Huang, Y. Ban, E. Y. Sherman, and X. Chen, Machine-learning-assisted quantum control in a random environment, Phys. Rev. Appl., 17 (2022), p. 024040.
  * [17] , Machine-learning-assisted quantum control in a random environment, Phys. Rev. Applied, 17 (2022), p. 024040.
  * [18] R. K. Jackson and M. I. Weinstein, Geometric analysis of bifurcation and symmetry breaking in a Gross-Pitaevskii equation, J. Statist. Phys., 116 (2004), pp. 881–905.
  * [19] G. Jäger, D. M. Reich, M. H. Goerz, C. P. Koch, and U. Hohenester, Optimal quantum control of Bose-Einstein condensates in magnetic microtraps: Comparison of gradient-ascent-pulse-engineering and krotov optimization schemes, Phys. Rev. A, 90 (2014), p. 033628.
  * [20] J. Joannopoulos, S. Johnson, J. Winn, and R. Meade, Photonic Crystals: Molding the Flow of Light, Princeton University Press, Princeton, 2008\.
  * [21] T. Kapitula, P. G. Kevrekidis, and Z. Chen, Three is a crowd: solitary waves in photorefractive media with three potential wells, SIAM J. Appl. Dyn. Syst., 5 (2006), pp. 598–633.
  * [22] Y. Kawaguchi and M. Ueda, Spinor Bose–Einstein condensates, Physics Reports, 520 (2012), pp. 253–381.  Spinor Bose–Einstein condensates.
  * [23] P. Kevrekidis, D. Frantzeskakis, and R. Carretero-González, The defocusing nonlinear Schrödinger equation, SIAM, Philadelphia, PA, 2015\.  From dark solitons to vortices and vortex rings.
  * [24] P. G. Kevrekidis, Z. Chen, B. A. Malomed, D. J. Frantzeskakis, and M. I. Weinstein, Spontaneous symmetry breaking in photonic lattices : Theory and experiment, Phys. Lett. A, 340 (2005), pp. 275–280.
  * [25] D. E. Kirk, Optimal control theory: An introduction, Courier Corporation, 2004.
  * [26] E. Kirr, P. G. Kevrekidis, and D. E. Pelinovsky, Symmetry-breaking bifurcation in the nonlinear Schrödinger equation with symmetric potentials, Comm. Math. Phys., 308 (2011), pp. 795–844.
  * [27] E. W. Kirr, P. G. Kevrekidis, E. Shlizerman, and M. I. Weinstein, Symmetry-breaking bifurcation in nonlinear Schrödinger/Gross-Pitaevskii equations, SIAM J. Math. Anal., 40 (2008).
  * [28] Y. S. Kivshar and G. P. Agrawal, Optical Solitons: From Fibers to Photonic Crystals, Academic Press, 2003.
  * [29] F. Lederer, G. I. Stegeman, D. N. Christodoulides, G. Assanto, M. Segev, and Y. Silberberg, Discrete solitons in optics, Physics Reports, 463 (2008), pp. 1–126.
  * [30] J. L. Marzuola and M. I. Weinstein, Long time dynamics near the symmetry breaking bifurcation for nonlinear Schrödinger/Gross-Pitaevskii equations, Discrete Contin. Dyn. Syst., 28 (2010), pp. 1505–1554.
  * [31] D. Masiello and W. Reinhardt, Symmetry-broken many-body excited states of the gaseous atomic double-well Bose–Einstein condensate, J. Phys. Chem. A, 123 (2019), pp. 1962–1967.
  * [32] J. Mennemann, D. Matthes, R. Weishaupl, and T. Langen, Optimal control of Bose-Einstein condensates in three dimensions, New J. Phys., 17 (2015), p. 113027.
  * [33] S. I. Mistakidis, A. G. Volosniev, R. E. Barfknecht, T. Fogarty, T. Busch, A. Foerster, P. Schmelcher, and N. T. Zinner, Cold atoms in low dimensions–a laboratory for quantum dynamics, arXiv:2202.11071, (2022).
  * [34] B. J. Mommers, A. Pritchard, T. A. Bell, R. N. Kohn, S. E. Olson, M. Baker, and M. W. J. Bromley, Atomtronic many-body transport using husimi driving, 2021.
  * [35] O. Morsch and M. Oberthaler, Dynamics of Bose-Einstein condensates in optical lattices, Rev. Mod. Physics, 78 (2006), pp. 179–215.
  * [36] P. Mujal, B. Juliá-Díaz, and A. Polls, Quantum properties of a binary bosonic mixture in a double well, Phys. Rev. A, 93 (2016), p. 043619.
  * [37] C. J. Pethick and H. Smith, Bose-Einstein condensation in dilute gases, Cambridge University Press, Cambridge, 2003.
  * [38] L. Pitaevskii and S. Stringari, Bose-Einstein condensation, Oxford University Press, Oxford, 2003.
  * [39] S. Raghavan, A. Smerzi, S. Fantoni, and S. Shenoy, Coherent oscillations between two weakly coupled Bose-Einstein condensates: Josephson effects, $\pi$ oscillations, and macroscopic quantum self-trapping, Phys. Rev. A, 59 (1999), pp. 620–633.
  * [40] A. Sacchetti, Universal critical power for nonlinear Schrödinger equations with a symmetric double well potential, PRL, 103 (209), p. 194101.
  * [41] R. W. Spekkens and J. E. Sipe, Spatial fragmentation of a Bose-Einstein condensate in a double-well potential, Phys. Rev. A, 59 (1999), pp. 3868–3877.
  * [42] R. Storn and K. Price, Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces, J. Global Optim., 11 (1997), pp. 341–359.
  * [43] I. Stroescu, D. Hume, and M. Oberthaler, Dissipative double-well potential for cold atoms: Kramers rate and stochastic resonance, PRL, 117 (2016), p. 243005.
  * [44] G. Theocharis, P. G. Kevrekidis, D. J. Frantzeskakis, and P. Schmelcher, Symmetry breaking in symmetric and asymmetric double-well potentials, Phys. Rev. E, 74 (2006), p. 056608.
  * [45] L. N. Trefethen, Spectral methods in MATLAB, SIAM, New York, NY, 2000\.
  * [46] J. Tuan and H. Qiu, A dynamical phase transition of binary species becs mixtures in a double well potential, Int. J. Theor. Phys., 56 (2017), pp. 321–327.
  * [47] J. Viana-Gomes and N. M. R. Peres, Solution of the quantum harmonic oscillator plus a delta-function potential at the origin: the oddness of its even-parity solutions, Eur. J. Phys, 32 (2011), p. 1377.
  * [48] C. Wang, P. Kevrekidis, N. Whitaker, and B. Malomed, Two-component nonlinear Schrödinger models with a double-well potential, Physica D: Nonlinear Phenomena, 237 (2008), pp. 2922–2932.
  * [49] C. Wang, G. Theocharis, P. G. Kevrekidis, N. Whitaker, K. J. H. Law, D. J. Frantzeskakis, and B. A. Malomed, Two-dimensional paradigm for symmetry breaking: The nonlinear Schrödinger equation with a four-well potential, Phys. Rev. E, 80 (2009), p. 046611.
  * [50] T. Zibold, E. Nicklas, C. Gross, and M. Oberthaler, Classical bifurcation at the transition from Rabi to Josephson dynamics, PRL, 105 (2010), p. 204101.

##  Appendix A Numerical Method for Solving the Gross-Pitaevskii Equation

It is necessary to solve the GPE (2.1) in order to validate the Galerkin
reduction (2.5) and to evaluate the performance of the dimensionally reduced
optimal control problem. The boundary conditions are assumed to be periodic so
that the use of spectral methods is straightforward. We use an operator
splitting method, and, to this end, rewrite Equation (2.1) in the form

(A.1) $i\partial_{t}\psi=\mathcal{L}\psi+\mathcal{N}(\psi),$

where the linear operator is given by
$\mathcal{L}=-\frac{1}{2}\partial_{x}^{2}$ and the nonlinear and inhomogeneous
operator $\mathcal{N}(\psi)$ incorporates the remaining terms.

We choose to use a second order in time operator splitting, often referred to
as Strang splitting [9], to approximate the resulting matrix exponential by

(A.2)
$e^{\left(\mathcal{L}+\mathcal{N}\right)Mh}=e^{\mathcal{L}h/2}e^{\mathcal{N}h}e^{\mathcal{L}h}\ldots
e^{\mathcal{L}h}e^{\mathcal{N}h}e^{\mathcal{L}h/2}+\mathcal{O}\left(h^{2}\right),$

where $h=T/M$ is the time discretization for a given number of time steps $M$.
The solution of the linear equation resulting from the matrix exponential
$e^{\mathcal{L}h}$ is facilitated by the Fourier transform and is given by

(A.3)
$\psi_{n+1}=\mathcal{F}^{-1}\left\\{\mathcal{F}\\{\psi_{n}\\}e^{-\frac{ihk^{2}}{2}}\right\\},$

where $\mathcal{F}$ and $\mathcal{F}^{-1}$ denote the Fourier and inverse
Fourier transformation, respectively. Fourier and inverse Fourier transforms
are computed via the fast Fourier transform functions in MATLAB, fftw and
ifftw, with discretized wavenumbers
$k\in\left[\frac{-N/2+1}{2l},\frac{N/2-1}{2l}\right],$ where $2\pi l$ units of
length are assumed in the truncated spatial domain and $N$ is the number of
spatial discretization points.

The nonlinear equation resulting from the matrix exponential
$e^{\mathcal{N}h}$ is quite simple. Since the nonlinearity does not involve
spatial derivatives, we are simply tasked with solving ODEs. Using polar
coordinates, i.e., letting

(A.4) $\psi=\rho(x,t)e^{i\theta(x,t)},$

results in the system

(A.5) $\displaystyle\partial_{t}\rho$ $\displaystyle=0,$
$\displaystyle\partial_{t}{\theta}$ $\displaystyle=-\rho^{2}-V(x,t).$

The first of these equations is a statement about conservation of the mass
$||\psi||_{L^{2}({\mathbb{C}})}$. The second equation, governing the phase
$\theta,$ can be solved via any number of standard numerical ODE techniques;
we simply use the second-order accurate midpoint method. The update for the
phase $\theta,$ in this case, is

(A.6)
$\theta(x,t_{n+1})=\theta(x,t_{n})-h\rho(x,t_{n})^{2}-hV\left(x,t_{n}+\frac{h}{2}\right)$

##  Appendix B Details of the Three-Mode Model

The Hamiltonian for the three-mode model is found to be

(B.1) $\displaystyle\mathcal{H}$
$\displaystyle=\alpha\absolutevalue{c_{0}}^{2}+\beta\absolutevalue{c_{1}}^{2}+\sigma\absolutevalue{c_{2}}^{2}+2\Delta\imaginary\left\\{c_{0}\bar{c}_{2}\right\\}+\frac{\gamma_{0}}{2}\absolutevalue{c_{0}}^{4}+\frac{\gamma_{1}}{2}\absolutevalue{c_{1}}^{4}+\frac{\gamma_{2}}{2}\absolutevalue{c_{2}}^{4}$
$\displaystyle+\gamma_{3}\left(\real\left\\{c_{0}^{2}\bar{c}_{1}^{2}\right\\}+2\absolutevalue{c_{0}}^{2}\absolutevalue{c_{1}}^{2}\right)+\gamma_{4}\left(\real\left\\{c_{0}^{2}\bar{c}_{2}^{2}\right\\}+2\absolutevalue{c_{0}}^{2}\absolutevalue{c_{2}}^{2}\right)+\gamma_{5}\left(\real\left\\{c_{1}^{2}\bar{c}_{2}^{2}\right\\}+2\absolutevalue{c_{1}}^{2}\absolutevalue{c_{2}}^{2}\right)$
(B.2)
$\displaystyle+2\left(\gamma_{6}\absolutevalue{c_{0}}^{2}+\gamma_{7}\absolutevalue{c_{2}}^{2}\right)\real\left\\{c_{0}\bar{c}_{2}\right\\}+2\gamma_{8}\left(\real\left\\{c_{0}\bar{c_{1}}^{2}c_{2}\right\\}+2\real\left\\{c_{0}\absolutevalue{c_{1}}^{2}\bar{c}_{2}\right\\}\right),$

where the projection coefficients are given by

$\displaystyle\alpha=\left\langle\mathcal{L}\varphi_{0},\varphi_{0}\right\rangle,\
\beta=\langle{\mathcal{L}\varphi_{1},\varphi_{1}}\rangle,\
\sigma=\left\langle\mathcal{L}\varphi_{2},\varphi_{2}\right\rangle,\
\Delta=\left\langle\varphi_{0},\dot{\varphi}_{2}\right\rangle,$
$\displaystyle\gamma_{0}=\|\varphi_{0}^{4}\|,\
\gamma_{1}=\|\varphi_{1}^{4}\|,\
\gamma_{2}=\|\varphi_{2}^{4}\|,\gamma_{3}=\left\langle\varphi_{0}^{2},\varphi_{1}^{2}\right\rangle,\
\gamma_{4}=\left\langle\varphi_{0}^{2},\varphi_{2}^{2}\right\rangle,$ (B.3)
$\displaystyle\gamma_{5}=\left\langle\varphi_{1}^{2},\varphi_{2}^{2}\right\rangle,\
\gamma_{6}=\left\langle\varphi_{0}^{3},\varphi_{2}\right\rangle,\
\gamma_{7}=\left\langle\varphi_{0},\varphi_{2}^{3}\right\rangle,\gamma_{8}=\left\langle\varphi_{0},\varphi_{1}^{2}\varphi_{2}\right\rangle.$

from which the analogous system to the two-mode system (2.6) can be derived
easily. The appearance of the term with coefficient $\Delta$ is due to the
time-dependent nature of the basis functions $\varphi_{n}$, and did not appear
in the two-mode system due to the parity of the first two modes. This, as well
as the implicit claim that
$\left\langle\varphi_{0},\dot{\varphi}_{2}\right\rangle=-\left\langle\varphi_{2},\dot{\varphi}_{0}\right\rangle$
may be verified numerically.

As was done with the two-mode system, we use changes of variables and the, now
three-mode, discrete mass to reduce dimensionality. Because we may only reduce
the number of degrees of freedom to two for the three-mode system, we lose the
comfort of visualizing the dynamics via phase portraits as before.
Nevertheless, we pursue a reduction since this helps make identifying fixed
points of the dynamical system simpler. To this end, we begin by using

(B.4) $c_{0}=\eta_{0}e^{i\theta},\quad c_{1}=re^{i\theta},\quad
c_{2}=\eta_{1}e^{i\theta},$

where $\eta_{j}\in{\mathbb{C}}$ and $r\in{\mathbb{R}}.$ We can easily
eliminate $r$ by using the discrete mass once again: $r=\sqrt{M_{\rm
d}-\absolutevalue{\eta_{0}}^{2}-\absolutevalue{\eta_{1}}^{2}}.$ Now,
converting to polar coordinates via

(B.5)
$\eta_{0}=\sqrt{\rho_{0}}e^{i\theta_{0}},\quad\eta_{1}=\sqrt{\rho_{1}}e^{i\theta_{1}},$

we find the following, two-degree of freedom Hamiltonian

$\displaystyle\mathcal{H}=$ $\displaystyle\frac{M_{\rm
d}}{2}\left(2\beta+\gamma_{1}M_{\rm
d}\right)+\frac{\rho_{0}^{2}}{2}\left(\gamma_{0}+\gamma_{1}-2\gamma_{5}\left(\cos\left(2\theta_{0}\right)-2\right)\right)-\frac{\rho_{1}^{2}}{2}\left(\gamma_{1}+\gamma_{2}-2\gamma_{7}\left(\cos\left(2\theta_{1}\right)-2\right)\right)$
$\displaystyle\rho_{0}\left(\alpha-\beta-\gamma_{1}M_{\rm d}+\gamma_{5}M_{\rm
d}\left(2+\cos\left(2\theta_{0}\right)\right)\right)+\rho_{1}\left(\beta+\gamma_{1}M_{\rm
d}-\gamma_{7}M_{\rm
d}\left(2+\cos\left(2\theta_{1}\right)\right)-\sigma\right)$
$\displaystyle+\rho_{0}\rho_{1}\left(\gamma_{1}-\gamma_{5}\left(2+\cos\left(2\theta_{0}\right)\right)+\gamma_{6}\left(2+\cos\left(2\left(\theta_{0}-\theta_{1}\right)\right)\right)-\gamma_{7}\left(2+\cos\left(2\theta_{1}\right)\right)\right)$
$\displaystyle+2\rho_{0}^{3/2}\rho_{1}^{1/2}\left(\gamma_{3}\sin\left(\theta_{0}\right)\sin\left(\theta_{1}\right)-\gamma_{8}\sin\left(\theta_{0}\right)\sin\left(\theta_{1}\right)+\gamma_{3}\cos\left(\theta_{0}\right)\cos\left(\theta_{1}\right)-3\gamma_{8}\cos\left(\theta_{0}\right)\cos\left(\theta_{1}\right)\right)$
$\displaystyle+2\rho_{0}^{1/2}\rho_{1}^{3/2}\left(\gamma_{4}\sin\left(\theta_{0}\right)\sin\left(\theta_{1}\right)-\gamma_{8}\sin\left(\theta_{0}\right)\sin\left(\theta_{1}\right)+\gamma_{4}\cos\left(\theta_{0}\right)\cos\left(\theta_{1}\right)-3\gamma_{8}\cos\left(\theta_{0}\right)\cos\left(\theta_{1}\right)\right)$
(B.6) $\displaystyle+2M_{\rm
d}\gamma_{8}\sqrt{\rho_{0}\rho_{1}}\left(\sin\left(\theta_{0}\right)\sin\left(\theta_{1}\right)+3\cos\left(\theta_{0}\right)\cos\left(\theta_{1}\right)\right).$

##  Appendix C Optimization via Differential Evolution

DE is a stochastic optimization method used to search for candidate solutions
to non-convex optimization problems. The idea behind DE is a so-called genetic
algorithm that draws inspiration from evolutionary genetics. DE searches the
space of candidate solutions by initializing a population set of vectors,
known as agents, within some chosen region of the search space. These vectors
are then randomly mutated into a new population set, or generation.

Result: A vector $z$ mutated from agents in a given generation as required by
the DE Algorithm (2).

Input : 4 distinct members $a,b,c,d$ from the current generation of agents
each with $N$ components, the crossover ratio $R_{C}\in(0,1)$, and weight
$F\in(0,2)$.

for  _j=1:N_ do

Compute a random variable $\mathtt{rand}$;

if  _$\mathtt{rand} <R_{C}$_ then

$z[j]\leftarrow a[j]+F*(b[j]-c[j])$

else

$z[j]\leftarrow d[j]$

end if

end for

Algorithm 1 Differential Evolution Mutation

Result: A vector likely to be globally optimal with respect to an objective
$J$.

Input : A maximum number of iterations $\mathtt{Nmax}$, crossover ratio
$R_{C}\in(0,1)$ and weight $F\in(0,2)$

while  _$\mathtt{counter} <\mathtt{Nmax}$_ do

Generate a population pop of $N_{\rm pop}$ vectors.

for  _$i=1:N_{\rm pop}$_ do

$\mathtt{CurrentMember}\leftarrow\mathtt{Pop}_{i}$;

Choose three distinct vectors $a_{i},b_{i},c_{i}$ different from the vector
$\mathtt{Pop}_{i}$;

Mutate $a_{i},b_{i},c_{i}$, and the $\mathtt{CurrentMember}$ into the mutated
vector $z$ using the mutation parameters $R_{C},F$ and Algorithm 1;

if  _$J(z) <J(\mathtt{CurrentMember})$_ then

$\mathtt{TemporaryPop}_{i}=z$;

end if

end for

$\mathtt{Pop}\leftarrow\mathtt{TemporaryPop}$;

$\mathtt{counter}\leftarrow\mathtt{counter}+1$;

end while

Algorithm 2 Differential Evolution

The mutation operates via two mechanisms: a weighted combination and a
“crossover” which randomly exchanges “traits”, or vector elements, between
agents. The method requires three parameters; the weight $F\in(0,2)$, the
crossover parameter $R_{C}\in(0,1)$, and the size of the population $N_{\rm
pop}$, which, by Algorithm (1), is required to be an integer greater than
three. A pseudo-code illustrating the implementation of the relevant
algorithms is given in Algorithms (1) and (2). Through trial and error, we
find the parameters $F=0.8,\ R_{C}=0.9,\ N_{\rm pop}=20$ work well. DE ensures
that the objective functional $\mathcal{J}$ decreases monotonically with each
generation. As each iteration “evolves” into the next, inferior vectors
“inherit” optimal traits from superior vectors via mutations. DE only allows
mutations that are more optimal with respect to $\mathcal{J}$ to pass to the
next generation. After a sufficient number of iterations, the best vector in
the final generation is chosen as the candidate solution most likely to be
globally optimal with respect to an objective functional.

Genetic algorithms, which require very few assumptions about the objective
functional, are part of a wider class of optimization methods called
metaheuristics. Although metaheuristics are useful for non-convex optimization
problems, they do not provide guarantees about the global optimality of
candidate solutions. Since the algorithm is stopped after a finite number of
iterations, different random realizations return different candidate
optimizers. The results we show are the best among five different
realizations.

In practice, we do not recommend taking fewer realizations since one runs the
risk of computing highly sub-optimal controls, which, indeed, is a generic
issue when solving non-convex optimization problems using stochastic methods.
Of course, taking more realization is expensive, but we find that with the
optimization and physical parameters used throughout this work, five
realizations are sufficient to guarantee the discovery of, at least, a couple
of optimal control policies which are extremely competitive with regards to
the objective functional (3.1). From a computation of 30 realizations, about a
third of the control policies are within 1% of the infidelity returned by the
best control.