# Finite Dimensional Control of Gross-Pitaevskii Equations in Two-Dimensions

Consider the dimensionless Gross-Pitaevskii equation (GPE)
$$
i\partial_t\psi+\frac{1}{2}\nabla^2\psi-V(x,u(t))\psi-|\psi|^2\psi=0, \ \ \
\psi(x,0)=\psi_0(x)\in H^{1}\left(\mathbb{R}^n\right)
$$
with normalization condition $\|\psi_0\|_{L^2(\mathbb{R}^n)}=1.$ We formulate the design problem of finding the time-dependent parameter $u$ that evolves an initial state $\psi_0$ through the GPE into some desired state $\psi_d$ as follows.

Define an admissible class of controls $\mathcal{U}$ such that $u\in H^1([0,T])$ with fixed initial and terminal conditions $u(0)$ and $u(T)$. In work due to Mennemann, and others, they solve the optimization problem
$$
\inf_{u\in \mathcal{U}}{J}=\frac{1}{2}\inf_{u\in \mathcal{U}} \left [J^{\rm infidelity}(u) + J^{\rm regular}(u)\right]
$$
subject to the $n=3$ GPE where
$J^{\rm infidelity}(u)= 1-\left|{\left\langle\psi_d(x),\psi(x,T)\right\rangle}\right|^2_{L^2(\mathbb{R}^3)}$ and $J^{\rm regular}(u)=\gamma\int_0^T|\dot{u}|^2 dt.$ The infidelity penalizes mismatches in the realized state $\psi(x,T)$ and the desired state $\psi_d$ and has a minimum of 0. Sparber, and others, proved that this optimal control problem is well-posed upon adding the Tikhononv regularization $J^{\rm regular}.$

The main computational drawback in this approach is that 3D GPE solves are expensive. Roy and I used a Galerkin reduction to overcome this. We use a Galerkin expansion of the form
$$
\psi(x,t)=\sum_{j=0}^{N}c_n(t)\varphi_j(x;u(t)),
$$
where each of the basis functions $\varphi_j(x;u(t))$ is an instantaneously normalized eigenfunction of the time-independent linear Schrodinger equation. We find a finite-dimensional Hamiltonian system for the time-dependent coefficients
$$
i\dot{c}_j=\partial_{c^{\dagger}_j}\mathcal{H},\qquad i\dot{c}^{\dagger}_j=-\partial_{c_j}\mathcal{H},\qquad j=0,1,...,N
$$
where the Hamiltonian is really only depedent on the form of the potential $V(x,u(t))$ through the projection coefficients involved in the derivation.

We choose an initial condition $c_0$ which is a fixed-point of the system at $t=0$ and choose a desired condition $c_d$ which is a fixed point at $t=T$. The optimal control problem is now
$$
\min_{u(t)\in\mathcal{U}}\mathcal{J} =
\min_{u\in\mathcal{U}}\left\{1-\left|{\langle{c}_{\rm d}},{c_T}\rangle\right|^2\right\}+J^{\rm regular},
$$
subject to the above Hamiltonian system. Roy and I numerically solved this control problem for the quadratic potential $V=\frac{1}{2}ux^2$ and splitting potential $V=\delta(x)u.$ Together with Panos, we solved the problem of pushing a wavefunction from one side of a double well to another via combining the quadratic and splitting potential. Taking $N=2$ worked sufficiently well in all three cases.

We are interested in extending the above strategy to higher spatial dimensions and we propose to begin with dimension two. A schematic for the type of problem we want to start with is one where we squeeze and elongate the condensate along two different directions. Something that acheives something like this
In this case, let us assume the potential is an anisotropic quadratic well of the form $V(x,u(t),v(t))=\frac{u(t)}{2}x^2+\frac{v(t)}{2}y^2.$ The question we ask is the following: How many and which basis states do we use?

After some thought, I think this question has a very simple answer. Let us revisit the 1D case, with a quadratic potential $V=\frac{1}{2}ux^2$. Set $T=2.5,$ $u(0)=u_0=1,$ $u(T)=u_T=100,$ with a linear in time control $u(t)=u_0+(u_T-u_0)t/T$ which interpolates between the final and intial values. The inital condition $\psi_0$ is mostly made up of the ground-state, but to be more precise, it is chosen as a superposition of the first two even Hermite states which correspond to the single fixed point of the Galerkin-reduced system with $N=2$ (we can ignore the odd state by symmetry of the potential). Here is what a numerical simulation of the GPE looks like, in absolute value squared, with this setup:

Let us define a projected wavefunction by
$$\psi_{\rm proj}(x,t)=\sum_{n=0}^1\left\langle\psi_{\rm GPE}(x,t),\varphi_n(x;u(t))\right\rangle\varphi_n(x;u(t))$$
It turns out, that by $t=T, ||\psi_{\rm GPE}(x,T)-\psi_{\rm proj}||_{L^2(\mathbb{C})}$ is only about 0.03. That is, although we've certainly excited higher modes, a Galerkin strategy with two even modes should get the job done. In fact, if we take a look at $|\psi_{\rm proj}(x,t)|^2$
we see that we capture the Rabi oscillation pretty nicely despite discarding so many modes.

The first step in this project would be to repeat this numerical experiment in two dimensions. Let's double the time to $T=5$. In the context of the anisotropic trap, set $u(0)=1,\ u(T/2)=100,\ v(T/2)=100,\ v(T)=1.$ On $t\in[0,T/2],$ $v(t)=v(T/2)$ is constant, and on $t\in[T/2,T]$ $u(t)=u(T/2),$ is constant. During the time intervals the controls vary, they're linear functions which interpolate between their respective inital and final values. We perform a numerical simulation of the 2D GPE, and we quantify the projected infidelity
$$
\mathcal{E}(\tau)=1-\left|{\left\langle\psi_{\rm GPE}(\cdot,\tau),\psi_{\rm proj}(\cdot,\tau)\right\rangle}\right|^2_{L^2(\mathbb{C}^2)}$$
where $\tau\in[T,T+S]$ and $S$ could be something small like one. By trial and error, we find a combination of basis states which achieves a maximum infidelity over this time window which is below some desired value (something like 0.1 meaning at worst during the time $\tau$ we capture 90% of the information). This result would inform of us of where to go next.

## Results of this Numerical Experiment

Carrying out this experiment, with $S=3$, I see that I can capture about 99% of the dynamics using three modes: $\varphi_{00},\ \varphi_{020},\ \varphi_{02}$ where I'm defining the index by
$$
\varphi_{ij}(x,y):=\varphi_i(x)\varphi_j(y)
$$
with each $\varphi_j(z)$ being the $j^{\rm th}$ Hermite function appropriate for this problem. These three modes look like this:

Here's what the dynamics over the control period looks like:


And for when the trapping potential is held constant:


And we see over the constant control phase, the infidelity is quite small
where increasing index in $\mathcal{E}_j$ means we include more modes in the order that they are visually presented.

It seems obvious to me that the optimal control strategy that we used before is going to work just fine in higher dimensions. This seems like a pretty low risk/low reward project though. Is there something more interesting we could do that builds on these ideas? Maybe moving straight to 3D will help sell the merit of the methodology. Something to think about.

## Adding a Twist

It seems like the above problem is too easy for us to be working on. An interesting direction suggested by Panos is to add a rotation to see if there is an optimal way to excite a vortex state such as, using the above notation, the $\varphi_{01}+i\varphi_{10}$ state

To even begin studying this problem, we need to be able to numerically solve the following GPE
\begin{aligned}
i \partial_t \psi&=-\frac{1}{2} \Delta \psi+V(r,t) \psi+|\psi|^2 \psi+i a(t) \partial_\theta \psi \\
& :=L \psi+N \psi \\
\end{aligned}
where the Laplacian in polar coordinates is given by
$$
\Delta \psi=\partial_r^2 \psi+\frac{1}{r} \partial_ r \psi+\frac{1}{r^2} \partial^2_{\theta} \psi \\
$$
and where we assume the potential $V$ is radially symmetric (whether we want to include a control here or not is still to be decided and is not critical in what follows). Note that we have a time-dependent control in front of the term modeling rotation.

For purposes of an operator splitting, for example,

$$
e^{-it(L+N)}\psi_0=e^{-itN/2}e^{-itL}e^{-itN/2}\psi_0
$$
let us focus on
$$
\begin{aligned}
i \partial_t \psi & =L \psi \\
& =-\frac{1}{2} \Delta \psi+i a(t) \partial_\theta \psi
\end{aligned}
$$
since the remaining terms in the nonlinear update are handled easily by use of the Madelung transform.

We can begin considering the numerical evaluation of the operator $e^{-itL}$ by analyzing each harmonic in the azimuthal direction, i.e., by letting
$$
\psi=\sum_{m=-N / 2+1}^{N / 2} e^{i m \theta} \varphi(r, t)
$$

Then for each $m$ we have
$$
\begin{aligned}
& \partial_t \varphi-i m a(t) \varphi=\frac{i}{2}\left(\partial_r^2 \varphi+\frac{1}{r} \partial_r \varphi-\frac{m^2}{2 r^2} \varphi\right)
\end{aligned}
$$

We can remove the last term by performing the change of variables
$$
\phi=e^{-i m \int^t a(s) d s} \varphi
$$

so that
$$
i \partial_t \phi=-\frac{1}{2} \partial_r^2 \phi-\frac{1}{2 r} \partial_r \phi+\frac{m^2}{2 r^2} \phi.
$$
This is great because we can parallelize over the wavenumber $m$ for a substantial speedup. To solve this equation, we have three options. At this point, I'm not sure which to pursue.

### Option 1: Use an integral transform
Define the order $v$ Hankel transform pair by
$$
\begin{aligned}
& \hat{\phi}_v(k)=\int_0^{\infty} \phi(r) J_v(k r) r d r, \\
& \phi(r)=\int_0^{\infty} \hat{\phi}_v(k) J_v(k r) k d k,
\end{aligned}
$$

and denote the Hankel transform/transformation by
$$
{\phi}_v(k)=\mathcal{H}_v\{\phi(r)\}.
$$

It is a fact that
$$
\begin{aligned}
\mathcal{H}_0 & \left\{\partial_r^2 \phi+\frac{1}{r} \partial_r \phi\right\}=-k^2 \hat{\phi}_0(k).
\end{aligned}
$$

Therefore, suppressing the notational dependence on $\nu$, we see that
$$
i \partial_t \hat{\phi}=\frac{k^2}{2} \hat{\phi}+\frac{m^2}{2}\mathcal{H}\left\{\frac{\phi}{r^2}\right\}.
$$

Now the problem has been completely diagonalized, and we are in a position to use yet another integrating factor to finish things off. To this end, let
$$
\hat{g}=e^{\frac{i k^2 t}{2}} \hat{\phi}.
$$
We find that
$$
\partial_t \hat{g}=-\frac{i m^2}{2} e^{\frac{i k^2 t}{2}}\mathcal{H}\left\{\frac{\phi}{r^2}\right\}
$$

which upon an integration gives
$$
\hat{g}_{n+1}=\hat{g}_n-\frac{i m^2}{2} \int_{t_n}^{t_{n+1}}e^{\frac{i k^2 t}{2}}\mathcal{H}\left\{\frac{\phi}{r^2}\right\} d t.
$$

Defining the uniform timestep $h:=t_{n+1}-t_n$ and undoing all of the changes of variables
$$
\varphi_{n+1}=e^{i m\int_{t_n}^{t_{n+1}} a(t) d t} \mathcal{H}^{-1}\left\{e^{-\frac{i k^2 h}{2}} \mathcal{H}\{\varphi_n\}-\frac{im^2}{2}e^{i m\int^{t_{n+1}} a(t) d t}
\int_{t_n}^{t_{n+1}} e^{\frac{i k^2 t}{2}} e^{-i m \int^t a(s) d s} \mathcal{H}\left\{\frac{\varphi}{r^2}\right\} d t\right\}
$$
Thus, the problem has been reduced to quadrature and numerical computation of $2+p$ Hankel transformations ($p$ being the order of the numerical integration method) and one inverse Hankel transformation. Assuming we use a second-order in time method, one time step of the $2D$ linear operator costs $2N$ numerical Hankel transformations and $N$ inverse Hankel transformations. This can all be done very quickly because of the quasi fast Hankel transform: https://opg.optica.org/viewmedia.cfm?r=1&rwjcode=ol&uri=ol-1-1-13&html=true

What's the drawback? How do we even make sense of $\mathcal{H}\left\{\frac{\varphi}{r^2}\right\}$? This is really singular if $\varphi$ doesn't vanish quickly enough at the origin. This seems really restrictive, and doesn't seem like the way to go. I can only imagine using this if the initial condition is an odd excited state of high enough order and as we rotate, the mass concentrates in the vortex ring away from the origin. Is there some other way to save this approach?

Note that if we revisit the original operator splitting and consider moving rotational derivatives to the nonlinear part, this doesn't fix things. The resulting nonlinear ODE will have a singular coefficient due to the factor of $r^{-2}$.

### Option 2: Solve the Problem Exactly

The equation
$$
i \partial_t \phi=-\frac{1}{2} \partial_r^2 \phi-\frac{1}{2 r} \partial_r \phi+\frac{m^2}{2 r^2} \phi
$$
can in principle be solved exactly by separation of variables on a large disc of radius $R$ with vanishing Dirichlet boundary conditions. Let
$$
\phi=e^{-i E t} f(r)
$$
so that
$$
r^2 f^{\prime \prime}+r f^{\prime}+\left(2 E r^2-m^2\right) f=0.
$$

Letting $x=\sqrt{2 E r}$ we have Bessel's equation
$$
x^2 f^{\prime \prime}+x f+\left(x^2-m^2\right) f=0
$$
Because of the boundary conditions, the solution to this equation is given by the Bessel function of the first kind of order $m$ along with a quantization condition
$$
J_m(\sqrt{2E}R)=0
$$
determing the energy level $E_{k,m}$ as proportional to the $k^{\rm th}$ positive root of $J_m.$ By superposition, the truncated solution is given by
$$
\psi(r,\theta,t_{n+1})=\sum_{m=-N / 2+1}^{N / 2} e^{i m \left(\theta+\int_{t_n}^{t_{n+1}}a(s)ds\right)}\sum_{k=1}^{m}A_{m,k} e^{-iE_{k,m}t_{n+1}}J_m(\sqrt{2E_{m,k}}r)
$$
where the coefficients $A_{m,k}$ are determined by orthogonality and the wavefunction $\psi_n$. These coefficients can be computed via Gaussian quadrature.

What's the drawback? This seems like too obvious of an approach, and so I'm left to think that there is something inherently inefficient about computing this many Bessel functions. Of course, I don't have to do this on the fly. I could precompute and store them. I can also parallelize the computation of $A_{m,k}$ as much as possible for each update. What's wrong with this approach?

### Option 3: Direct Numerical Solution
Instead of being clever about it, we can just use a finite element method. See section 3.2 in Bao, et al.: https://scholarsmine.mst.edu/cgi/viewcontent.cgi?article=1154&context=math_stat_facwork
Note that they don't bother with the integrating factor at all, that is, they solve the equation for $\varphi$ and not $\phi$.

Drawback: I would need to implement a lot.

## Preliminary Numerics Done