<figure><img src="image_147.png" /><figcaption>(a) The first and second component of a minimizer of the energy $E_{\varepsilon}^{\Delta}$ using EPG (*).</figcaption></figure>

<figure><img src="image_148.png" /><figcaption>(b) The vortex indices of the first and second component of the minimizer obtained using EPG (*).</figcaption></figure>

<figure><img src="image_149.png" /><figcaption>(c) The first and second component of a minimizer of the energy $E_{\varepsilon}^{\Delta}$ using EPG (**).</figcaption></figure>

<figure><img src="image_150.png" /><figcaption>(d) The vortex indices of the first and second component of the minimizer obtained using EPG (**).</figcaption></figure>

<figure><img src="image_151.png" /><figcaption>(e) The first and second component of a minimizer of the energy $E_{\varepsilon}^{\Delta}$ using GPELab.</figcaption></figure>

<figure><img src="image_152.png" /><figcaption>(f) The vortex indices of the first and second component of the minimizer obtained using GPELab.</figcaption></figure>

<table><thead><tr><th></th><th>E<sub>ε,δ</sub><sup>Δ</sup></th><th>G<sub>0</sub></th><th>K<sup>Δ</sup></th><th>Iteration #</th><th>E<sub>ε,δ</sub><sup>Δ</sup>(u<sub>n</sub>) - E<sub>ε,δ</sub><sup>Δ</sup>(u<sub>n-1</sub>)</th><th>Time(s)</th></tr></thead><tbody><tr><td>GPELab</td><td>-1396</td><td>10<sup>-5</sup></td><td>0.07</td><td>508 797</td><td>-10<sup>-6</sup></td><td>1 526 529</td></tr><tr><td>EPG (**)</td><td>-1389</td><td>10<sup>-5</sup></td><td>0.01</td><td>321 882</td><td>-1.7 × 10<sup>-5</sup></td><td>62 745</td></tr><tr><td>EPG (*)</td><td>-1375</td><td>10<sup>-5</sup></td><td>0.015</td><td>349 265</td><td>-3.7 × 10<sup>-5</sup></td><td>75 224</td></tr></tbody></table>

Figure 3.49 – Comparison between the results of GPELab and EPG for $\Omega = 3$ and $\delta = 1.5$.

We notice that the execution time of GPELab is 20 times greater than EPG and yet the second criterion ($K^{\Delta}$) is 4 times smaller for EPG than for GPELab.

We see more vortices in the first component of the minimizer computed with EPG (*) than in the corresponding component computed with GPELab.

## 5.3 Conclusion

The numerical results displayed in Figures 3.45, 3.46, 3.47, 3.48 and 3.49 show that the methods may converge to different minimizers, with different energies (see Figure 3.45), as well as they may converge to similar minimizers with similar energies (see Figure 3.47), taking into account suitable re-scaling of the discrete energies (as
indicated in Remark 24). The qualitative behaviour of both methods is correct in
the sense that higher rotations end up in creating more vortices, and the beginning
of a hole when Ω = 6.0 (see Figure 3.47) for the one component case. For the two
components case, the qualitative behaviour of both methods is correct in the sense
that they provide coexisting minimizers when δ = 0.7 and segregated minimizers
when δ = 1.5. The EPG method uses more iterations in all cases. Since the EPG
method is explicit and the GPELab method is linearly implicit, the EPG method
is faster in all cases. The order of magnitude of the speed is roughly 3 for all the
simulations except for the two components case with δ = 1.5 where it is roughly 25.
Of course, the stopping criterion is much smaller in minimizers obtained by EPG
than in that computed using GPELab.

In total, there were 7 tests varying from one component BEC with low and
moderate number of points in the discretization to two components BEC in both
segregation and coexistence regimes. These many tests show a better performance
on EPG with roughly 3 times faster computation time and 5 times smaller stopping
criterion ($K^\Delta$) than GPELab.
# Chapter 4

## Theoretical part

In this chapter we adapt a few theorems that are known in the continuous setting to a discrete setting. Using the same notations as before, we prove, first, the existence of a minimizer for the energy defined in (2.1) for only one component with some of its properties. Second we prove a symmetric result related to numerical results that we obtained in Chapter 3.

### 1 Theory for the finite difference discretization

In this section, we consider the two components discrete energy (see equation (2.1)) using a finite difference scheme. We prove in this section some theorems that are known in the continuous setting, as well as symmetric property related to what we observed in some of the simulations.

First of all, we consider the one component discrete energy (with $N_2 = 0$), and we prove the existence of a minimizer for the discrete energy with no rotation in Theorem 29. Then, in Theorem 30, we prove the existence of a real minimizer for the energy without rotation. We proceed to prove in Theorem 31 that in the case of strong segregation, the minimizer for the energy without rotation tends to the confinement function $\rho$. Moreover, we prove also in Theorem 32 that this real minimizer has mass everywhere in the discretization problem, i.e. it is positive. We prove the existence of a global minimizer for the energy with rotation in Theorem 33.

Finally, we consider the two components energy and we prove in Theorem 39, that if the initial datum is symmetric with respect to the several axes, then all the iterations of our gradient method share the same symmetry property.

#### 1.1 Theoretical studies of the energy $E_{\varepsilon}^{\Delta}$

##### Notation used

Let $\tilde{D}$ be the square $[-L, L] \times [-L, L]$ containing the disk $D$ of radius $R < L$ centered at the origin (following the definition used in Chapter 1).

The energy (2.1) is split as follows:

$$ 
\begin{aligned}
E_{\varepsilon}^{\Delta}(\psi) &= (E_{kin})_{\varepsilon}^{\Delta}(\psi) + (E_r)_{\varepsilon}^{\Delta}(\psi) + (E_W)_{\varepsilon}^{\Delta}(\psi) \\
&= F_{\varepsilon}(\psi) + L_{\varepsilon}(\psi),
\end{aligned} 
\quad (4.1) $$
where

* $F_{\varepsilon}(\psi) = (E_{kin})_{\varepsilon}^{\Delta}(\psi) + (E_{W})_{\varepsilon}^{\Delta}(\psi)$, is defined as

$$F_{\varepsilon}(\psi) = \delta_x \delta_y \sum_{n=0}^{N} \sum_{k=0}^{K} \left( \frac{|\psi_{n+1,k} - \psi_{n,k}|^2}{\delta_x^2} + \frac{|\psi_{n,k+1} - \psi_{n,k}|^2}{\delta_y^2} \right) + \frac{\delta_x \delta_y}{4\varepsilon^2} \sum_{n=0}^{N} \sum_{k=0}^{K} (\rho(r_{n,k}) - |\psi_{n,k}|^2)^2,$$

with $\rho$ the non-negative part of the function (3.1) defined in Chapter 3.

* $L_{\varepsilon}(\psi) = (E_r)_{\varepsilon}^{\Delta}(\psi)$, is defined as

$$L_{\varepsilon}(\psi) = -\Omega_{\varepsilon} \delta_x \delta_y \sum_{n=0}^{N} \sum_{k=0}^{K} \Re \left( i \psi_{n,k} [-y_k, x_n] \cdot \left[ \frac{\overline{\psi_{n+1,k} - \psi_{n,k}}}{\delta_x} \frac{\overline{\psi_{n,k+1} - \psi_{n,k}}}{\delta_y} \right] \right).$$

Notice that both $F_{\varepsilon}$ and $L_{\varepsilon}$ are polynomial expressions of $(\psi_{n,k})_{\substack{1 \le n \le N \\ 1 \le k \le K}}$ and $(\overline{\psi_{n,k}})_{\substack{1 \le n \le N \\ 1 \le k \le K}}$ thus continuous over the matrix space $\mathbb{C}^{N \times K}$ with the norm $\|\cdot\|_{\Delta}$ (see (2.2)) defined in Chapter 2.

Our main result in this section is to prove the existence of a minimizer for the energy (4.1) with and without rotation under the following constraints

$$\|\psi\|_{\Delta}^2 = \delta_x \delta_y \sum_{n=1}^{N} \sum_{k=1}^{K} \rho(r_{n,k}) = M. \qquad (4.2)$$

Then, we will prove a symmetric problem we encountered in some of the simulations in Chapter 3 Section 2 (one component condensate with a high rotational speed).

**Definition 27** Let $\mathcal{C}$ be the discrete analogue of the constraints manifold using the norm $\|\cdot\|_{\Delta}$ (see (2.2)), i.e.

$$\mathcal{C} = \{ \psi \in \mathbb{C}^{N \times K} \text{ such that } \|\psi\|_{\Delta}^2 = M \}.$$

**Definition 28** Let the following norm be the discrete analogue of $H_0^1$-norm in the continuous setting

$$\|u\|_{h_0^1}^2 = \delta_x \delta_y \sum_{n=0}^{N} \sum_{k=0}^{K} \left( \frac{|u_{n+1,k} - u_{n,k}|^2}{\delta_x^2} + \frac{|u_{n,k+1} - u_{n,k}|^2}{\delta_y^2} \right),$$

for all $u \in \mathbb{C}^{N \times K}$ with the convention that $u_{n,k} = 0$ as soon as $n = 0$ or $n = N + 1$ or $k = 0$ or $k = K + 1$.

## Theoretical studies of the function $F_{\varepsilon}$

**Theorem 29** For all $\varepsilon > 0$, $F_{\varepsilon}(\psi)$ has a global minimizer on $\mathcal{C}$.
**Proof:** The following function

$$g: \mathbb{C}^{N \times K} \to \mathbb{R}$$
$$U \mapsto \delta_x \delta_y \sum_{n=1}^{N} \sum_{k=1}^{K} |u_{n,k}|^2 - M$$

is continuous over $\mathbb{C}^{N \times K}$ and $\{0\}$ is a closed set in the real normed vector space $\mathbb{R}$. Hence, $g^{-1}(\{0\})$ is a closed set in $\mathbb{C}^{N \times K}$. In addition, the set $g^{-1}(\{0\}) = \mathcal{C}$ is non-empty and bounded. Therefore, $\mathcal{C}$ is compact subset of $\mathbb{C}^{N \times K}$. Since $F_{\varepsilon}$ is continuous over $\mathbb{C}^{N \times K}$, $F_{\varepsilon}$ has a global minimizer over $\mathcal{C}$. $\square$

**Theorem 30** For all $\varepsilon > 0$, $F_{\varepsilon}$ admits a real non-negative minimizer.

**Proof:** Let $\psi^* \in \mathcal{C}$. Using that for all $a, b \in \mathbb{C}$, we have $||a| - |b|| \le |a - b|$, we infer

$$\sum_{n=0}^{N} \sum_{k=0}^{K} \left( \frac{||\psi_{n+1,k} - |\psi_{n,k}|||^2}{\delta_x^2} + \frac{||\psi_{n,k+1} - |\psi_{n,k}|||^2}{\delta_y^2} \right) \\ \le \sum_{n=0}^{N} \sum_{k=0}^{K} \left( \frac{|\psi_{n+1,k} - \psi_{n,k}|^2}{\delta_x^2} + \frac{|\psi_{n,k+1} - \psi_{n,k}|^2}{\delta_y^2} \right).$$

Hence $F_{\varepsilon}(|\psi^*|) \le F_{\varepsilon}(\psi^*)$. Moreover, $|\psi^*| \in \mathcal{C}$. This means there exists a real non-negative minimizer for $F_{\varepsilon}$ in $\mathcal{C}$, using Theorem 29. $\square$

**Theorem 31** Let $\psi^* \in \mathcal{C}$ be a real non-negative minimizer of $F_{\varepsilon}$ over $\mathcal{C}$, then
$$\psi^* \xrightarrow{\|\cdot\|_{\Delta}} \sqrt{\rho} \text{ when } \varepsilon \to 0.$$

**Proof:** For all $\phi \in \mathcal{C}$, $F_{\varepsilon}(\psi^*) \le F_{\varepsilon}(\phi)$. In particular, it is also true for $\phi = \sqrt{\rho}$ since $\|\sqrt{\rho}\|_{\Delta}^2 = M$. Moreover,

$$F_{\varepsilon}(\psi_{\varepsilon}^*) - F_{\varepsilon}(\sqrt{\rho}) = \frac{\delta_x \delta_y}{2} \sum_{n=0}^{N} \sum_{k=0}^{K} \frac{|\psi_{\varepsilon,n+1,k}^* - \psi_{\varepsilon,n,k}^*|^2}{\delta_x^2} + \frac{|\psi_{\varepsilon,n,k+1}^* - \psi_{\varepsilon,n,k}^*|^2}{\delta_y^2} \\ + \frac{\delta_x \delta_y}{4\varepsilon^2} \sum_{n=0}^{N} \sum_{k=0}^{K} (\rho(r_{n,k}) - |\psi_{\varepsilon,n,k}^*|^2)^2 \\ - \frac{\delta_x \delta_y}{2} \sum_{n=0}^{N} \sum_{k=0}^{K} \frac{|\sqrt{\rho(r_{n+1,k})} - \sqrt{\rho(r_{n,k})}|^2}{\delta_x^2} + \frac{|\sqrt{\rho(r_{n,k+1})} - \sqrt{\rho(r_{n,k})}|^2}{\delta_y^2} \\ \le 0,$$
hence,

$$ 
\begin{aligned}
& \frac{\delta_x \delta_y}{4\epsilon^2} \sum_{n=0}^{N} \sum_{k=0}^{K} (\rho(r_{n,k}) - |\psi_{\epsilon,n,k}^*|^2)^2 \\
& \leq -\frac{\delta_x \delta_y}{2} \sum_{n=0}^{N} \sum_{k=0}^{K} \frac{|\psi_{\epsilon,n+1,k}^* - \psi_{\epsilon,n,k}^*|^2}{\delta_x^2} + \frac{|\psi_{\epsilon,n,k+1}^* - \psi_{\epsilon,n,k}^*|^2}{\delta_y^2} \\
& \quad + \frac{\delta_x \delta_y}{2} \sum_{n=0}^{N} \sum_{k=0}^{K} \frac{|\sqrt{\rho(r_{n+1,k})} - \sqrt{\rho(r_{n,k})}|^2}{\delta_x^2} + \frac{|\sqrt{\rho(r_{n,k+1})} - \sqrt{\rho(r_{n,k})}|^2}{\delta_y^2} \\
& \leq \frac{\delta_x \delta_y}{2} \sum_{n=0}^{N} \sum_{k=0}^{K} \frac{|\sqrt{\rho(r_{n+1,k})} - \sqrt{\rho(r_{n,k})}|^2}{\delta_x^2} + \frac{|\sqrt{\rho(r_{n,k+1})} - \sqrt{\rho(r_{n,k})}|^2}{\delta_y^2} \\
& \leq \frac{1}{2} \|\sqrt{\rho}\|_{h_0^1}^2.
\end{aligned}
$$

We conclude with

$$ \|\rho - |\psi_{\epsilon}^*|^2\|_{\Delta}^2 = \delta_x \delta_y \sum_{n=0}^{N} \sum_{k=0}^{K} (\rho(r_{n,k}) - |\psi_{\epsilon,n,k}^*|^2)^2 \leq 2\epsilon^2 \|\sqrt{\rho}\|_{h_0^1}^2. $$

Then $\psi_{\epsilon}^* \xrightarrow{\|\cdot\|_{\Delta}} \sqrt{\rho}$ as $\epsilon \to 0$. □

**Theorem 32** For $\epsilon > 0$ small enough, a non-negative minimizer $\psi^* \in \mathcal{C}$ of $F_\epsilon$ satisfies for all $(n, k) \in \{1, \dots, N\} \times \{1, \dots, N\}$, $\psi_{n,k}^* > 0$.

**Proof:** Let $\psi^* \in \mathcal{C}$ be a non-negative real minimizer of $F_\epsilon$ under the constraints (4.2). Its existence is granted by Theorem 30. We prove the result by contradiction. Let us suppose there exists a couple $(m, l) \in \{1, \dots, N\} \times \{1, \dots, K\}$ such that $\psi_{m_0, l_0}^* = 0$. The set $\mathcal{A} = \{(n, k) \in \{1, \dots, N\} \times \{1, \dots, K\} \text{ such that } \psi_{n,k}^* = 0\}$ contains $(m_0, l_0)$ and has at most $N \times K - 1$ elements since $\psi^* \neq 0$. Therefore, one can choose a couple $(m, l) \in \mathcal{A}$ such that

$$ \psi_{m+1,l}^* > 0 \text{ or } \psi_{m-1,l}^* > 0 \text{ or } \psi_{m,l+1}^* > 0 \text{ or } \psi_{m,l-1}^* > 0. \quad (4.3) $$

Since $\psi^*$ is a minimizer of $F_\epsilon$ over $\mathcal{C}$, it satisfies, for all $v \in \mathbb{C}^{N \times K}$,

$$ \left. \frac{dF_{\epsilon}\left(\frac{\psi^{*}+tv}{\|\psi^{*}+tv\|_{\Delta}}\sqrt{M}\right)}{dt} \right|_{t=0} = 0. $$

Similarly to (2.8) in Chapter 2, we have, for all the $v \in \mathbb{C}^{N \times K}$, the following equation

$$ \left. \frac{dF_{\epsilon}\left(\frac{\psi^{*}+tv}{\|\psi^{*}+tv\|_{\Delta}}\sqrt{M}\right)}{dt} \right|_{t=0} = \langle \nabla F_{\epsilon}(\psi^{*}), v \rangle - \frac{1}{\|\psi^{*}\|_{\Delta}^{2}} \langle \psi^{*}, v \rangle \langle \nabla F_{\epsilon}(\psi^{*}), \psi^{*} \rangle. $$

In particular, for $v = e_{m,l}$ (the canonical base of $\mathbb{C}^{N \times K}$), we have

$$ \langle \nabla F_{\epsilon}(\psi^*) - \frac{1}{\|\psi^*\|_{\Delta}^2} \psi^* \langle \nabla F_{\epsilon}(\psi^*), \psi^* \rangle, e_{m,l} \rangle = \frac{\partial F_{\epsilon}(\psi^*)}{\partial \psi_{m,l}} - \frac{\psi_{m,l}^*}{\|\psi^*\|_{\Delta}^2} \langle \nabla F_{\epsilon}(\psi^*), \psi^* \rangle = 0. $$
Since $\psi_{m,l}^* = 0$, we have $\frac{\partial F_\epsilon(\psi^*)}{\partial \psi_{m,l}} = 0$. According to Section 1.2 in Chapter 2, we have

$$
\begin{align*}
\frac{\partial F_{\epsilon}(\psi^*)}{\partial \psi_{m,l}} &= \frac{\delta_y}{\delta_x}(2\psi_{m,l}^* - \psi_{m-1,l}^* - \psi_{m+1,l}^*) + \frac{\delta_x}{\delta_y}(2\psi_{m,l}^* - \psi_{m,l-1}^* - \psi_{m,l+1}^*) \\
&\quad - \frac{\delta_x \delta_y}{\epsilon^2}\psi_{m,l}^*(\rho(r_{n,k}) - \psi_{n,k}^{*2}) \\
&= \frac{\delta_y}{\delta_x}(-\psi_{m-1,l}^* - \psi_{m+1,l}^*) + \frac{\delta_x}{\delta_y}(-\psi_{m,l-1}^* - \psi_{m,l+1}^*) \\
&= 0.
\end{align*}
$$

Using the non-negativity of $\psi^*$ we infer $\psi_{m+1,l}^* = \psi_{m-1,l}^* = \psi_{m,l+1}^* = \psi_{m,l-1}^* = 0$. This is in contradiction with (4.3). □

## Theoretical studies of the function $E_\epsilon^\Delta$

**Theorem 33** For all $\epsilon > 0$, $E_\epsilon^\Delta(\psi)$ has a global minimizer in $\mathcal{C}$.

*Proof:* Similarly to the proof of Theorem 30, $\mathcal{C}$ is compact. Since $E_\epsilon^\Delta$ is continuous over $\mathbb{C}^{N \times K}$, then $E_\epsilon^\Delta$ has a global minimizer over $\mathcal{C}$. □

## 1.2 Symmetric theory

In this Section we study a particular symmetric property related to the gradient method used in EPG (see Section 4 of Chapter 2). Therefore $N$ is equal to $K$.

The gradient of the continuous energy (1.13) conserves any orthogonal symmetry with respect to any line that crosses the origin in the following sense. Let $S$ be an orthogonal symmetry of $\mathbb{R}^2$, and $u_1, u_2 \in H_0^1(D, \mathbb{C})$, we have $\overline{\nabla (E_{\epsilon,\delta}^\Omega(u_1, u_2))} \circ S = \nabla E_{\epsilon,\delta}^\Omega(\overline{u_1 \circ S}, \overline{u_2 \circ S})$. Therefore if $u_1, u_2$ satisfy $\overline{u_\ell \circ S} = u_\ell$ for all $\ell = 1, 2$, then we have

$$ \overline{\nabla (E_{\epsilon,\delta}^\Omega(u_1, u_2)) \circ S} = \nabla E_{\epsilon,\delta}^\Omega(u_1, u_2), $$

which gives exactly the preservation of the orthogonal symmetry for the gradient of the energy. In the following of this section, we investigate the preservation of this property with respect to four axes (x and y-axes, first and second diagonals) for the numerical approaches. We start with a proof of the aforementioned result.

**Proposition 34** Let $S$ be an orthogonal symmetry of $\mathbb{R}^2$, and $u_1, u_2 \in H_0^1(D, \mathbb{C})$. If $u_1, u_2$ satisfy $\overline{u_\ell \circ S} = u_\ell$ for all $\ell = 1, 2$ then we have

$$ \overline{\nabla (E_{\epsilon,\delta}^\Omega(u_1, u_2)) \circ S} = \nabla E_{\epsilon,\delta}^\Omega(u_1, u_2), $$

with $E_{\epsilon,\delta}^\Omega$ being the continuous energy defined in (1.13).

*Proof:* Let $S$ be an orthogonal symmetry of $\mathbb{R}^2$ defined as follows

$$ S = \mathcal{R}^t \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \mathcal{R}, $$
where $\mathcal{R} = \begin{pmatrix} \cos(\theta) & \sin(\theta) \\ \sin(\theta) & -\cos(\theta) \end{pmatrix}$ and $\theta \in [-\pi, \pi]$ is the rotation angle. Therefore, we can rewrite $S$ as follows

$$S = \begin{pmatrix} \cos^2(\theta) - \sin^2(\theta) & 2 \cos(\theta) \sin(\theta) \\ 2 \cos(\theta) \sin(\theta) & -\cos^2(\theta) + \sin^2(\theta) \end{pmatrix} = \begin{pmatrix} a & b \\ b & -a \end{pmatrix},$$

where $a = \cos^2(\theta) - \sin^2(\theta)$ and $b = 2 \cos(\theta) \sin(\theta)$.

According to the hypothesis, we have $\overline{u_{\ell} \circ S} = u_{\ell}$ for $\ell = 1, 2$. Our goal is to prove that this symmetry holds for the gradient of the energy under this hypothesis.

## Gradient of the kinetic energy:

We have for the kinetic gradient and for $\ell = 1, 2$

$$\begin{aligned}
\nabla E_{k,\varepsilon,\delta}^{\Omega}(\overline{u_{\ell} \circ S}) &= -\overline{\Delta(u_{\ell} \circ S)} \\
&= -(a^2 + b^2)\overline{\partial_x^2 u_{\ell}} \circ S - 2(ab - ab)\overline{\partial_x \partial_y u_{\ell}} \circ S - (a^2 + b^2)\overline{\partial_y^2 u_{\ell}} \circ S \\
&= -((\cos^2(\theta) - \sin^2(\theta))^2 + (2 \cos(\theta) \sin(\theta))^2)\overline{((\partial_x^2 + \partial_y^2)(u_{\ell}) \circ S)} \\
&= -(\cos^2(\theta) + \sin^2(\theta))\overline{((\partial_x^2 + \partial_y^2)(u_{\ell}) \circ S)} \\
&= -\overline{((\partial_x^2 + \partial_y^2)(u_{\ell}) \circ S)} \\
&= -\overline{\Delta u_{\ell} \circ S} \\
&= \overline{\nabla E_{k,\varepsilon,\delta}^{\Omega}(u_{\ell}) \circ S}
\end{aligned}$$

The symmetry property holds for the gradient of the kinetic energy.

## Gradient of the rotational energy:

Note that $\nabla E_{r,\varepsilon,\delta}^{\Omega}(u_{\ell}) = -2\Omega i \mathbf{x}^{\perp} \cdot \nabla u_{\ell}$ with $\mathbf{x}^{\perp} = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} \mathbf{x}$. Let $P = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}$, then $SP = -PS$, and $S^t = S$ since $S$ is symmetric. We have for the rotational gradient and for $\ell = 1, 2$

$$\begin{aligned}
\nabla E_{r,\varepsilon,\delta}^{\Omega}(\overline{u_{\ell} \circ S}) &= -2\Omega i \mathbf{x}^{\perp} \cdot \overline{\nabla(u_{\ell} \circ S)} \\
&= -2\Omega i \mathbf{x}^{\perp} \cdot S(\overline{\nabla u_{\ell} \circ S}) \\
&= -2\Omega i P \begin{pmatrix} x \\ y \end{pmatrix} \cdot S(\overline{\nabla u_{\ell} \circ S}) \\
&= -2\Omega i S^t P \begin{pmatrix} x \\ y \end{pmatrix} \cdot (\overline{\nabla u_{\ell} \circ S}) \\
&= 2\Omega i P S \begin{pmatrix} x \\ y \end{pmatrix} \cdot (\overline{\nabla u_{\ell} \circ S}) \\
&= 2\Omega i \left( S \begin{pmatrix} x \\ y \end{pmatrix} \right)^{\perp} \cdot (\overline{\nabla u_{\ell} \circ S}) \\
&= \overline{-2\Omega i \left( S \begin{pmatrix} x \\ y \end{pmatrix} \right)^{\perp} \cdot (\overline{\nabla u_{\ell} \circ S})} \\
&= \overline{\nabla E_{r,\varepsilon,\delta}^{\Omega}(u_{\ell}) \circ S}
\end{aligned}$$

The symmetry holds for the gradient of the rotational energy.
**Gradient of the confinement and interaction energy:**

Note that, because $\overline{u_{\ell} \circ S} = u_{\ell}$, $|\cdot|^2$ is invariant under rotations of the plane around the origin, i.e. $|u \circ S|^2 = |u|^2$. We have for the confinement and interaction gradient and for $\ell = 1, 2$

$$ 
\begin{align*}
\nabla E_{W,\varepsilon,\delta}^{\Omega}(\overline{u_{\ell} \circ S}) &= -\frac{1}{\varepsilon^2} \overline{u_{\ell} \circ S}(\rho(r) - |\overline{u_{\ell} \circ S}|^2 - |\overline{u_{3-\ell} \circ S}|^2) + \frac{\delta-1}{\varepsilon^2} \overline{u_{\ell} \circ S} |\overline{u_{3-\ell} \circ S}|^2 \\
&= -\frac{1}{\varepsilon^2} \overline{u_{\ell} \circ S}(\rho(r) - |u_{\ell}|^2 - |u_{3-\ell}|^2) + \frac{\delta-1}{\varepsilon^2} \overline{u_{\ell} \circ S} |u_{3-\ell}|^2 \\
&= \left( -\frac{1}{\varepsilon^2}(\rho(r) - |u_{\ell}|^2 - |u_{3-\ell}|^2) + \frac{\delta-1}{\varepsilon^2} |u_{3-\ell}|^2 \right) \overline{u_{\ell} \circ S} \\
&= \left( -\frac{1}{\varepsilon^2}(\rho(r) - |u_{\ell}|^2 - |u_{3-\ell}|^2) + \frac{\delta-1}{\varepsilon^2} |u_{3-\ell}|^2 \right) \overline{u_{\ell} \circ S} \\
&= \overline{\nabla E_{W,\varepsilon,\delta}^{\Omega}(u_{\ell}) \circ S}
\end{align*}
$$

The symmetry holds for the gradient of the confinement and interaction energy.

Putting together all the previous results, we obtain $\nabla E_{\varepsilon,\delta}^{\Omega}(u_{\ell}) = \nabla E_{\varepsilon,\delta}^{\Omega}(\overline{u_{\ell} \circ S}) = \overline{\nabla E_{\varepsilon,\delta}^{\Omega}(u_{\ell}) \circ S}$ under the hypothesis that $\overline{u_{\ell} \circ S} = u_{\ell}$ for $\ell = 1, 2$, thus proving the proposition. □

In the following of this section, we investigate the preservation of this property with respect to four axes (x and y-axes, first and second diagonals) for the numerical approaches.

**Definition 35** Let $\psi \in \mathbb{C}^{N \times N}$. We say $\psi$ is symmetric (resp. skew-symmetric) with respect to first diagonal if $\psi_{n,k} = \psi_{k,n}$ (resp. $\psi_{n,k} = \bar{\psi}_{k,n}$) for all $n, k \in \{1, \dots, N\}$.

**Definition 36** Let $\psi \in \mathbb{C}^{N \times N}$. We say $\psi$ is symmetric (resp. skew-symmetric) with respect to second diagonal if $\psi_{n,k} = \psi_{N+1-k,N+1-n}$ (resp. $\psi_{n,k} = \bar{\psi}_{N+1-k,N+1-n}$) for all $n, k \in \{1, \dots, N\}$.

**Remark 37** For Proposition 38 and for Theorems 39 and 41, we suppose that the given confinement function $\rho$ is symmetric with respect to the first and second diagonal, x-axis and y-axis. This is the case if the function $\rho$ depends only on the distance from the origin as for example the function (3.1).

**Proposition 38** Let $\psi^1, \psi^2 \in \mathbb{C}^{N \times N}$ and $\rho$ satisfying the hypothesis of Remark 37. If $\psi^1$ and $\psi^2$ share the same symmetry property (symmetric or skew-symmetric with respect to the first or second diagonal, x-axis or y-axis.) then the gradient of the discrete confinement energy shares the same symmetry property as $\psi^l$ for $l = 1, 2$.

**Proof:** Assume $\psi^l$ share both the same symmetry property for $l = 1, 2$. If both of them are symmetric or skew-symmetric with respect to $d^{th}$ diagonal (for $d = 1, 2$) then $|\psi^l|^2$ is symmetric with respect to the $d^{th}$ diagonal. According to Remark 37, the function $\rho$ is already symmetric with respect to the first and second diagonal.

**First diagonal case:**
Assume $\psi^\ell$ is symmetric (resp. skew-symmetric) with respect to first diagonal for $\ell = 1, 2$, then for all $n, k \in \{1, \dots, N\}$, according to Section 1.2, we have

$$ \begin{aligned} \frac{\partial(E_W)^{\Delta}_{\varepsilon, \delta}}{\partial p_{k,n}^{\ell}}(\psi^1, \psi^2) &= -\frac{\delta_x^2}{\varepsilon^2} p_{k,n}^{\ell} \left(\rho_{k,n} - |\psi_{k,n}^1|^2 - |\psi_{k,n}^2|^2\right) + \frac{\delta_x^2(\delta-1)}{\varepsilon^2} p_{k,n}^{\ell} |\psi_{k,n}^{3-\ell}|^2 \\ &= -\frac{\delta_x^2}{\varepsilon^2} p_{n,k}^{\ell} \left(\rho_{n,k} - |\psi_{n,k}^1|^2 - |\psi_{n,k}^2|^2\right) + \frac{\delta_x^2(\delta-1)}{\varepsilon^2} p_{k,n}^{\ell} |\psi_{n,k}^{3-\ell}|^2 \\ &= \frac{\partial(E_W)^{\Delta}_{\varepsilon, \delta}}{\partial p_{n,k}^{\ell}}(\psi^1, \psi^2), \end{aligned} $$

$$ \begin{aligned} \frac{\partial(E_W)^{\Delta}_{\varepsilon, \delta}}{\partial q_{k,n}^{\ell}}(\psi^1, \psi^2) &= -\frac{\delta_x^2}{\varepsilon^2} q_{k,n}^{\ell} \left(\rho_{k,n} - |\psi_{k,n}^1|^2 - |\psi_{k,n}^2|^2\right) + \frac{\delta_x^2(\delta-1)}{\varepsilon^2} q_{k,n}^{\ell} |\psi_{k,n}^{3-\ell}|^2 \\ &= -\frac{\delta_x^2}{\varepsilon^2} q_{k,n}^{\ell} \left(\rho_{n,k} - |\psi_{n,k}^1|^2 - |\psi_{n,k}^2|^2\right) + \frac{\delta_x^2(\delta-1)}{\varepsilon^2} q_{k,n}^{\ell} |\psi_{n,k}^{3-\ell}|^2 \\ &= (\pm 1) \frac{\partial(E_W)^{\Delta}_{\varepsilon, \delta}}{\partial q_{n,k}^{\ell}}(\psi^1, \psi^2), \end{aligned} $$

where the last sign is +1 if $\psi^\ell$ is symmetric and -1 if $\psi^\ell$ is skew-symmetric.

## Second diagonal case:

Assume $\psi^\ell$ is symmetric (resp. skew-symmetric) with respect to second diagonal, then we have, for all $n, k \in \{1, \dots, N\}$

$$ \begin{aligned} \frac{\partial(E_W)^{\Delta}_{\varepsilon, \delta}}{\partial p_{N+1-k,N+1-n}^{\ell}}(\psi^1, \psi^2) &= -\frac{\delta_x^2}{\varepsilon^2} p_{N+1-k,N+1-n}^{\ell} \left(\rho_{N+1-k,N+1-n} - \sum_{j=1}^2 |\psi_{N+1-k,N+1-n}^j|^2\right) \\ &\quad + \frac{\delta_x^2(\delta-1)}{\varepsilon^2} p_{N+1-k,N+1-n}^{\ell} |\psi_{N+1-k,N+1-n}^{3-\ell}|^2 \\ &= -\frac{\delta_x^2}{\varepsilon^2} p_{n,k}^{\ell} \left(\rho_{n,k} - |\psi_{n,k}^1|^2 - |\psi_{n,k}^2|^2\right) + \frac{\delta_x^2(\delta-1)}{\varepsilon^2} p_{n,k}^{\ell} |\psi_{n,k}^{3-\ell}|^2 \\ &= \frac{\partial(E_W)^{\Delta}_{\varepsilon, \delta}}{\partial p_{n,k}^{\ell}}(\psi^1, \psi^2), \\ \frac{\partial(E_W)^{\Delta}_{\varepsilon, \delta}}{\partial q_{N+1-k,N+1-n}^{\ell}}(\psi^1, \psi^2) &= -\frac{\delta_x^2}{\varepsilon^2} q_{N+1-k,N+1-n}^{\ell} \left(\rho_{N+1-k,N+1-n} - \sum_{j=1}^2 |\psi_{N+1-k,N+1-n}^j|^2\right) \\ &\quad + \frac{\delta_x^2(\delta-1)}{\varepsilon^2} q_{N+1-k,N+1-n}^{\ell} |\psi_{N+1-k,N+1-n}^{3-\ell}|^2 \\ &= q_{N+1-k,N+1-n}^{\ell} \left(-\frac{\delta_x^2}{\varepsilon^2} \left(\rho_{n,k} - |\psi_{n,k}^1|^2 - |\psi_{n,k}^2|^2\right) + \frac{\delta_x^2(\delta-1)}{\varepsilon^2} |\psi_{n,k}^{3-\ell}|^2\right) \\ &= (\pm 1) \frac{\partial(E_W)^{\Delta}_{\varepsilon, \delta}}{\partial q_{n,k}^{\ell}}(\psi^1, \psi^2), \end{aligned} $$

where the last sign is +1 if $\psi^\ell$ is symmetric and -1 if $\psi^\ell$ is skew-symmetric.

## x-axis case:
Assume $\psi^\ell$ is symmetric (resp. skew-symmetric) with respect to the x-axis, then we have, for all $n, k \in \{1, \dots, N\}$

$$ 
\begin{aligned}
\frac{\partial(E_W)^\Delta_{\varepsilon,\delta}}{\partial p_{N+1-n,k}^\ell}(\psi^1, \psi^2) &= -\frac{\delta_x^2}{\varepsilon^2} p_{N+1-n,k}^\ell \left( \rho_{N+1-n,k} - \left|\psi_{N+1-n,k}^1\right|^2 - \left|\psi_{N+1-n,k}^2\right|^2 \right) \\
&\quad + \frac{\delta_x^2(\delta-1)}{\varepsilon^2} p_{N+1-n,k}^\ell \left|\psi_{N+1-n,k}^{3-\ell}\right|^2 \\
&= -\frac{\delta_x^2}{\varepsilon^2} p_{n,k}^\ell \left( \rho_{n,k} - \left|\psi_{n,k}^1\right|^2 - \left|\psi_{n,k}^2\right|^2 \right) + \frac{\delta_x^2(\delta-1)}{\varepsilon^2} p_{n,k}^\ell \left|\psi_{n,k}^{3-\ell}\right|^2 \\
&= \frac{\partial(E_W)^\Delta_{\varepsilon,\delta}}{\partial p_{n,k}^\ell}(\psi^1, \psi^2),
\end{aligned}
$$

$$ 
\begin{aligned}
\frac{\partial(E_W)^\Delta_{\varepsilon,\delta}}{\partial q_{N+1-n,k}^\ell}(\psi^1, \psi^2) &= -\frac{\delta_x^2}{\varepsilon^2} q_{N+1-n,k}^\ell \left( \rho_{N+1-n,k} - \left|\psi_{N+1-n,k}^1\right|^2 - \left|\psi_{N+1-n,k}^2\right|^2 \right) \\
&\quad + \frac{\delta_x^2(\delta-1)}{\varepsilon^2} q_{N+1-n,k}^\ell \left|\psi_{N+1-n,k}^{3-\ell}\right|^2 \\
&= -\frac{\delta_x^2}{\varepsilon^2} q_{N+1-n,k}^\ell \left( \rho_{n,k} - \sum_{j=1}^2 \left|\psi_{n,k}^j\right|^2 \right) + \frac{\delta_x^2(\delta-1)}{\varepsilon^2} q_{N+1-n,k}^\ell \left|\psi_{n,k}^{3-\ell}\right|^2 \\
&= (\pm 1) \frac{\partial(E_W)^\Delta_{\varepsilon,\delta}}{\partial q_{n,k}^\ell}(\psi^1, \psi^2),
\end{aligned}
$$

where the last sign is +1 if $\psi^\ell$ is symmetric and -1 if $\psi^\ell$ is skew-symmetric.

**y-axis case:**

Similarly to the x-axis case, we can prove that

$$ 
\begin{aligned}
\frac{\partial(E_W)^\Delta_{\varepsilon,\delta}}{\partial p_{n,N+1-k}^\ell}(\psi^1, \psi^2) &= \frac{\partial(E_W)^\Delta_{\varepsilon,\delta}}{\partial p_{n,k}^\ell}(\psi^1, \psi^2), \\
\frac{\partial(E_W)^\Delta_{\varepsilon,\delta}}{\partial q_{n,N+1-k}^\ell}(\psi^1, \psi^2) &= (\pm 1) \frac{\partial(E_W)^\Delta_{\varepsilon,\delta}}{\partial q_{n,k}^\ell}(\psi^1, \psi^2),
\end{aligned}
$$

where the last sign is +1 if $\psi^\ell$ is symmetric and -1 if $\psi^\ell$ is skew-symmetric.

This proves that the gradient of the discrete confinement energy shares the same symmetry property as $\psi^\ell$ for $\ell = 1, 2$. □

**Theorem 39 (Finite difference scheme)** If the initial datum of the gradient descent algorithm is such that the two components are either skew-symmetric with respect to the first or skew-symmetric with respect to the second diagonal or skew-symmetric with respect to the x-axis or skew-symmetric with respect to the y-axis and $\rho$ satisfies the hypothesis of Remark 37, then all of the iterations of the method share the same symmetry property.

**Proof:** It is sufficient to prove the result for one iteration. Let $\psi_{n,k}^\ell$ for $\ell \in \{1, 2\}$ be the initial data of one iteration and let us denote by $\phi_{n,k}^\ell$ for $\ell = 1, 2$ the result of the iteration. For the convenience of the reader, we introduce $g_1, g_2, g_3$ and $g_4$ as functions from $\{1, \dots, N\}^2$ to itself defined as

$$ 
\begin{aligned}
g_1(n, k) &= (k, n), \\
g_2(n, k) &= (N + 1 - k, N + 1 - n), \\
g_3(n, k) &= (N + 1 - n, k), \\
g_4(n, k) &= (n, N + 1 - k).
\end{aligned}
$$
Notice that $g_1$ correspond to the first diagonal, $g_2$ corresponds to the second diagonal, $g_3$ corresponds to the y-axis and $g_4$ corresponds to the x-axis. Depending on the skew-symmetry property we have $\psi_{n,k}^\ell = \psi_{g_j(n,k)}^\ell$ for all $n, k \in \{1, \dots, N\}$ and in particular

$$p_{n,k}^{\ell} = p_{g_j(n,k)}^{\ell} \quad \text{and} \quad q_{n,k}^{\ell} = -q_{g_j(n,k)}^{\ell}. \qquad (4.4)$$

According to Sections 1.2 and 4 of Chapter 2, we have

$$\phi^{\ell} = \psi^{\ell} - h \nabla E_{\varepsilon, \delta}^{\Delta}(\psi^{\ell}),$$

with

$$\left(\nabla E_{\varepsilon, \delta}^{\Delta}(\psi^{\ell})\right)_{n,k} = \left( \frac{\partial(E_{kin})_{\varepsilon}^{\Delta}}{\partial p_{n,k}^{\ell}}(\psi^{\ell}) + \frac{\partial(E_r)_{\varepsilon}^{\Delta}}{\partial p_{n,k}^{\ell}}(\psi^{\ell}) + \frac{\partial(E_W)_{\varepsilon, \delta}^{\Delta}}{\partial p_{n,k}^{\ell}}(\psi^1, \psi^2) \right. \\ \left. \frac{\partial(E_{kin})_{\varepsilon}^{\Delta}}{\partial q_{n,k}^{\ell}}(\psi^{\ell}) + \frac{\partial(E_r)_{\varepsilon}^{\Delta}}{\partial q_{n,k}^{\ell}}(\psi^{\ell}) + \frac{\partial(E_W)_{\varepsilon, \delta}^{\Delta}}{\partial q_{n,k}^{\ell}}(\psi^1, \psi^2) \right).$$

It is sufficient to prove that each one of these partial derivatives shares the same symmetry property as $\psi^\ell$, i.e. for all $n, k \in \{1, \dots, N\}$ and for $\ell = 1, 2$, we have to prove one of the corresponding property below

$$\frac{\partial E_{\varepsilon, \delta}^{\Delta}(\psi^{\ell})}{\partial p_{g_j(n,k)}} = \frac{\partial E_{\varepsilon, \delta}^{\Delta}(\psi^{\ell})}{\partial p_{n,k}} \quad \text{and} \quad \frac{\partial E_{\varepsilon, \delta}^{\Delta}(\psi^{\ell})}{\partial q_{g_j(n,k)}} = -\frac{\partial E_{\varepsilon, \delta}^{\Delta}(\psi^{\ell})}{\partial q_{n,k}}.$$

## Partial derivative of the discrete kinetic energy:

Using (4.4) and Section 1.2 of Chapter 2, we have

$$\begin{aligned} \frac{\partial(E_{kin})_{\varepsilon}^{\Delta}(\psi^{\ell})}{\partial p_{g_j(n,k)}^{\ell}} &= 4p_{g_j(n,k)}^{\ell} - p_{g_j(n,k)+(1,0)}^{\ell} - p_{g_j(n,k)-(1,0)}^{\ell} - p_{g_j(n,k)+(0,1)}^{\ell} - p_{g_j(n,k)-(0,1)}^{\ell} \\ &= 4p_{n,k}^{\ell} - p_{n-1,k}^{\ell} - p_{n+1,k}^{\ell} - p_{n,k-1}^{\ell} - p_{n,k+1}^{\ell} \\ &= \frac{\partial(E_{kin})_{\varepsilon}^{\Delta}(\psi^{\ell})}{\partial p_{n,k}^{\ell}}. \end{aligned}$$

In the same way we can obtain

$$\frac{\partial(E_{kin})_{\varepsilon}^{\Delta}(\psi^{\ell})}{\partial q_{g_j(n,k)}^{\ell}} = -\frac{\partial(E_{kin})_{\varepsilon}^{\Delta}(\psi^{\ell})}{\partial q_{n,k}^{\ell}}.$$

## Partial derivative of the discrete rotational energy:

For the partial derivative of the discrete rotational energy we have to distinguish two cases. The first one corresponds to $g_1$, $g_2$ for which the proof relies on switching the roles of $x$ and $y$. The second one corresponds to $g_3$, $g_4$ which share some similarities.

Let us first consider the partial derivative with respect to the real part. The proof being the same for $g_1$ we only consider the case of $g_2$. Using (4.4) and the fact that for all $n, k \in \{1, \dots, N\}$,

$$\begin{aligned} y_{N+1-k} &= -L + (N + 1 - k)\delta_y = -L + 2L - k\delta_y = -y_k, \\ x_{N+1-n} &= -L + (N + 1 - n)\delta_x = -L + 2L - n\delta_x = -x_n, \\ x_n &= y_n, \end{aligned}$$
we have

$$
\begin{align*}
\frac{\partial(E_r)_\varepsilon^\Delta(\psi^\ell)}{\partial p_{g_2(n,k)}^\ell} &= -\Omega_\varepsilon \delta_x \left( y_{N+1-n}(q_{g_2(n,k)}^\ell-(1,0)) - q_{g_2(n,k)}^\ell+(1,0)) \right. \\
&\qquad \left. + x_{N+1-k}(q_{g_2(n,k)}^\ell+(0,1)) - q_{g_2(n,k)}^\ell-(0,1)) \right) \\
&= -\Omega_\varepsilon \delta_x \left( -y_n(-q_{n,k+1}^\ell + q_{n,k-1}^\ell) - x_k(-q_{n-1,k}^\ell + q_{n+1,k}^\ell) \right) \\
&= -\Omega_\varepsilon \delta_x \left( x_n(q_{n,k+1}^\ell - q_{n,k-1}^\ell) + y_k(q_{n-1,k}^\ell - q_{n+1,k}^\ell) \right) \\
&= \frac{\partial(E_r)_\varepsilon^\Delta(\psi^\ell)}{\partial p_{n,k}^\ell}.
\end{align*}
$$

This proves the result for $g_1$ and $g_2$. Let us now consider only $g_3$, since the proof is similar for $g_4$. We have

$$
\begin{align*}
\frac{\partial(E_r)_\varepsilon^\Delta(\psi^\ell)}{\partial p_{g_3(n,k)}^\ell} &= -\Omega_\varepsilon \delta_x \left( y_k(q_{g_3(n,k)}^\ell-(1,0)) - q_{g_3(n,k)}^\ell+(1,0)) \right. \\
&\qquad \left. + x_{N+1-n}(q_{g_3(n,k)}^\ell+(0,1)) - q_{g_3(n,k)}^\ell-(0,1)) \right) \\
&= -\Omega_\varepsilon \delta_x \left( y_k(-q_{n+1,k}^\ell + q_{n-1,k}^\ell) - x_n(-q_{n,k+1}^\ell + q_{n,k-1}^\ell) \right) \\
&= -\Omega_\varepsilon \delta_x \left( y_k(q_{n-1,k}^\ell - q_{n+1,k}^\ell) + x_n(q_{n,k+1}^\ell - q_{n,k-1}^\ell) \right) \\
&= \frac{\partial(E_r)_\varepsilon^\Delta(\psi^\ell)}{\partial p_{n,k}^\ell}.
\end{align*}
$$

In the same way we can prove that for $j = 1, 2, 3, 4$, we have

$$
\frac{\partial(E_r)_\varepsilon^\Delta(\psi^\ell)}{\partial q_{g_j(n,k)}^\ell} = -\frac{\partial(E_r)_\varepsilon^\Delta(\psi^\ell)}{\partial q_{n,k}^\ell}.
$$

Partial derivative of the discrete confinement energy:

Last, using Proposition 38 and under the assumption $\psi_{n,k}^\ell = \bar{\psi}_{g_j(n,k)}^\ell$, we can
prove that

$$
\begin{align*}
\frac{\partial(E_W)_\varepsilon^\Delta(\psi^1, \psi^2)}{\partial p_{g_j(n,k)}^\ell} &= \frac{\partial(E_W)_\varepsilon^\Delta(\psi^1, \psi^2)}{\partial p_{n,k}^\ell}, \\
\frac{\partial(E_W)_\varepsilon^\Delta(\psi^1, \psi^2)}{\partial q_{g_j(n,k)}^\ell} &= -\frac{\partial(E_W)_\varepsilon^\Delta(\psi^1, \psi^2)}{\partial q_{n,k}^\ell}.
\end{align*}
$$

Putting together all the previous results we obtain $\psi_{n,k}^{\ell} = \overline{\psi_{g_j(n,k)}^{\ell}}$, for the corre-
sponding $j = 1, 2, 3, 4$ (depending on the hypothesis) and for all $n, k \in \{1, \dots, N\}$
and $\ell = 1, 2$. This completes the proof of Theorem 39.
$\square$

**Example 40** In this example, we display the results of the gradient descent method (see Section 4 of Chapter 2) with one component $\psi$ using 19 points in the discretization after 1 iteration. We consider a real initial datum depending only on the distance from the center hence satisfying all the conditions of Theorem 39 (similarly to Remark 37). As we can see in Figures 4.1c and 4.1a, the squared modulus of $\psi$ and its real part both conserve the symmetry with respect to all axes, while in Figure 4.1b we can see that the imaginary part of $\psi$ conserves the skew-symmetry property with respect to all axes.
<figure><img src="image_153.png" /><figcaption>Figure 4.1 – Example of symmetry conservation using the finite difference approach.</figcaption></figure>

# 2 Theory for the Fast Fourier transform discretization

In this section, we consider the two components discrete energy (see Equation (2.4)) using fast Fourier discretization. We prove in this section a symmetry property related to some of the simulations of Chapter 3. Moreover we prove a Plancherel’s equality.

First, we prove in Theorem 41, that if the initial datum is skew-symmetric with respect to the second diagonal i.e. $\overline{\psi_{N+1-k,N+1-n}^\ell} = \psi_{n,k}^\ell$, then all of our gradient method’s iterations share the same skew-symmetry property. Then, in Theorem 44, we prove a Plancherel’s equality.

**Theorem 41 (Fast Fourier transform scheme)** If the initial data’s matrices of the gradient descent algorithm are skew-symmetric with respect to the second diagonal (see definition (36)) and $\rho$ satisfies the hypothesis of Remark 37, then all the iterations of the EPG algorithm 4 share the same skew-symmetry property.

**Proof** : It is sufficient to prove the result for one iteration. Let $\psi_{n,k}^\ell$ for $\ell \in \{1, 2\}$ be the initial data of one iteration and let us denote by $\phi_{n,k}^\ell$ for $\ell = 1, 2$ the result of the iteration. According to Sections 2.3 and 4 of Chapter 2, we have the following relation

$$ \phi^\ell = \psi^\ell - h \nabla E_{\varepsilon, \delta}^{\Delta}(\psi^\ell), $$

with

$$ \left(\nabla E_{\varepsilon, \delta}^{\Delta}(\psi^{\ell})\right)_{n,k} = \left( \frac{\partial(E_{kin})_{\varepsilon}^{\Delta}}{\partial p_{n,k}^{\ell}}(\psi^{\ell}) + \frac{\partial(E_{r})_{\varepsilon}^{\Delta}}{\partial p_{n,k}^{\ell}}(\psi^{\ell}) + \frac{\partial(E_{W})_{\varepsilon, \delta}^{\Delta}}{\partial p_{n,k}^{\ell}}(\psi^{1}, \psi^{2}) \right) \left( \frac{\partial(E_{kin})_{\varepsilon}^{\Delta}}{\partial q_{n,k}^{\ell}}(\psi^{\ell}) + \frac{\partial(E_{r})_{\varepsilon}^{\Delta}}{\partial q_{n,k}^{\ell}}(\psi^{\ell}) + \frac{\partial(E_{W})_{\varepsilon, \delta}^{\Delta}}{\partial q_{n,k}^{\ell}}(\psi^{1}, \psi^{2}) \right). $$

It is sufficient to prove that each one of these partial derivatives shares the same skew-symmetry property as $\psi^\ell$, i.e. for all $n, k \in \{1, \dots, N\}$ and for $\ell = 1, 2$, we have to prove

$$ \frac{\partial E_{\varepsilon, \delta}^{\Delta}(\psi^{\ell})}{\partial p_{N+1-k, N+1-n}} = \frac{\partial E_{\varepsilon, \delta}^{\Delta}(\psi^{\ell})}{\partial p_{n,k}} \quad \text{and} \quad \frac{\partial E_{\varepsilon, \delta}^{\Delta}(\psi^{\ell})}{\partial q_{N+1-k, N+1-n}} = -\frac{\partial E_{\varepsilon, \delta}^{\Delta}(\psi^{\ell})}{\partial q_{n,k}}. $$
Partial derivative of the discrete kinetic energy:

Using the fact that $\xi_n = \lambda_n$ for all $n \in \{0, \dots, N+1\}$ (see (2.3)) and the definition of $g_2$ in the proof of Theorem 39, we have for the kinetic partial derivative:

$$
\begin{align*}
& \frac{\partial E_{kin}^{\Delta}(\psi^{\ell})}{\partial p_{N+1-j,N+1-i}^{\ell}} \\
&= \delta_x^2 \Re \left( \left[ \mathbf{iff}_x (\Xi^2 \mathbf{fft}_x(\psi^{\ell})) \right]_{N+1-k, N+1-n} + \left[ \mathbf{iff}_y (\Lambda^2 \mathbf{fft}_y(\psi^{\ell})) \right]_{N+1-k, N+1-n} \right) \\
&= \delta_x^2 \Re \left( \frac{1}{N+1} e^{-i\pi(N+1-k)} \sum_{p=0}^{N+1} \xi_p^2 \sum_{m=0}^{N+1} \psi_{m,N+1-n}^{\ell} e^{i\pi m} e^{-2\pi i \frac{mp}{N+1}} e^{2\pi i \frac{(N+1-k)p}{N+1}} \right) \\
&\quad + \delta_x^2 \Re \left( \frac{1}{N+1} e^{-i\pi(N+1-n)} \sum_{q=0}^{N+1} \lambda_q^2 \sum_{l=0}^{N+1} \psi_{N+1-k,l}^{\ell} e^{i\pi l} e^{-2\pi i \frac{lq}{N+1}} e^{2\pi i \frac{(N+1-n)q}{N+1}} \right) \\
&= \delta_x^2 \Re \left( \frac{1}{N+1} e^{-i\pi(N+1-k)} \sum_{p=0}^{N+1} \xi_p^2 \sum_{m=N+1}^{0} \psi_{g_2(n,m)}^{\ell} e^{i\pi(N+1-m)} e^{-2\pi i \frac{(N+1-m)p}{N+1}} e^{-2\pi i \frac{kp}{N+1}} e^{2\pi iq} \right) \\
&\quad + \delta_x^2 \Re \left( \frac{1}{N+1} e^{-i\pi(N+1-n)} \sum_{q=0}^{N+1} \lambda_q^2 \sum_{l=N+1}^{0} \psi_{g_2(l,k)}^{\ell} e^{i\pi(N+1-l)} e^{-2\pi i \frac{(N+1-l)q}{N+1}} e^{-2\pi i \frac{nq}{N+1}} e^{2\pi iq} \right) \\
&= \delta_x^2 \Re \left( \frac{1}{N+1} e^{-i\pi k} \sum_{p=0}^{N+1} \xi_p^2 \sum_{m=0}^{N+1} \psi_{n,m}^{\ell} e^{i\pi m} e^{-2\pi i \frac{mp}{N+1}} e^{2\pi i \frac{kp}{N+1}} \right) \\
&\quad + \delta_x^2 \Re \left( \frac{1}{N+1} e^{-i\pi n} \sum_{q=0}^{N+1} \lambda_q^2 \sum_{l=0}^{N+1} \psi_{l,k}^{\ell} e^{i\pi l} e^{-2\pi i \frac{lq}{N+1}} e^{2\pi i \frac{nq}{N+1}} \right) \\
&= \delta_x^2 \Re \left( \frac{1}{N+1} e^{-i\pi k} \sum_{p=0}^{N+1} \lambda_p^2 \sum_{m=0}^{N+1} \psi_{n,m}^{\ell} e^{i\pi m} e^{-2\pi i \frac{mp}{N+1}} e^{2\pi i \frac{kp}{N+1}} \right) \\
&\quad + \delta_x^2 \Re \left( \frac{1}{N+1} e^{-i\pi n} \sum_{q=0}^{N+1} \xi_q^2 \sum_{l=0}^{N+1} \psi_{l,k}^{\ell} e^{i\pi l} e^{-2\pi i \frac{lq}{N+1}} e^{2\pi i \frac{nq}{N+1}} \right) \\
&= \frac{\partial E_{kin}^{\Delta}(\psi^{\ell})}{\partial p_{n,k}^{\ell}}.
\end{align*}
$$

Similarly, we can prove for the partial derivative of the discrete kinetic energy with respect to $q_{n,k}$ that, for all $n, k \in \{0, \dots, N+1\}$, we have

$$
\frac{\partial E_{kin}^{\Delta}(\psi^{\ell})}{\partial q_{N+1-k,N+1-n}^{\ell}} = - \frac{\partial E_{kin}^{\Delta}(\psi^{\ell})}{\partial q_{n,k}^{\ell}}.
$$

This is due to the presence of the imaginary unit $i$ inside the real part of $\frac{\partial E_{kin}^{\Delta}(\psi^{\ell})}{\partial q_{n,k}}$ (see Proposition 11 of Chapter 2).
### Partial derivative of the discrete rotational energy:

Using the fact that $y_n = x_n$ and $x_n = -x_{N+1-n}$ for all $n \in \{0, \dots, N+1\}$, we have for the rotational partial derivative:

$$ 
\begin{aligned}
& \frac{\partial(E_r)_\varepsilon^\Delta(\psi^\ell)}{\partial p_{N+1-k,N+1-n}^\ell} \\
&= -\frac{2}{\Omega_\varepsilon} \delta_x^2 \Re \left( [X * \text{ifft}_y (\Lambda * \text{fft}_y(\psi^\ell)) - Y * \text{ifft}_x (\Xi * \text{fft}_x(\psi^\ell))]_{N+1-k,N+1-n} \right) \\
&= -\frac{2\Omega_\varepsilon \delta_x^2}{N+1} \Re \left( x_{N+1-k} e^{-i\pi(N+1-n)} \sum_{q=0}^{N+1} \lambda_q \sum_{l=0}^{N+1} \psi_{N+1-k,l}^\ell e^{i\pi l} e^{-2\pi i \frac{lq}{N+1}} e^{2\pi i \frac{(N+1-n)q}{N+1}} \right) \\
&\quad + \frac{2\Omega_\varepsilon \delta_x^2}{N+1} \Re \left( y_{N+1-n} e^{-i\pi(N+1-k)} \sum_{p=0}^{N+1} \xi_p \sum_{m=0}^{N+1} \psi_{m,N+1-n}^\ell e^{i\pi m} e^{-2\pi i \frac{mp}{N+1}} e^{2\pi i \frac{(N+1-k)p}{N+1}} \right) \\
&= -\frac{2\Omega_\varepsilon \delta_x^2}{N+1} \Re \left( -x_k e^{-i\pi(N+1-n)} \sum_{q=0}^{N+1} \sum_{l=N+1}^{0} \psi_{g_2(l,k)}^\ell e^{i\pi(N+1-l)} e^{-2\pi i \frac{(N+1-l)q}{N+1}} e^{2\pi i \frac{(N+1-n)q}{N+1}} \right) \\
&\quad + \frac{2\Omega_\varepsilon \delta_x^2}{N+1} \Re \left( -y_n e^{-i\pi(N+1-k)} \sum_{p=0}^{N+1} \sum_{m=N+1}^{0} \psi_{g_2(n,m)}^\ell e^{i\pi(N+1-m)} e^{-2\pi i \frac{(N+1-m)p}{N+1}} e^{2\pi i \frac{(N+1-k)p}{N+1}} \right) \\
&= \frac{2\Omega_\varepsilon \delta_x^2}{N+1} \Re \left( x_k e^{-i\pi n} \sum_{q=0}^{N+1} \sum_{l=0}^{N+1} \psi_{l,k}^\ell e^{i\pi l} e^{-2\pi i \frac{lq}{N+1}} e^{2\pi i \frac{nq}{N+1}} \right) \\
&\quad - \frac{2\Omega_\varepsilon \delta_x^2}{N+1} \Re \left( y_n e^{-i\pi k} \sum_{p=0}^{N+1} \sum_{m=0}^{N+1} \psi_{n,m}^\ell e^{i\pi m} e^{-2\pi i \frac{mp}{N+1}} e^{2\pi i \frac{kp}{N+1}} \right) \\
&= \frac{2\Omega_\varepsilon \delta_x^2}{N+1} \Re \left( y_k e^{-i\pi n} \sum_{q=0}^{N+1} \sum_{l=0}^{N+1} \psi_{l,k}^\ell e^{i\pi l} e^{-2\pi i \frac{lq}{N+1}} e^{2\pi i \frac{nq}{N+1}} \right) \\
&\quad - \frac{2\Omega_\varepsilon \delta_x^2}{N+1} \Re \left( x_n e^{-i\pi k} \sum_{p=0}^{N+1} \sum_{m=0}^{N+1} \psi_{n,m}^\ell e^{i\pi m} e^{-2\pi i \frac{mp}{N+1}} e^{2\pi i \frac{kp}{N+1}} \right) \\
&= \frac{\partial(E_r)_\varepsilon^\Delta(\psi^\ell)}{\partial p_{n,k}^\ell}
\end{aligned}
$$

Similarly, we can prove for the partial derivative of the discrete rotational energy with respect to $q_{n,k}$ that, for all $n, k \in \{0, \dots, N+1\}$, we have

$$ \frac{\partial E_r^\Delta(\psi^\ell)}{\partial q_{N+1-k,N+1-n}^\ell} = - \frac{\partial E_r^\Delta(\psi^\ell)}{\partial q_{n,k}^\ell}. $$

Again, this is due to the presence of the imaginary unit $i$ inside the real part of $\frac{\partial E_r^\Delta(\psi^\ell)}{\partial q_{n,k}}$ (see Proposition 12 of Chapter 2).
**Partial derivative of the discrete confinement energy:**

Last, using Proposition 38 of this chapter, we can prove that

$$ \frac{\partial(E_W)^{\Delta}_{\varepsilon,\delta}(\psi^1, \psi^2)}{\partial p_{N+1-k,N+1-n}^{\ell}} = \frac{\partial(E_W)^{\Delta}_{\varepsilon,\delta}(\psi^1, \psi^2)}{\partial p_{n,k}^{\ell}}, \\ \frac{\partial(E_W)^{\Delta}_{\varepsilon,\delta}(\psi^1, \psi^2)}{\partial q_{N+1-k,N+1-n}^{\ell}} = -\frac{\partial(E_W)^{\Delta}_{\varepsilon,\delta}(\psi^1, \psi^2)}{\partial q_{n,k}^{\ell}}. $$

Putting together all the previous results, we obtain $\psi_{n,k}^{\ell} = \overline{\psi_{N+1-k,N+1-n}^{\ell}}$ for all $(n,k) \in \{1, \dots, N\}^2$ and $\ell = 1, 2$. This completes the proof of Theorem 41. □

**Example 42** In this example, we display the results of the gradient descent method (see Section 4 of Chapter 2) with one component $\psi$ using 19 points in the discretization after 1 iteration. We consider a real initial datum depending only on the distance from the center hence satisfying the condition of Theorem 41. As we can see in the table of Figure 4.2, the symmetry property of the real and imaginary parts of $\phi$ is preserved with respect to the second diagonal. The other three possible symmetries are not preserved by the numerical discretization.

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

<table><caption>Figure 4.2 – Example of symmetry conservation using the fast Fourier transformation approach.</caption><thead><tr><th></th><th>j = 1</th><th>j = 2</th><th>j = 3</th><th>j = 4</th></tr></thead><tbody><tr><th>max<sub>n,k</sub> |R(φ<sub>n,k</sub> - φ&#773;<sub>g<sub>j</sub>(n,k)</sub>)|</th><td>3 × 10<sup>-5</sup></td><td>4.7 × 10<sup>-16</sup></td><td>1.7 × 10<sup>-5</sup></td><td>1.7 × 10<sup>-5</sup></td></tr><tr><th>max<sub>n,k</sub> |Im(φ<sub>n,k</sub> - φ&#773;<sub>g<sub>j</sub>(n,k)</sub>)|</th><td>10<sup>-5</sup></td><td>5.3 × 10<sup>-20</sup></td><td>10<sup>-5</sup></td><td>10<sup>-5</sup></td></tr><tr><th>max<sub>n,k</sub> ||φ<sub>n,k</sub>||<sup>2</sup> - ||φ<sub>g<sub>j</sub>(n,k)</sub>||<sup>2</sup></th><td>5.6 × 10<sup>-5</sup></td><td>3 × 10<sup>-16</sup></td><td>3.5 × 10<sup>-5</sup></td><td>3.5 × 10<sup>-5</sup></td></tr></tbody></table>

**Remark 43** Note that visually, the first diagonal and the second diagonal are inverted. In fact, at the bottom left of each figure we have the couple $(x_0, y_0) = (-7, -7)$ and at the top right we have the couple $(x_{N+1}, y_{N+1}) = (7, 7)$. Since we define the first diagonal as the definition of a matrix diagonal (see Definition 35), visually it is associated with the diagonal connecting the couples $(x_0, y_0)$ and $(x_{N+1}, y_{N+1})$.
**Theorem 44** Let $v \in \mathbb{C}^{N^2}$ and let $\text{fft}_x, \text{fft}_y$ denote the Fourier transform defined in Section 2.1 of Chapter 2, then we have

$$
\begin{align*}
\delta_x^2 \sum_{p=0}^{N} \sum_{k=0}^{N} |(\text{fft}_x(v))_{p,k}|^2 &= (N+1)\delta_x^2 \sum_{n=0}^{N} \sum_{k=0}^{N} |v_{n,k}|^2, \\
\delta_x^2 \sum_{q=0}^{N} \sum_{n=0}^{N} |(\text{fft}_y(v))_{n,q}|^2 &= (N+1)\delta_x^2 \sum_{n=0}^{N} \sum_{k=0}^{N} |v_{n,q}|^2, \\
\delta_x^2 \sum_{n=0}^{N} \sum_{k=0}^{N} |(\text{ifft}_x(v))_{n,k}|^2 &= \frac{\delta_x^2}{N+1} \sum_{n=0}^{N} \sum_{k=0}^{N} |v_{n,k}|^2, \\
\delta_x^2 \sum_{k=0}^{N} \sum_{n=0}^{N} |(\text{ifft}_y(v))_{n,k}|^2 &= \frac{\delta_x^2}{N+1} \sum_{n=0}^{N} \sum_{k=0}^{N} |v_{n,k}|^2.
\end{align*}
$$

*Proof:* We have, with $\delta$ the Kronecker symbol,

$$
\begin{align}
\delta_x^2 \sum_{p=0}^{N} |(\text{fft}_x(v))_{p,k}|^2 &= \delta_x^2 \sum_{p=0}^{N} \left( \sum_{m=0}^{N} v_{m,k} e^{i\pi m} e^{-\frac{2i\pi}{N+1}mp} \sum_{l=0}^{N} \overline{v_{l,k}} e^{-i\pi l} e^{\frac{2i\pi}{N+1}lp} \right) \nonumber \\
&= \delta_x^2 \sum_{p=0}^{N} \sum_{m=0}^{N} \sum_{l=0}^{N} v_{m,k} \overline{v_{l,k}} e^{i\pi(m-l)} e^{-\frac{2i\pi}{N+1}(m-l)p} \nonumber \\
&= \delta_x^2 \sum_{m=0}^{N} \sum_{l=0}^{N} v_{m,k} \overline{v_{l,k}} e^{i\pi(m-l)} \sum_{p=0}^{N} e^{-\frac{2i\pi}{N+1}(m-l)p} \tag{4.5} \\
&= (N+1)\delta_x^2 \sum_{m=0}^{N} \sum_{l=0}^{N} v_{m,k} \overline{v_{l,k}} e^{i\pi(m-l)} \delta_{m,l} \nonumber \\
&= (N+1)\delta_x^2 \sum_{m=0}^{N} v_{m,k} \overline{v_{m,k}}. \nonumber
\end{align}
$$

Finally, by summing (4.5) over $k \in \{0, \dots, N\}$, we obtain

$$
\delta_x^2 \sum_{p=0}^{N} \sum_{k=0}^{N} |(\text{fft}_x(v))_{p,k}|^2 = (N+1)\delta_x^2 \sum_{m=0}^{N} \sum_{k=0}^{N} |v_{m,k}|^2.
$$

Similarly, we can prove the rest of the equalities of Theorem 44. $\square$
Conclusion

I have considered in this thesis a dimensionless Gross-Pitaevskii energy as a model for rotating one component and two components Bose-Einstein condensates in two dimensions in a strong confinement regime.

First, I have introduced a new discretization of this energy using both finite
difference and FFT approaches which allowed using predefined sparse matrices to
speed up the computation.

Second, in contrast to the literature, I proposed a minimization method for this discrete energy using an explicit $L^2$ gradient method with projection (EPG). This method allowed for the derivation of a stopping criterion.

Third, I introduced two post processing algorithms for the numerical minimizers. One is aimed for the single vortices and the other for vortex sheets. Both allow to detect these structures and compute their indices.

Fourth, I have ran the EPG method and algorithms for different physical regimes from one component with high rotation to two components in coexistence and segregation regimes. This allowed to validate recent theoretical results as well as to support conjectures (as for example the existence of vortex sheets in the segregation regime) and cover different physical cases.

Fifth, I have compared the efficiency of EPG to the GPELab method [11] in both cases of one component and two components condensates. On all the numerical tests presented with a moderate number of points, EPG appears to be roughly speaking 3 times faster than GPELab.

Last, I have adapted some theorems that are known in the continuous setting to a discrete setting, like the existence of a global minimizer of the energy (2.1). Then I have studied a symmetry property appearing in some of the simulations of Chapter 3 using both finite difference and FFT approaches.
Bibliography

[1] Ground states of two-component bose-einstein condensates with an internal atomic josephson junction. *East Asian Journal on Applied Mathematics*, 1(1):49-81, 2018.

[2] A. Aftalion and Q. Du. Vortices in a rotating Bose-Einstein condensate: Critical angular velocities and energy diagrams in the Thomas-Fermi regime. *Phys. Rev. A*, 64:063603, Nov 2001.

[3] A. Aftalion, B. Noris, and C. Sourdis. Thomas-Fermi approximation for coexisting two component Bose-Einstein condensates and nonexistence of vortices for small rotation. *Com. Math.Physics*, page a paraitre, 2014.

[4] A. Aftalion and T. Riviere. Vortex energy and vortex bending for a rotating Bose-Einstein condensate. *Physical Review A*, 64:043611, 2001. 7 pages, 2 figures. submitted.

[5] A. Aftalion and J. Royo-Letelier. A minimal interface problem arising from a two component Bose-Einstein condensate via $\Gamma$-convergence. *Calc. Var. Partial Differential Equations*, 52(1-2):165-197, 2015.

[6] A. Aftalion and E. Sandier. Vortex patterns and sheets in segregated two component Bose-Einstein condensates. *Calc. Var. Partial Differential Equations*, 59(1):Paper No. 19, 38, 2020.

[7] B. P. Anderson, P. C. Haljan, C. E. Wieman, and E. A. Cornell. Vortex precession in Bose-Einstein condensates: Observations with filled and empty cores. *Phys. Rev. Lett.*, 85:2857-2860, Oct 2000.

[8] X. Antoine and R. Duboscq. GPELab toolbox. http://gpelab.math.cnrs.fr/.

[9] X. Antoine and R. Duboscq. GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations i: computation of stationary solutions. *Computer Physics Communications*, 185:2969-2991, 11 2014.

[10] X. Antoine and R. Duboscq. Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates. *J. Comput. Phys.*, 258:509-523, 2014.

[11] X. Antoine and R. Duboscq. GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations ii: Dynamics and stochastic simulations. *Computer Physics Communications*, 193:95-117, 2015.
[12] W. Bao. Ground states and dynamics of multicomponent bose-einstein condensates. *Multiscale Modeling & Simulation*, 2(2):210-236, 2004.

[13] W. Bao and Y. Cai. Mathematical theory and numerical methods for Bose-Einstein condensation. *Kinet. Relat. Models*, 6(1):1-135, 2013.

[14] W. Bao, Y. Cai, and H. Wang. Efficient numerical methods for computing ground states and dynamics of dipolar bose-einstein condensates. *Journal of Computational Physics*, 229(20):7874-7892, 2010.

[15] W. Bao, I-L. Chern, and F. Y. Lim. Efficient and spectrally accurate numerical methods for computing ground and first excited states in Bose-Einstein condensates. *Journal of Computational Physics*, 219(2):836-854, 2006.

[16] W. Bao and Q. Du. Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow. *SIAM J. Sci. Comput.*, 25(5):1674-1697, 2004.

[17] W. Bao, D. Jaksch, and P. A. Markowich. Numerical solution of the Gross-Pitaevskii equation for Bose-Einstein condensation. *Journal of Computational Physics*, 187(1):318-342, may 2003.

[18] W. Bao, S. Jiang, Q. Tang, and Y. Zhang. Computing the ground state and dynamics of the nonlinear schrödinger equation with nonlocal interactions via the nonuniform fft. *Journal of Computational Physics*, 296:72-89, 2015.

[19] W. Bao and W. Tang. Ground-state solution of Bose-Einstein condensate by directly minimizing the energy functional. *Journal of Computational Physics*, 187(1):230-254, may 2003.

[20] W. Bao and W. Tang. Ground-state solution of bose-einstein condensate by directly minimizing the energy functional. *Journal of Computational Physics*, 187(1):230-254, 2003.

[21] W. Bao and H. Wang. An efficient and spectrally accurate numerical method for computing dynamics of rotating Bose-Einstein condensates. *Journal of Computational Physics*, 217(2):612-626, 2006.

[22] S. N. Bose. Plancks gesetz und lichtquantenhypothese. *Zeitschrift für Physik*, 26:178-181, 1924.

[23] T. Byrnes. Quantum computation using two component Bose-Einstein condensates. *International Journal of Nuclear and Quantum Engineering*, 6(3):363-368, 2012.

[24] T. Byrnes, S. Koyama, K. Yan, and Y. Yamamoto. Neural networks using two-component Bose-Einstein condensates. 2012.

[25] M. M. Cerimele, M. L. Chiofalo, F. Pistella, S. Succi, and M. P. Tosi. Numerical solution of the Gross-Pitaevskii equation using an explicit finite-difference scheme: An application to trapped Bose-Einstein condensates. *Phys. Rev. E*, 62:1382-1389, Jul 2000.
[26] E. A. Cornell and C. E. Wieman. Nobel lecture: Bose-Einstein condensation in a dilute gas, the first 70 years and some recent experiments. *Rev. Mod. Phys.*, 74:875-893, Aug 2002.

[27] M. Correggi, F. Pinsker, N. Rougerie, and J. Yngvason. Critical rotational speeds in the Gross-Pitaevskii theory on a disc with Dirichlet boundary conditions. *J. Stat. Phys.*, 143(2):261-305, 2011.

[28] I. Danaila and P. Kazemi. A new Sobolev gradient method for direct minimization of the Gross-Pitaevskii energy with rotation. *SIAM Journal on Scientific Computing*, 32(5):2447-2467, 2010.

[29] A. Einstein. *Quantentheorie des einatomigen idealen Gases*, pages 237-244. 2005.

[30] A. Fetter. Vortices and dynamics in trapped Bose-Einstein condensates. *Journal of Low Temperature Physics*, 161:445-459, 12 2010.

[31] M. Goldman and B. Merlet. Phase segregation for binary mixtures of Bose-Einstein condensates. *SIAM J. Math. Anal.*, 49(3):1947-1981, 2017.

[32] M. Goldman and J. Royo-Letelier. Sharp interface limit for two components Bose-Einstein condensates. *ESAIM Control Optim. Calc. Var.*, 21(3):603-624, 2015.

[33] E. P. Gross. Hydrodynamics of a superfluid condensate. *Journal of Mathematical Physics*, 4(2):195-207, 1963.

[34] D. S. Hall, M. R. Matthews, J. R. Ensher, C. E. Wieman, and E. A. Cornell. Dynamics of component separation in a binary mixture of Bose-Einstein condensates. *Phys. Rev. Lett.*, 81:1539-1542, Aug 1998.

[35] A. Hansen, L. Leslie, M. Bhattacharya, and N. Bigelow. Transfer and storage of optical information in a spinor Bose-Einstein condensate. 10 2010.

[36] T. L. Ho. Bose-Einstein condensates with large number of vortices. *Physical review letters*, 87:060403, 09 2001.

[37] R. Ignat and V. Millot. The critical velocity for vortex existence in a two-dimensional rotating Bose-Einstein condensate. *Journal of Functional Analysis*, 233(1):260-306, 2006.

[38] K. Kasamatsu and M. Tsubota. Vortex sheet in rotating two-component Bose-Einstein condensates. *Phys. Rev. A*, 79:023606, Feb 2009.

[39] K. Kasamatsu, M. Tsubota, and M. Ueda. Vortex phase diagram in rotating two-component Bose-Einstein condensates. *Phys. Rev. Lett.*, 91:150406, Oct 2003.

[40] K. Kasamatsu, M. Tsubota, and M. Ueda. Vortices in multicomponent Bose-Einstein condensates. *International Journal of Modern Physics B*, 19(11):1835-1904, apr 2005.
[41] J. R. Li, J. Lee, W. Huang, S. Burchesky, B. Shteynas, F. Ç. Top, A.O. Jamison, and W. Ketterle. A stripe phase with supersolid properties in spin-orbit-coupled Bose-Einstein condensates. *Nature*, 543(7643):91-94, mar 2017.

[42] O. V. Lounasmaa and E. Thuneberg. Vortices in rotating superfluid $^3$He. *Proceedings of the National Academy of Sciences*, 96(14):7760-7767, 1999.

[43] K. W. Madison, F. Chevy, W. Wohlleben, and J. Dalibard. Vortex formation in a stirred Bose-Einstein condensate. *Phys. Rev. Lett.*, 84:806-809, Jan 2000.

[44] B. Malomed. *Multi-Component Bose-Einstein Condensates: Theory*, pages 287-305. Springer Berlin Heidelberg, Berlin, Heidelberg, 2008.

[45] P. Mason and A. Aftalion. Classification of the ground states and topological defects in a rotating two-component Bose-Einstein condensate. *Phys. Rev. A*, 84:033611, Sep 2011.

[46] M. R. Matthews, B. P. Anderson, P. C. Haljan, D. S. Hall, C. E. Wieman, and E. A. Cornell. Vortices in a Bose-Einstein condensate. *Phys. Rev. Lett.*, 83:2498-2501, Sep 1999.

[47] T. Mizushima, N. Kobayashi, and K. Machida. Coreless and singular vortex lattices in rotating spinor Bose-Einstein condensates. *Physical Review A*, 70(4), oct 2004.

[48] L. P. Pitaevskii. Vortex linex in an imperfect Bose gas. *Soviet Physics-JETP*, 13(2):451-454, 1961.

[49] Y. Shin, M. Saba, M. Vengalattore, T. A. Pasquini, C. Sanner, A. E. Leanhardt, M. Prentiss, D. E. Pritchard, and W. Ketterle. Dynamical instability of a doubly quantized vortex in a Bose-Einstein condensate. *Physical Review Letters*, 93(16), oct 2004.

[50] W. F. Vinen. The physics of superfluid helium. 2004.

[51] Q. Z. Zhu and B. Wu. Superfluidity of Bose-Einstein condensates in ultracold atomic gases. *Chinese Physics B*, 24(5):050507, apr 2015.