# Numerical simulation of Bose-Einstein condensates

Anthony Nahas

► To cite this version:

Anthony Nahas. Numerical simulation of Bose-Einstein condensates. Mathematics [math]. Université de Lille, 2022. English. NNT: 2022ULILB026. tel-03992001

HAL Id: tel-03992001

https://hal.science/tel-03992001

Submitted on 16 Feb 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
École Doctorale MADIS-631
Université de Lille – Laboratoire Paul Painlevé

THÈSE DE DOCTORAT

présentée pour l'obtention du grade de DOCTEUR
DE L'UNIVERSITÉ DE LILLE

par

Anthony NAHAS

Simulation numérique de
condensats de Bose-Einstein

Numerical simulation of
Bose-Einstein condensates

Discipline: Mathématiques Appliquées

Soutenue le 11 Octobre 2022 devant le jury composé par:

Laurent DI MENZA

Professeur

Université de Reims

Rapporteur

Claire SCHEID

MCF HDR

Université de Nice

Rapportrice

Jean-Claude GARREAU

Directeur de recherche

Université de Lille

Examinateur

Christophe GEUZAINE

Professeur

Université de Liège

Examinateur

Guillaume DUJARDIN

CR HDR

Université de Lille

Co-directeur

Ingrid LACROIX-VIOLET

Professeure

Université de Lorraine

Co-directrice
BEC Art

Check out the art presented below that has been randomly created through-
out the work of my thesis. These are beautiful numerical minimizers obtained by
simulating a two components Bose-Einstein condensate in a segregation regime.

<figure><img src="image_2.png" /><figcaption>Segregated two components BEC fish</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Two components BEC version of Thalia & Melpomene masks</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Butterfly breaking free from a two components BEC</figcaption></figure>

Enjoy.
Acknowledgment

First of all, I would like to acknowledge both my supervisors Guil-
laume DUJARDIN and Ingrid LACROIX-VIOLET for their guidance,
enthusiasm and availability over the last few years, making my PhD
interesting and enjoyable. I also want to thank them for their patience
and for treating me as a part of their family.

I warmly thank Laurent DI MENZA and Claire SCHEID for ac-
cepting the role of reviewers for this thesis. I would also like to thank
Jean-Claude GARREAU and Christophe GEUZAINE to have taken
part in the jury.

My gratitude is extended to Laboratoire Paul Painlevé and to Inria
for their hospitality and kindness. This could not be complete without
mentioning Labex CEMPI for funding my 2nd year of masters as well
as co-funding my thesis with Région Hauts de France.

I want to say a huge thank you to my parents Elias and Nawal for
encouraging me to always pursue my dreams and for being next to
me in the hardest of times and supporting me both financially and
morally. To my brother and sister, Elie and Maria, for always being
by my side.

And finally, to my one and only Rita MSAWBAH, for her endless
patience, understanding, support and love.
Abstract

In this thesis, we consider a Gross-Pitaevskii (GP) energy functional as a model for rotating one component and two components Bose-Einstein condensates (BEC) in two dimensions. This model can be non-dimensionalized to highlight a strong confinement regime with strong interaction between the two components. We introduce a new discretization of this energy, featuring both finite difference and fast Fourier transform approaches, in a bounded domain of $\mathbb{R}^2$ using Dirichlet boundary conditions. We develop an explicit gradient method algorithm with adaptive step and projection (EPG) over the constraints manifold for the minimization of the discrete energy. This method allows for the derivation of a stopping criterion. Moreover, we propose two post processing algorithms for the numerical minimizers. One is aimed for single vortices while the other is aimed for vortex sheets. Both algorithms detect these structures and compute their indices.

In a recent article [6], the authors study the behaviour of a segregated two components BEC set into rotation. They are able to prove that for large rotation, the interface between the components gets long, conjecturing the possibility towards vortex sheets. They also study the vortex structures of BEC in a segregated regime. In this thesis, we produce numerical simulations using EPG, validating these recent theoretical results, supporting conjectures, and covering different physical cases (the cases of one component and two components in coexistence regime [45, 39]). We also illustrate the efficiency of EPG compared to the well-known GPELab method [11] which requires to solve a linearly implicit system at each step.

Finally, we adapt a few theorems found in the literature to our discrete setting. We prove the existence of a global minimizer of the GP energy functional for the finite difference scheme and study some of its properties. We also work on symmetrical problems we encountered in some of the numerical simulations.
Résumé

Dans cette thèse, nous considérons une fonctionnelle d'énergie Gross-Pitaevskii (GP) comme modèle d'une espèce et de deux espèces de condensats de Bose-Einstein (BEC) en deux dimensions mise en rotation. Ce modèle peut être non dimensionné pour mettre en évidence un régime de fort confinement avec une forte interaction entre les deux espèces. Nous introduisons une nouvelle discrétisation de cette énergie, comportant à la fois des approches par différences finies et par transformée de Fourier, dans un domaine borné de $\mathbb{R}^2$ en utilisant des conditions de Dirichlet nulles aux bords. Nous développons ainsi un algorithme de méthode de gradient explicite avec pas adaptatif et projection (EPG) sur la variété de contraintes pour la minimisation de l'énergie discrète. Cette méthode permet de dériver un critère d'arrêt. De plus, nous proposons deux algorithmes de post-traitement pour les minimiseurs numériques. L'un est destiné aux vortex simples tandis que l'autre est destiné aux nappes de vortex. Les deux algorithmes détectent ces structures et calculent leurs indices.

Dans un article récent [6], les auteurs étudient le comportement d'un BEC à deux espèces en ségrégation et mis en rotation. Ils sont capables de prouver que pour une forte rotation, l'interface entre les espèces devient longue, conjecturant la possibilité vers des nappes de vortex. Ils étudient également les structures des vortex du BEC en régime de ségrégation. Dans cette thèse, nous produisons des simulations numériques à l'aide d'EPG, validant ces résultats théoriques récents, supportant des conjectures, et couvrant différents cas physiques (les cas d'une espèce et de deux espèces en régime de co-existence [45, 39]). Nous illustrons également l'efficacité d'EPG par rapport à la méthode bien connue GPELab [11] qui nécessite de résoudre un système linéaire implicite à chaque étape.

Enfin, nous adaptons quelques théorèmes trouvés dans la littérature à notre cadre discret. Nous prouvons l'existence d'un minimiseur global de la fonctionnelle d'énergie GP pour le schéma aux différences finies et étudions certaines de ses propriétés. Nous travaillons également sur des problèmes symétriques rencontrés dans certaines simulations numériques.
<table>
  <tbody>
    <tr>
      <td>3</td>
      <td>Numerical results for a discretization of the energy using Fast Fourier Transform</td>
      <td>49</td>
    </tr>
    <tr>
      <td>3.1</td>
      <td>One component condensate without rotation</td>
      <td>49</td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>One component condensate with rotation</td>
      <td>50</td>
    </tr>
    <tr>
      <td>3.3</td>
      <td>Two components condensate without rotation (Ω = 0)</td>
      <td>56</td>
    </tr>
    <tr>
      <td>3.4</td>
      <td>Two components condensate with rotation (Ω ≠ 0) in a seg-regation regime (δ<sub>ε</sub> > 1)</td>
      <td>58</td>
    </tr>
    <tr>
      <td>3.5</td>
      <td>Two components condensate in the coexistence regime δ<sub>ε</sub> < 1</td>
      <td>62</td>
    </tr>
    <tr>
      <td>4</td>
      <td>Comparison between finite difference and Fourier transform approaches</td>
      <td>68</td>
    </tr>
    <tr>
      <td>5</td>
      <td>GPELab comparison</td>
      <td>69</td>
    </tr>
    <tr>
      <td>5.1</td>
      <td>One component condensate comparison</td>
      <td>69</td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>Two components condensate comparison</td>
      <td>76</td>
    </tr>
    <tr>
      <td>5.3</td>
      <td>Conclusion</td>
      <td>78</td>
    </tr>
    <tr>
      <td>4</td>
      <td>Theoretical part</td>
      <td>81</td>
    </tr>
    <tr>
      <td>1</td>
      <td>Theory for the finite difference discretization</td>
      <td>81</td>
    </tr>
    <tr>
      <td>1.1</td>
      <td>Theoretical studies of the energy E<sub>ε</sub><sup>Δ</sup></td>
      <td>81</td>
    </tr>
    <tr>
      <td>1.2</td>
      <td>Symmetric theory</td>
      <td>85</td>
    </tr>
    <tr>
      <td>2</td>
      <td>Theory for the Fast Fourier transform discretization</td>
      <td>92</td>
    </tr>
    <tr>
      <td></td>
      <td>Conclusion</td>
      <td>97</td>
    </tr>
  </tbody>
</table>
# Chapter 1

## Introduction

In condensed matter physics, a Bose-Einstein condensate (BEC) is a state of matter that is typically formed when a gas of bosons at very low densities is cooled to temperatures very close to absolute zero. Under these conditions, a large fraction of bosons occupy the lowest quantum state, at which point microscopic quantum mechanical phenomena, particularly wave function interference, become apparent macroscopically. Thus, a single wave function is sufficient to describe each component of the BEC. When the so called Bose-Einstein condensates are set to rotation, under a strong confinement regime, topological defects often manifest themselves as vortices that correspond to a zero of the wave function with phase circulation. This phenomenon was first observed in two components BEC [46].

When a BEC is set to a high rotation in a strong confinement regime, the vortices align and form unique structures. In the case of a single component condensate, we observe singly quantized vortices, forming triangular lattices once they are numerous [36, 45]. In the case of a two components condensate, depending on the interaction between the two components, many structures may appear. For example, we can observe coreless vortices, which refer to having singly quantized vortices in one component while having a corresponding peak in the second component [47], or vortex sheets [45, 40].

Different models of BEC with one or several components have already been studied in the mathematical literature. For example, the minimization of the Gross-Pitaevskii functional in $R^2$ is studied theoretically in [37] (see also references therein). In [6, 3], the authors study different structures of a BEC in a strong confinement and coupling regime in a bounded domain of $R^2$. Several methods have been developed for the numerical computation of approximations of minimizers of Gross-Pitaevskii energies. For example, in [20, 12, 16, 21, 15, 14, 1, 13, 10, 18], the authors develop numerical methods which require solving a linear system at each time-step. Another option is to use Sobolev gradients, as opposed to $L^2$ gradients, as developed in [28].

The outline of this thesis is as follows. The first chapter is devoted to defining the problematic of this thesis. In Chapter 1 Section 1, we introduce the physical problem behind BEC as well as its uses in the technological intervention while recalling the main advances to it. We then introduce in Section 2 the different continuous models for describing two components rotating BEC that has been studied in the latest mathematical literature. In Section 4.1, we introduce the continuous model used in the rest of this work for rotating Bose-Einstein condensates with two components.
Finally in Section 4, we recall the different regimes for single and two components condensates.

In order to numerically study these regimes, we discretize in Chapter 2 Sections 1 and 2 (using finite difference and FFT approaches respectively) the continuous Gross-Pitaevskii energy. For the finite difference scheme, we re-write the energy in Section 1.3 using predefined sparse matrices so that we optimize the computation time of the numerical simulations. For the FFT approach, our discretization uses the predefined FFT algorithm of Python for a faster computation of the gradient. We introduce in Section 4 a gradient method for the minimization of the energy with an adaptive step and a projection step to take the constraints into account. One of the interests of this method is that it allows for the derivation of a stopping criterion which we developed in Section 3. We also develop in Chapter 2 Section 4 two post processing algorithms for the computation of indices. The first deals with vortices and the second deals with vortex sheets.

In Chapter 3, we present numerical results for the regimes described in Sections 4.2, 4.3 and 4.4 of Chapter 1. In Section 2 we present the numerical results of the finite difference scheme developed in Chapter 1 Section 1. We can see in some of the simulations (especially in the one component case with a huge rotational speed) a grid orientation effect. However, in Section 3, we present the numerical results of the Fast Fourier Transformation scheme developed in Chapter 1 Section 2. Looking at the numerical simulations, we can clearly see that there is no longer a grid orientation effect. In particular, we validate numerically, using both approaches, recent theoretical results and we support some conjectures as for example the existence of vortex sheets in a segregation regime. The last Section 5 is devoted to the comparison of the efficiency of our method to that of GPELab [9, 11].

Finally, in Chapter 4, we adapt few theorems we already have for the continuous problem into our discrete analogue for the finite difference scheme in Section 1. We also prove a symmetrical problem we encountered in some of the simulations of Chapter 3 Section 2. In Section 2, we prove another symmetrical problem encountered using the FFT scheme.

## 1 The physical problem

Bose-Einstein condensate (BEC), a state of matter in which separate atoms or subatomic particles, cooled to near absolute zero (0 Kelvin), coalesce into a single quantum mechanical entity—that is, one that can be described by a single wave function—on a near-macroscopic scale. This form of matter was predicted in 1925 by Albert Einstein [29] on the basis of the quantum formulations of the Indian physicist Satyendra Nath Bose [22].

Although it had been predicted for decades, the first atomic BEC was realized only in 1995, when Eric Cornell and Carl Wieman of JILA, a research institution jointly operated by the National Institute of Standards and Technology (NIST) and the University of Colorado at Boulder, cooled a gas of rubidium atoms to $1.7 \times 10^{-7}$ K above absolute zero (Figure 1.1). Along with Wolfgang Ketterle of the Massachusetts Institute of Technology (MIT), who created a BEC with sodium atoms, these researchers received the 2001 Nobel Prize for Physics. Research on BECs has expanded the understanding of quantum physics and has led to the discovery of
new physical effects such as turning a metal into an insulator (an effect known as
Anderson localization) and solids that flow through themselves [41].

<figure><img src="image_5.png" /><figcaption>Figure 1.1 – Precessing two components vortex (a) direct images, at 50 ms intervals (b) smoothed images of the first Bose Einstein condensate experiment ([30]).</figcaption></figure>

BEC theory traces back to 1924, when Bose considered how groups of photons behave. Photons belong to one of the two great classes of elementary or submicroscopic particles defined by whether their quantum spin is a non-negative integer (0, 1, 2, ... ) or an odd half integer (1/2, 3/2, ... ). The former type, called bosons, includes photons, whose spin is 1. The latter type, called fermions, includes electrons, whose spin is 1/2.

As Bose noted, the two classes behave differently. According to the Pauli exclusion principle, fermions tend to avoid each other, for which reason each electron in a group occupies a separate quantum state (indicated by different quantum numbers, such as the electron’s energy). In contrast, an unlimited number of bosons can have the same energy state and share a single quantum state.

Einstein soon extended Bose's work to show that at extremely low temperatures "bosonic atoms" with even spins would coalesce into a shared quantum state at the lowest available energy. The requisite methods to produce temperatures low enough to test Einstein's prediction did not become attainable, however, until the 1990s. One of the breakthroughs depended on the novel technique of laser cooling and trapping, in which the radiation pressure of a laser beam cools and localizes atoms by slowing them down. The second breakthrough depended on improvements in magnetic confinement in order to hold the atoms in place without a material container. Using these techniques, Cornell and Wieman [26] succeeded in merging about 2000 individual atoms into a "superatom," a condensate large enough to observe with a microscope, that displayed distinct quantum properties (Figure 1.2).

Superfluidity, a remarkable macroscopic quantum phenomenon, was first discovered in the study of liquid helium 4 in 1938. It characterizes the property of a fluid with zero viscosity which therefore flows without any loss of kinetic energy. When stirred, a superfluid forms vortices that continue to rotate indefinitely, meaning it cannot fade away or disappear, it is only allowed to move out of the superfluid or annihilate with a vortex of opposite direction [42, 49]. Although it is found theoretically that superfluidity is a general phenomenon for interacting boson systems, liquid helium had been the only bosonic superfluid available in experiments until the first experimental realization of BEC of dilute alkali atomic gases. This addition to the family of superfluids is highly nontrivial as BECs offer various aspects
<figure><img src="image_6.png" /><figcaption>Figure 1.2 – Large triangular vortex array in a rotating BEC ([30]).</figcaption></figure>

of advantages over liquid helium that can greatly enrich our understanding of superfluidity [51, 50]. Almost all the parameters of a BEC can be controlled easily in experiments: its kinetic energy, density, and the interaction between atoms can all be tuned easily by engineering the atom-laser interaction, magnetic or optical traps, and the Feshbach resonance.

One of the most intriguing property of BECs is that they can slow down light. In 1998 Lene Hau of Harvard University and her colleagues slowed light traveling through a BEC from its speed in vacuum of $3 \times 10^8$ meters per second to a mere 17 meters per second, or about 60 kilometers per hour. Since then, Hau and others have completely halted and stored a light pulse within a BEC, later releasing the light unchanged or sending it to a second BEC. These manipulations hold promise for new types of light-based telecommunications, optical storage of data [35], and quantum computation using qubits made of two components Bose-Einstein condensates [23], though the low-temperature requirements of BECs offer practical difficulties. Also it can be used in the field of computer vision and deep learning since the properties of a BEC allow the creation of a neural network (see [24]).

Among various fascinating features demonstrated by BECs, much interest has been drawn to experiments performed in two components condensates composed of $^{87}$Rb atoms in two different spin states. In addition to this most familiar type of the two components condensate, more sophisticated experimental techniques open a way to experiments with hetero nuclear mixtures, such as $^{41}$K-$^{85}$Rb, $^{39}$K-$^{85}$Rb, and $^{85}$Rb-$^{87}$Rb. Actual experimental results reported thus far for two components BEC mixtures are limited. They include the demonstration of segregation between two species in immiscible BECs [34], and various vortex configurations in $^{87}$Rb, with one spin-state species filling an empty core of a vortex created in the other species [7] (see [44] for more details).

## 2 Gross-Pitaevskii equations

In this section, we discuss the different physical and mathematical models describing a two components BEC under rotation and confinement.
At a very low temperature $\mathcal{T}$, below a certain critical temperature $\mathcal{T}_c$ given by

$$ \mathcal{T}_c \approx 3.3125 \frac{\hbar^2 n_\beta^2}{m k_B}, $$

where $m$ represents the mass per boson, $n$ the particle density, $\hbar$ the reduced Planck constant and $k_B$ the Boltzmann constant, the dilute multi-component Bose-Einstein condensates can be represented by their macroscopic wave functions (each component is represented by a single wave function) which satisfies a nonlinear Schrödinger equation, called the Gross-Pitaevskii equation (known as GPE), based on the work of Eugene P. Gross [33] and Lev P. Pitaevskii [48]. For $\ell \in \{1, 2\}$, the two dimensional two components GPE is given by

$$ i\hbar\partial_t u_\ell(\mathbf{x}, t) = \left[ -\frac{\hbar^2}{2m_\ell} \nabla^2 - i\hbar\Omega\mathbf{x}^\perp\nabla + V_\ell(r) + U_\ell|u_\ell|^2 + U_{12}|u_{3-\ell}|^2 \right] u_\ell(\mathbf{x}, t), \quad (1.1) $$

under the mass $^1$ constraints

$$ \int_D |u_\ell|^2 = N_\ell \quad \text{where} \quad N_1 + N_2 = 1, $$

where $m_\ell$ is the mass of the bosons in the $\ell^{th}$ component, $D \subset \mathbb{R}^2$ is the bounded physical domain of interest, $\Omega$ is the angular velocity of the condensate along the $z$-axis, $i$ is the purely imaginary unit, $u_\ell$ is the macroscopic wave function for component $\ell$ which depends on the spatial variable $\mathbf{x} = (x, y) \in \mathbb{R}^2$ and time $t > 0$. Note that $|u_\ell|^2/N_\ell$ is a real number interpreted as the probability density of measuring the particle as being at a given place. We define $\mathbf{x}^\perp = (-y, x)$, $r = |\mathbf{x}| = \sqrt{x^2 + y^2}$, and we denote by $V_\ell(r)$ the confinement function. In the case of harmonic confinement we have $V_\ell(r) = m_\ell\omega_\ell^2\frac{r^2}{2}$ with trapping frequencies $w_\ell$. The quantity $U_\ell$ defined by

$$ U_\ell = \frac{4\pi\hbar^2 a_\ell}{m_\ell}, $$

describes the interaction between the atoms of the condensate, $a_\ell$ being the scattering length of component $\ell$ which is positive for a repulsive interaction and negative for an attractive interaction and $U_{12}$ defined by

$$ U_{12} = \frac{2\pi(m_1 + m_2)\hbar^2 a_{12}}{m_1 m_2}, $$

represents the interaction between the two components. The energy functional associated to (1.1) is given by

$$ \mathcal{E}_0(u_1, u_2) = \int_D \sum_{\ell=1,2} \left( \frac{\hbar^2}{2m_\ell} |\nabla u_\ell|^2 - \hbar \Omega \Re(i u_\ell \mathbf{x}^\perp \overline{\nabla u_\ell}) + V_\ell(r) |u_\ell|^2 + \frac{1}{2} U_\ell |u_\ell|^4 \right) \\ + \int_D U_{12} |u_1|^2 |u_2|^2. \qquad (1.2) $$

<footnote n="1">¹The mass here is not related to physical masses $m_1$ nor $m_2$, instead it indicates the particles density of the normalized wave function.</footnote>
One can obtain (1.2) from (1.1) (or the opposite) using the following variational procedure (see [38])

$$i\hbar \frac{\partial u_{\ell}}{\partial t} = \frac{1}{2} \frac{\delta \mathcal{E}}{\delta u_{\ell}^{*}}. \qquad (1.3)$$

From now on, we will consider the case of harmonic confinement. Following the physical literature [19, 17, 25], let $\tilde{\omega} = \frac{\omega_1 + \omega_2}{2}$ be the average of the trapping frequencies of the two components and introduce the reduced mass $m_{12}$ such that $m_{12}^{-1} = m_1^{-1} + m_2^{-1}$. The coupled GP equations can be undimensionalised by choosing

$$\tilde{w}^{-1}, \quad \hbar\tilde{\omega} \quad \text{and} \quad r_0 = \sqrt{\hbar/(2m_{12}\tilde{\omega})} \qquad (1.4)$$

as units of time, energy and length respectively. On defining the undimensional intra-component coupling parameters $g_{\ell} = 2U_{\ell}m_{12}/\hbar^2$ and the inter-component coupling parameter $g_{12} = 2U_{12}m_{12}/\hbar^2$ and by plugging the new time and length from (1.4) into (1.1) and dividing the result by the energy unit $\hbar\tilde{\omega}$, the dimensionless coupled GP equations read

$$i\partial_t u_{\ell}(\mathbf{x}, t) = \left[ -\frac{m_{12}}{m_{\ell}}\nabla^2 - i\frac{\Omega}{\tilde{\omega}}\mathbf{x}^{\perp}\nabla + V_{\ell}(r) + g_{\ell}|u_{\ell}|^2 + g_{12}|u_{3-\ell}|^2 \right] u_{\ell}(\mathbf{x}, t), \quad (1.5)$$

where the new harmonic potential is $V_{\ell}(r) = \frac{m_{\ell}}{4m_{12}}\frac{\omega_{\ell}^2}{\tilde{\omega}^2}r^2$. The energy functional associated to (1.5) is given by

$$\mathcal{E}_1(u_1, u_2) = \int_D \sum_{\ell=1,2} \left( \frac{m_{12}}{m_\ell} |\nabla u_\ell|^2 - \Re(i u_\ell \frac{\Omega}{\tilde{\omega}} \mathbf{x}^\perp \overline{\nabla u_\ell}) + V_\ell(r) |u_\ell|^2 + \frac{1}{2} g_\ell |u_\ell|^4 \right) \\ + g_{12} |u_1|^2 |u_2|^2. \qquad (1.6)$$

One can obtain (1.6) from (1.5) (or the opposite) using the variational procedure (1.3).

In the case of strong confinement regime, a different undimensionalization method for GPE is used. We assume here that $m_1 = m_2 = m$, $\omega_1 = \omega_2 = \omega$ and $a_1 = a_2 = a$ (i.e. $U_1 = U_2 = U)^2$. Following articles [2, 6], we change the coupled GP equations to dimensionless GPE by choosing $d = \sqrt{\hbar/m\omega}$ as the characteristic length unit and

$$\varepsilon^2 = \frac{\hbar^2}{2Um}.$$

We rescale the distance by $R = \frac{d}{\sqrt{\varepsilon}}$ and define $\phi(\mathbf{r}) = Ru(\mathbf{x})$ where $\mathbf{x} = R\mathbf{r}$ and we set $\Omega_{\varepsilon} = \frac{\Omega}{\varepsilon\omega}$. The velocity $\Omega_{\varepsilon}$ is chosen such that $\Omega_{\varepsilon} < \frac{1}{\varepsilon}$, that is the trapping potential is stronger than the inertial potential. This time we choose $\hbar\omega\varepsilon$ as a unit of energy. The energy can be rewritten as

$$\mathcal{E}_2(\phi_1, \phi_2) = \int_D \sum_{\ell=1,2} \left( \frac{1}{2} |\nabla \phi_\ell|^2 - \Re(i \phi_\ell \Omega_\varepsilon \mathbf{x}^\perp \overline{\nabla \phi_\ell}) + \frac{1}{2\varepsilon^2} \tilde{V}(r) |\phi_\ell|^2 + \frac{1}{4\varepsilon^2} |\phi_\ell|^4 \right) \\ + \int_D \frac{\delta}{2\varepsilon^2} |\phi_1|^2 |\phi_2|^2,$$

<footnote n="2">²The case where the two components are the same atoms.</footnote>
where $\delta = \frac{U_{12}}{U}$ is a positive parameter determining the intra-component strength and $\tilde{V}(r) = r^2$ the confinement function in the Thomas-Fermi limit. According to the experimental values of [4, 43], $\epsilon$ is very small. This is why we are interested in the asymptotic behaviour of the minimizers of $\mathcal{E}_2$ when $\epsilon \to 0$. Notice that, since $\int |\phi_\ell|^2 = N_\ell$, adding to the energy $\mathcal{E}_2$ any function of $\int |\phi_1|^2$ and $\int |\phi_2|^2$ won't change the minimization problem.

## 3 The model

For the rest of this thesis, we will consider for a single component or two components rotating Bose-Einstein condensate in the limit of strong confinement and strong rotation, the following GPE

$$i \frac{\partial u_{\ell}}{\partial t} = -\frac{1}{2} (\nabla - i\Omega \mathbf{x}^{\perp})^2 u_{\ell} - \frac{1}{2\epsilon^2} \rho(r) u_{\ell} + \frac{1}{2\epsilon^2} |u_{\ell}|^2 u_{\ell} + \frac{\delta}{2\epsilon^2} |u_{3-\ell}|^2 u_{\ell}, \quad (1.7)$$

where $u_1, u_2 \in H^1(D \times \mathbb{R}^+, \mathbb{C})$ are two wave functions, each one representing one component of the condensate, $D \subset \mathbb{R}^2$ is the bounded physical domain of interest, $\Omega \in \mathbb{R}$ the rotation velocity, $\rho$ is a function of $r = \sqrt{x^2 + y^2}$ to be defined later, $\mathbf{x}^{\perp} = (-y, x)$, $\epsilon > 0$ such that $\frac{1}{\epsilon^2}$ measures the interaction strength in each component of the condensate (we are interested in the case where $\epsilon$ tends to 0), and $\delta > 0$ measures the interaction strength between the two components.

As studied in [6], the Gross-Pitaevskii energy functional of the rotating two components, two-dimensional BEC is given by

$$\mathcal{E}(u_1, u_2) = \sum_{\ell=1}^{2} \frac{1}{2} \int_D \|\nabla u_{\ell} - i\Omega u_{\ell} \mathbf{x}^{\perp}\|^2 dx dy + W_{\varepsilon, \delta}(u_1, u_2), \quad (1.8)$$

where $W_{\varepsilon, \delta}(u_1, u_2)$ is the confining part of the energy defined as

$$W_{\varepsilon, \delta}(u_1, u_2) = \frac{1}{4\varepsilon^2} \int_D (\rho(r) - |u_1|^2)^2 dx dy + \frac{1}{4\varepsilon^2} \int_D (\rho(r) - |u_2|^2)^2 dx dy \\ + \frac{\delta}{2\varepsilon^2} \int_D |u_1|^2 |u_2|^2 dx dy - \frac{1}{4\varepsilon^2} \int_D (\rho(r))^2 dx dy.$$

One can obtain (1.7) from (1.8) (or the opposite) using the variational procedure (1.3). The energy $\mathcal{E}$ in 1.8 will be used in all of this work. Starting from 1.9, We will change its notation to $\mathcal{E}_{\varepsilon, \delta}^{\Omega}$ to highlight the dependency on the parameters.

## 4 Aim of the thesis

Bose-Einstein condensation is a state of matter that is formed when a gas of bosons is cooled to very low temperatures. Under such conditions, a large number of bosons occupy the lowest quantum state which associates to the lowest "energy level" (also known as the ground states). In order to find the ground states, one has to find the minimum of the energy (1.8) under some constraints. My work in this thesis consist in numerically minimizing the energy (1.8) under these constraints.

With a strong confinement and depending on the rotation velocity, minimizing the energy (1.8) reveals vortices as well as different regimes of either coexistence or
segregation in a two components condensate. The aim of this thesis is to numerically simulate these vortices and to classify the other regimes as mentioned in [6]. There exists already a well implemented algorithm in MatLab called GPELab [8, 9, 11] that solves (1.7). In particular, the GPELab algorithm requires solving a linearly implicit system at each time-step, which is time consuming. Our approach consists in minimizing the energy (1.8) with an adaptive step gradient method and projection, which requires solving an explicit system at each time-step. For reference, we are also going to compare the execution time between GPELab and our approach later on.

The outline of this section is as follows. We introduce in Section 4.1 the Gross-Pitaevskii energy describing 2-dimensional, two components rotating Bose-Einstein condensates in a strong confinement regime as studied in [6]. This model is used for the computation of the minimizers for the rest of our work. In the case of one component BEC, according to [27], we can identify four regimes depending on the rotational speed. Section 4.2 is devoted to introducing the four critical rotational speeds which help us to identify the four different regimes. For the two components BEC, we distinguish two cases. The first one being the coexistence regime as developed in [3, 45]. In Section 4.3, we present the 4 different cases depending on the rotational speed and interaction strength. The second one is the segregation regime as developed in [6]. In Section 4.4, we present the segregation case both with and without rotation.

## 4.1 The Gross-Pitaevskii Energy

For the rest of this work, we consider the two components model for the energy of a two components rotating Bose-Einstein condensate in the limit of a strong confinement and a strong rotation studied in [6]

$$ \mathcal{E}_{\varepsilon, \delta}^{\Omega}(u_1, u_2) = \sum_{\ell=1}^{2} \frac{1}{2} \int_D \|\nabla u_{\ell} - i\Omega u_{\ell} \mathbf{x}^{\perp}\|^2 dx dy + W_{\varepsilon, \delta}(u_1, u_2), \quad (1.9) $$

where $D \subset \mathbb{R}^2$ is the bounded physical domain of interest (in our case $D = \{(x, y) \in \mathbb{R}^2 \text{ s.t. } x^2 + y^2 < R^2\}$), $\Omega \in \mathbb{R}$ is the rotation velocity, $\varepsilon$ and $\delta$ are positive constants, $\mathbf{x}^{\perp} = (-y, x)$, $u_1, u_2 \in H^1(D, \mathbb{C})$ are the wave functions related to each component of the condensate and $W_{\varepsilon, \delta}$ is the confining part of the energy defined as

$$ W_{\varepsilon, \delta}(u_1, u_2) = \frac{1}{4\varepsilon^2} \int_D (\rho(r) - |u_1|^2)^2 dx dy + \frac{1}{4\varepsilon^2} \int_D (\rho(r) - |u_2|^2)^2 dx dy \\ + \frac{\delta}{2\varepsilon^2} \int_D |u_1|^2 |u_2|^2 dx dy - \frac{1}{4\varepsilon^2} \int_D \rho^2(r) dx dy, $$

where $\rho$ is a function of $r = \sqrt{x^2 + y^2}$ to be defined later ($\rho \equiv 1$ in [6]), $\delta \ge 0$ measures the interaction strength between the two components, $1/\varepsilon^2$ measures the interaction strength in each component of the condensate.

As mentioned before, in order to find the lowest "energy level", we have to compute the minimizers of the energy in $H_0^1(D, \mathbb{C})$ under the constraints

$$ \int_D |u_1|^2 dx dy = M N_1, \quad \text{and} \quad \int_D |u_2|^2 dx dy = M N_2, \quad (1.10) $$

where $M = \int_D \rho(r) dx dy$, and $N_1, N_2 \ge 0$ with $N_1 + N_2 = 1$, in the regime $\varepsilon \to 0$.
**Remark 1** We can use the same two components model (1.9) for the one component simulations by taking $N_1 = 1$ and $N_2 = 0$.

Note that one has also

$$W_{\varepsilon,\delta}(u_1, u_2) = \frac{1}{4\varepsilon^2} \int_D (\rho(r) - |u_1|^2 - |u_2|^2)^2 dxdy + \frac{\delta-1}{2\varepsilon^2} \int_D |u_1|^2 |u_2|^2 dxdy. \quad (1.11)$$

The sign of the second term in (1.11) plays a role in the physical structure of the two components condensates. We can identify 3 cases:

* if $\delta = 0$ then there is no interaction between the two components of the condensate and each component is going to act like a one component condensate,

* if $0 < \delta \le 1$ then the second term in (1.11) is negative which is why both components tend to coexist (coexistence regime),

* if $\delta > 1$ then the second term in (1.11) is non-negative which is why the two components tend to avoid each other<fnref n="3" /> (segregation regime).

Observe also that for all $\ell \in \{1, 2\}$.

$$\frac{1}{2} \int_D \|\nabla u_\ell - i\Omega u_\ell \mathbf{x}^\perp\|^2 dxdy = \frac{1}{2} \int_D \|\nabla u_\ell\|^2 dxdy \quad (1.12) \\ - \Omega \int_D \Re(iu_\ell \mathbf{x}^\perp \cdot \overline{\nabla u_\ell}) dxdy + \frac{\Omega^2}{2} \int_D r^2 |u_\ell|^2 dxdy.$$

**Remark 2** Following Remark 1 of [6] and the one component analysis carried out in [27], in the regime $|\Omega| \ll 1/\varepsilon$, the contribution of the third term in (1.12) plays no role in the asymptotic $\varepsilon \to 0$. Therefore, from now on and unless stated otherwise, we shall consider in the regime $|\Omega| \ll 1/\varepsilon$ the minimization of the energy

$$E_{\varepsilon,\delta}^{\Omega}(u_1, u_2) = \sum_{\ell=1}^{2} \left[ \frac{1}{2} \int_D \|\nabla u_\ell\|^2 dxdy - \Omega \int_D \Re(iu_\ell \mathbf{x}^\perp \cdot \overline{\nabla u_\ell}) dxdy \right] + W_{\varepsilon,\delta}(u_1, u_2) \\ = \mathcal{E}_{\varepsilon,\delta}^{\Omega}(u_1, u_2) - \frac{\Omega^2}{2} \int_D r^2 (|u_1|^2 + |u_2|^2) dxdy, \quad (1.13)$$

## 4.2 One component regimes

According to [27], for one component condensates ($N_1 = 1, N_2 = 0$), the behaviour of the structure of the vortices of the minimizers of (1.13) with constraints (1.10) in the regime $\varepsilon \to 0$ depends on the dependency of the rotation speed $\Omega = \Omega_\varepsilon$ with respect to $\varepsilon$. Namely, there exist three critical rotation speeds $\Omega_\varepsilon^1 \ll \Omega_\varepsilon^2 \ll \Omega_\varepsilon^3$ such that

* if $\Omega_\varepsilon < \Omega_\varepsilon^1$, then the minimizers $u_1^\varepsilon$ have no vortices,

* if $\Omega_\varepsilon^1 < \Omega_\varepsilon < \Omega_\varepsilon^2$, then the minimizers $u_1^\varepsilon$ appear on a hexagonal lattice and are singly quantized,

<footnote n="3">³The two components of the condensate will be spatially separated.</footnote>
• if $\Omega_{\varepsilon}^2 < \Omega_{\varepsilon} < \Omega_{\varepsilon}^3$, the centrifugal force comes into play and a hole appears in the minimizers $u_1^{\varepsilon}$, the condensate looks like an annulus and the vortices appear on a lattice in the annulus,

• if $\Omega_{\varepsilon}^3 < \Omega_{\varepsilon}$, the centrifugal force is so important that the vortices retreat in the central hole of the condensate, creating a central giant vortex with high index.

Moreover, one has

$$ \Omega_{\varepsilon}^{1} \sim \log(1/\varepsilon), \quad \Omega_{\varepsilon}^{2} \sim \frac{1}{\varepsilon}, \quad \text{and} \quad \Omega_{\varepsilon}^{3} \sim \frac{1}{\varepsilon^{2} \log(1/\varepsilon)}. \qquad (1.14) $$

Note that $\delta$ plays no role in this case since $u_2 \equiv 0$ because $N_2 = 0$.

According to Remark 2, the regimes $\Omega_{\varepsilon} < \Omega_{\varepsilon}^2$ are similar for minimizers of $\mathcal{E}_{\varepsilon, \delta}^{\Omega_{\varepsilon}}$ and $E_{\varepsilon, \delta}^{\Omega_{\varepsilon}}$. In the regimes $\Omega_{\varepsilon}^2 < \Omega_{\varepsilon}$, the annulus behaviour is due to the centrifugal force missing from the energy $E_{\varepsilon, \delta}^{\Omega_{\varepsilon}}$. Numerical simulations illustrating the results of this section can be found in Chapter 3 Sections 2.1,2.2, 3.1 and 3.2.

## 4.3 Two components coexistence regimes ($0 \le \delta \le 1$)

According to [45] and [39], for $\delta \le 1$ we are in the coexistence regime since the support of each of the components of the condensate tend to overlap and has mass in the region where $\rho > 0$. In this regime $\delta \le 1$, depending on the rotational speed $\Omega$, we should observe four different regimes for the minimizers when $\varepsilon \to 0$.

• The first regime, with a very low velocity to none and $\delta \in [0, 1]$, is when the coexistence occurs with no vortices in any component of the minimizer. Moreover the profile density of each component depends on the values $N_1$ and $N_2$.

• For the second regime, when $\delta = 0$ (which means there is no interactions between the two components), we should observe, depending on the value of the rotational speed $\Omega$, the existence of a triangular vortex lattice.

• As $\delta \in (0, 1]$ increases, the positions of vortex cores in one component gradually shift from those of the other component and the triangular lattices are distorted. Eventually in the third regime, after a certain value of $\delta$, the vortices in each component form a square lattice.

• The last regime corresponds to the case $\delta \to 1$, in which we should observe either stripe or double-core vortex lattice.

Numerical simulations illustrating the results of this section can be found in Chapter 3 Sections 2.5, and 3.5.

## 4.4 Two components segregation regimes ($\delta > 1$)

In the recent paper [6], the authors consider two components condensates in the segregation case $\delta > 1$. For $N_1 \in (0, 1)$, they introduce the minimizing perimeter

$$ \ell_{N_1} = \min_{\substack{\omega \subset D \\ |\omega| = N_1}} \text{per}(\omega). $$
First, we address the minimization of (1.13) with constraints (1.10) when $\Omega = 0$. In this case, the squared modulus of the minimizers $u_1^\varepsilon$ and $u_2^\varepsilon$ tend to 1 in two separate regions of $D$ and the authors of [6] prove that their sum $v_\varepsilon^2 = |u_1^\varepsilon|^2 + |u_2^\varepsilon|^2$ and the normalized energy $\varepsilon E_{\varepsilon, \delta_\varepsilon}^0$ have the following behaviours

• in the regime $\delta_\varepsilon \varepsilon^2 \to +\infty$ (strong segregation regime [5]), $\inf_D v_\varepsilon^2$ tends to 0 and $\varepsilon E_{\varepsilon, \delta_\varepsilon}^0$ tends to some constant times $\ell_{N_1}$,

• in the regime $\varepsilon \to 0$ with $\delta > 1$ fixed, $\inf_D v_\varepsilon^2$ tends to some number between 0 and 1 and $\varepsilon E_{\varepsilon, \delta_\varepsilon}^0$ tends to some constant (which depends on $\delta$) times $\ell_{N_1}$ [32],

• in the regime $\delta_\varepsilon \to 1$ with $\tilde{\varepsilon} = \varepsilon/\sqrt{\delta_\varepsilon - 1} \to 0$, the re-scaled energy $\tilde{\varepsilon} E_{\varepsilon, \delta_\varepsilon}^0$ tends to $\ell_{N_1}/2$ and it is expected that $\inf_D v^\varepsilon$ tends to 1 [31].

Second, we address the minimization of (1.13) with constraints (1.10) when $\Omega = \Omega_\varepsilon \to +\infty$ as $\varepsilon$ tends to 0. With $\tilde{\varepsilon}$ defined as above, we consider a regime where $\delta = \delta_\varepsilon \to 1$ and $\tilde{\varepsilon} \to 0$ as $\varepsilon$ tends to 0. In [6], the authors prove that

• there exists two constants $C_1, C_2 > 0$ such that, for all $i$, if $\Omega^\varepsilon < C_i \log(1/\tilde{\varepsilon})$, then the infimum of the limiting density $|u_i|^2$ vanishes as $\varepsilon$ tends to 0;

• in the regime of moderate rotational speed $\log(1/\tilde{\varepsilon}) \ll \Omega^\varepsilon \ll 1/\tilde{\varepsilon}^2$, the limiting density $|u_i|^2$ is uniform in each region, and hence does not depend on the shape of the region.

• for higher rotational speeds $1/(\tilde{\varepsilon} \log(1/\tilde{\varepsilon})) \ll \Omega^\varepsilon \ll 1/\tilde{\varepsilon}^2$, the authors of [6] prove that the leading order in the energy $E_{\varepsilon, \delta_\varepsilon}^{\Omega^\varepsilon}$ is the vortex energy, and conjecture the possibility of observing vortex sheets.

Numerical simulations illustrating the results of this section can be found in Chapter 3 Sections 2.3, 2.4, 3.3 and 3.4. In particular, for the last bullet point, we actually observe vortex sheets (see Figures 3.36 and 3.15), validating numerically the conjecture.
# Chapter 2

## Discretization and minimization
method

Minimizing the energy $E_{\epsilon, \delta}^{\Omega}$ under the constraints (1.10) is a continuous minimization problem over an infinite dimensional manifold that we replace in this chapter by a discrete minimization problem over a finite dimensional manifold. We consider the gradient descent method with projection over the constraints manifold, which is a first-order iterative optimization algorithm for finding a local minimum of a differentiable function. By calculating the gradient of the energy $E_{\epsilon, \delta}^{\Omega}$, we are able to approach more and more towards the solution of the discrete problem after each iteration of the method. In this chapter, we describe in Section 1 the discretization of the energy using finite difference scheme as a first approximation, then in Section 2 using the fast-Fourier transform as a second approximation followed by a computation of its gradient (respectively Subsections 1.2, 1.4 and 2.3). Then in Section 3, we establish a criterion for the minimization of the energy on the finite dimensional manifold and we conclude in Section 4 by a full description of the minimization method. We develop as well in Section 4 two post processing algorithms for the detection of vortices and vortex sheets. It allows for the computation of their indices as well.

Theoretically, the function $\rho$ is strictly positive over the disk $D$. It plays no role outside the disk. The non-existence of mass outside of $D$ is modeled by the fact that $u \in H_0^1(D)$. For discretization purposes, we decide to work in a square box of length $L > D$ centered at the origin. We impose homogeneous Dirichlet boundary conditions on that box as well as negative values of (an extension) $\rho$ outside the disk of radius $R$. In the small $\epsilon$ regime, this setting should avoid the presence of mass outside the disk in the minimizers and be consistent with the continuous problem.

Now let us consider a square $[-L, L]^2$ such that $L > R$. We discretize $[-L, L]^2$ into $N + 2$ equidistant points with respect to the x-axis and $K + 2$ equidistant points with respect to the y-axis. We consider a finite sequence of points $(x_n, y_k)$, $n \in \{0, ..., N+1\}$, $k \in \{0, ..., K+1\}$ such that $x_n = -L + \frac{2L}{N+1}n$ and $y_k = -L + \frac{2L}{K+1}k$. Denoting $\delta_x = \frac{2L}{N+1}$ and $\delta_y = \frac{2L}{K+1}$, we rewrite $x_n = -L + n\delta_x$ and $y_k = -L + k\delta_y$. We use the letter $\psi$ to denote the discrete counterpart to the continuous wave function denoted by $u$ in Chapter 1 Section 4.1. We think of $\psi_{n,k}^\ell$ as an approximation
of $u_\ell(x_n, y_k)$ for $\ell = 1, 2$. In matrix form, we use the notation:

$$ \psi^\ell = \begin{pmatrix} \psi_{1,1}^\ell & \psi_{1,2}^\ell & \cdots & \psi_{1,K}^\ell \\ \psi_{2,1}^\ell & \psi_{2,2}^\ell & \cdots & \psi_{2,K}^\ell \\ \vdots & \vdots & \ddots & \vdots \\ \psi_{N,1}^\ell & \psi_{N,2}^\ell & \cdots & \psi_{N,K}^\ell \end{pmatrix}. $$

with the convention that $\psi_{n,k}^\ell = 0$ if either of these conditions is true: $n = 0$ or $k = 0$ or $n = N + 1$ or $k = K + 1$. In order to separate real and imaginary part of the unknowns, we set $\psi_{n,k}^\ell = p_{n,k}^\ell + iq_{n,k}^\ell$ for $\ell = 1, 2$. For the convenience of the reader, we will replace $\psi^\ell$ by a vector of 2 elements, the first one being the real part and the second one being the imaginary part, without changing its norm, when it is necessary.

# 1 Finite difference scheme

## 1.1 Discretization of the energy

In this section we choose the forward finite difference discretization which approaches $\nabla u(x_n, y_k)$ by $\nabla \psi_{n,k}$ defined as follows

$$ \nabla \psi_{n,k} = \left( \frac{\frac{\psi_{n+1,k} - \psi_{n,k}}{\delta_x}}{\frac{\psi_{n,k+1} - \psi_{n,k}}{\delta_y}} \right) \quad \text{for all } n \in \{0, \dots, N\} \text{ and } k \in \{0, \dots, K\}. $$

We then define $E_{\varepsilon, \delta}^{\Delta}(\psi^1, \psi^2)$ as a discrete counterpart of $E_{\varepsilon, \delta}^{\Omega}(u_1, u_2)$ by setting

$$ E_{\varepsilon, \delta}^{\Delta}(\psi^1, \psi^2) = \sum_{\ell=1,2} \left( (E_{kin})_{\varepsilon}^{\Delta}(\psi^{\ell}) + (E_r)_{\varepsilon}^{\Delta}(\psi^{\ell}) \right) + (E_W)_{\varepsilon, \delta}^{\Delta}(\psi^1, \psi^2), \quad (2.1) $$

where:

* $(E_{kin})_{\varepsilon}^{\Delta}$, corresponding to the discretization of the kinetic energy, is defined as

$$ \begin{aligned} (E_{kin})_{\varepsilon}^{\Delta}(\psi^{\ell}) &= \frac{\delta_x \delta_y}{2} \sum_{n=0}^{N} \sum_{k=0}^{K} \left( \frac{|\psi_{n+1,k}^{\ell} - \psi_{n,k}^{\ell}|^2}{\delta_x^2} + \frac{|\psi_{n,k+1}^{\ell} - \psi_{n,k}^{\ell}|^2}{\delta_y^2} \right) \\ &= \frac{\delta_x \delta_y}{2} \sum_{n=0}^{N} \sum_{k=0}^{K} \left( \frac{|p_{n+1,k}^{\ell} - p_{n,k}^{\ell}|^2}{\delta_x^2} + \frac{|p_{n,k+1}^{\ell} - p_{n,k}^{\ell}|^2}{\delta_y^2} \right) \\ &\quad + \frac{\delta_x \delta_y}{2} \sum_{n=0}^{N} \sum_{k=0}^{K} \left( \frac{|q_{n+1,k}^{\ell} - q_{n,k}^{\ell}|^2}{\delta_x^2} + \frac{|q_{n,k+1}^{\ell} - q_{n,k}^{\ell}|^2}{\delta_y^2} \right). \end{aligned} $$

* $(E_r)_{\varepsilon}^{\Delta}$, corresponding to the discretization of the rotational energy, is defined as

$$ \begin{aligned} (E_r)_{\varepsilon}^{\Delta}(\psi^{\ell}) &= -\Omega_{\varepsilon} \delta_x \delta_y \sum_{n=0}^{N} \sum_{k=0}^{K} \Re \left( i \psi_{n,k}^{\ell} [-y_k, x_n] \cdot \begin{bmatrix} \frac{\psi_{n+1,k}^{\ell} - \psi_{n,k}^{\ell}}{\delta_x} \\ \frac{\psi_{n,k+1}^{\ell} - \psi_{n,k}^{\ell}}{\delta_y} \end{bmatrix} \right) \\ &= -\Omega_{\varepsilon} \delta_x \delta_y \sum_{n=0}^{N} \sum_{k=0}^{K} \left( \frac{y_k}{\delta_x} (q_{n,k}^{\ell} p_{n+1,k}^{\ell} - q_{n+1,k}^{\ell} p_{n,k}^{\ell}) + \frac{x_n}{\delta_y} (q_{n,k+1}^{\ell} p_{n,k}^{\ell} - q_{n,k}^{\ell} p_{n,k+1}^{\ell}) \right). \end{aligned} $$
• $(E_W)_{\varepsilon, \delta}^{\Delta}$, corresponding to the discretization of the confinement and interaction energy, is defined as

$$ (E_W)_{\varepsilon, \delta}^{\Delta}(\psi^1, \psi^2) = \frac{\delta_x \delta_y}{4\varepsilon^2} \sum_{n=0}^{N} \sum_{k=0}^{K} (\rho(r_{n,k}) - (p_{n,k}^1)^2 - (q_{n,k}^1)^2 - (p_{n,k}^2)^2 - (q_{n,k}^2)^2)^2 \\ + \frac{\delta_x \delta_y (\delta - 1)}{2\varepsilon^2} \sum_{n=0}^{N} \sum_{k=0}^{K} ((p_{n,k}^1)^2 + (q_{n,k}^1)^2) ((p_{n,k}^2)^2 + (q_{n,k}^2)^2), $$
where $r_{n,k} = \sqrt{x_n^2 + y_k^2}$.

**Remark 3** We define $(E_{cf})_{\varepsilon, \delta}^{\Delta}$, corresponding to the discretization of the centrifugal energy (see Remark 2), as

$$ (E_{cf})_{\varepsilon, \delta}^{\Delta}(\psi^1, \psi^2) = \frac{\delta_x \delta_y \Omega^2}{2} \sum_{n=0}^{N} \sum_{k=0}^{K} r_{n,k}^2 ((p_{n,k}^1)^2 + (q_{n,k}^1)^2 + (p_{n,k}^2)^2 + (q_{n,k}^2)^2), $$

such as $\mathcal{E}_{\varepsilon, \delta}^{\Delta}(\psi^1, \psi^2) = E_{\varepsilon, \delta}^{\Delta}(\psi^1, \psi^2) + (E_{cf})_{\varepsilon, \delta}^{\Delta}(\psi^1, \psi^2)$ is the discrete counter part of $\mathcal{E}_{\varepsilon, \delta}^{\Omega}$.

## 1.2 Computation of the gradient of the discrete energy (2.1)

Let $P^l, Q^l \in \mathbb{R}^{N \times K}$ be the real and imaginary part (respectively) of $\psi^l$ for $l = 1, 2$. In order to minimize the energy $E_{\varepsilon, \delta}^{\Delta}$ defined in (2.1) using the gradient method described later in Section 4, we have to compute at each iteration the gradient of the energy $E_{\varepsilon, \delta}^{\Delta}$ with respect to $P^1, Q^1, P^2, Q^2$. Therefore, let us express the gradient of each energy term (still using the convention $p_{n,k}^l = q_{n,k}^l = 0$ if one of these conditions is true: $n = 0$ or $k = 0$ or $n = N + 1$ or $k = K + 1$).

**Proposition 4** For all $n \in \{1, \dots, N\}$ and $k \in \{1, \dots, K\}$, and for $l \in \{1, 2\}$, the gradient of the energy $E_{\varepsilon, \delta}^{\Delta}$ is equal to

• $\frac{\partial(E_{kin})_{\varepsilon}^{\Delta}(\psi^l)}{\partial p_{n,k}^{l}} = \frac{\delta_y}{\delta_x}(2p_{n,k}^{l} - p_{n+1,k}^{l} - p_{n-1,k}^{l}) + \frac{\delta_x}{\delta_y}(2p_{n,k}^{l} - p_{n,k-1}^{l} - p_{n,k+1}^{l})$,

• $\frac{\partial(E_{kin})_{\varepsilon}^{\Delta}(\psi^l)}{\partial q_{n,k}^{l}} = \frac{\delta_y}{\delta_x}(2q_{n,k}^{l} - q_{n+1,k}^{l} - q_{n-1,k}^{l}) + \frac{\delta_x}{\delta_y}(2q_{n,k}^{l} - q_{n,k-1}^{l} - q_{n,k+1}^{l})$,

• $\frac{\partial(E_r)_{\varepsilon}^{\Delta}(\psi^l)}{\partial p_{n,k}^{l}} = -\Omega_{\varepsilon} \delta_x \delta_y \left( \frac{y_k}{\delta_x}(q_{n-1,k}^{l} - q_{n+1,k}^{l}) + \frac{x_n}{\delta_y}(q_{n,k+1}^{l} - q_{n,k-1}^{l}) \right)$,

• $\frac{\partial(E_r)_{\varepsilon}^{\Delta}(\psi^l)}{\partial q_{n,k}^{l}} = -\Omega_{\varepsilon} \delta_x \delta_y \left( \frac{y_k}{\delta_x}(p_{n+1,k}^{l} - p_{n-1,k}^{l}) + \frac{x_n}{\delta_y}(p_{n,k-1}^{l} - p_{n,k+1}^{l}) \right)$,

• $\frac{\partial(E_W)_{\varepsilon, \delta}^{\Delta}(\psi^1, \psi^2)}{\partial p_{n,k}^{l}} = -\frac{\delta_x \delta_y}{\varepsilon^2} p_{n,k}^{l} (\rho(r_{n,k}) - |\psi_{n,k}^1|^2 - |\psi_{n,k}^2|^2) + \frac{\delta_x \delta_y (\delta - 1)}{\varepsilon^2} p_{n,k}^{l} ((p_{n,k}^{3-l})^2 + (q_{n,k}^{3-l})^2)$,
• $\frac{\partial(E_W)^\Delta_{\varepsilon,\delta}(\psi^1, \psi^2)}{\partial q_{n,k}^\ell} = -\frac{\delta_x \delta_y}{\varepsilon^2} q_{n,k}^\ell (\rho(r_{n,k}) - |\psi_{n,k}^1|^2 - |\psi_{n,k}^2|^2) + \frac{\delta_x \delta_y (\delta - 1)}{\varepsilon^2} q_{n,k}^\ell ((p_{n,k}^{3-\ell})^2 + (q_{n,k}^{3-\ell})^2).$

• $\frac{\partial(E_{cf})^\Delta_{\varepsilon,\delta}(\psi^1, \psi^2)}{\partial p_{n,k}^\ell} = \delta_x \delta_y \Omega^2 r_{n,k}^2 p_{n,k}^\ell$.

• $\frac{\partial(E_{cf})^\Delta_{\varepsilon,\delta}(\psi^1, \psi^2)}{\partial q_{n,k}^\ell} = \delta_x \delta_y \Omega^2 r_{n,k}^2 q_{n,k}^\ell$.

## 1.3 Discretization of the energy in a matrix form

In this subsection, we rewrite the discretized energy as a combination of sums and products of matrices.

From now on, we consider the case where $N = K$, in particular $\delta_x = \delta_y$. Let $\psi^\ell = P^\ell + iQ^\ell$ where $P^\ell, Q^\ell \in \mathbb{R}^{N^2}$ for $\ell = 1, 2$ and let the multiplication symbol $*$ denote the element wise multiplication.

**Example 5** Let $P = \begin{pmatrix} p_1 \\ p_2 \\ \vdots \\ p_N \end{pmatrix} \in \mathbb{R}^N$, and let $Q = \begin{pmatrix} q_1 \\ q_2 \\ \vdots \\ q_N \end{pmatrix} \in \mathbb{R}^N$, then

$$P * Q = \begin{pmatrix} p_1 \\ p_2 \\ \vdots \\ p_N \end{pmatrix} * \begin{pmatrix} q_1 \\ q_2 \\ \vdots \\ q_N \end{pmatrix} = \begin{pmatrix} p_1 q_1 \\ p_2 q_2 \\ \vdots \\ p_N q_N \end{pmatrix}, \quad \text{and} \quad (P)^2 = P * P.$$

We define the following matrices:

• Let $A_k \in \mathbb{R}^{N^2 \times N^2}$ such as:

$$A_k = \begin{pmatrix} \tilde{A} & -I & 0 & \cdots & \cdots & 0 \\ -I & \tilde{A} & -I & \cdots & \cdots & \vdots \\ 0 & -I & \tilde{A} & -I & \cdots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & \ddots & \tilde{A} \\ 0 & \cdots & \cdots & 0 & -I & \tilde{A} \end{pmatrix}$$

with $I, \tilde{A} \in \mathbb{R}^{N \times N}$, $I$ being the identity matrix and

$$\tilde{A} = \begin{pmatrix} 4 & -1 & 0 & \cdots & \cdots & 0 \\ -1 & 4 & -1 & \cdots & \cdots & \vdots \\ 0 & -1 & 4 & -1 & \cdots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & \ddots & \tilde{A} \\ 0 & \cdots & \cdots & 0 & -1 & 4 \end{pmatrix}.$$
• Let $\tilde{X}_n, \tilde{Y} \in \mathbb{R}^{N \times N}$ be defined for all $n \in \{1, \dots, N\}$ as

$$ \tilde{X}_n = \begin{pmatrix} 0 & -x_n & 0 & \cdots & 0 \\ x_n & 0 & -x_n & \cdots & \vdots \\ 0 & x_n & 0 & -x_n & \vdots \\ \vdots & \ddots & \ddots & \ddots & 0 \\ \vdots & \ddots & x_n & 0 & -x_n \\ 0 & \cdots & \cdots & 0 & x_n \end{pmatrix}, \quad \tilde{Y} = \begin{pmatrix} y_1 & 0 & 0 & \cdots & 0 \\ 0 & y_2 & 0 & \cdots & \vdots \\ 0 & 0 & y_3 & 0 & \vdots \\ \vdots & \ddots & \ddots & \ddots & 0 \\ \vdots & \ddots & 0 & y_{N-1} & 0 \\ 0 & \cdots & \cdots & 0 & y_N \end{pmatrix} $$

where $x_n = -L + n\delta_x$ and $y_n = -L + n\delta_y$ for all $n \in \{1, \dots, N\}$. And let $A_r \in \mathbb{R}^{N^2 \times N^2}$ be defined as:

$$ A_r = \begin{pmatrix} \tilde{X}_1 & \tilde{Y} & 0 & \cdots & 0 \\ -\tilde{Y} & \tilde{X}_2 & \tilde{Y} & \cdots & \vdots \\ 0 & -\tilde{Y} & \tilde{X}_3 & \tilde{Y} & \vdots \\ \vdots & \ddots & \ddots & \ddots & 0 \\ \vdots & \ddots & -\tilde{Y} & \tilde{X}_{N-1} & \tilde{Y} \\ 0 & \cdots & 0 & -\tilde{Y} & \tilde{X}_N \end{pmatrix} $$

• Finally let $V_W \in \mathbb{R}^{N^2}$, such as:

$$ V_W = \begin{pmatrix} \rho(r_{1,1}) \\ \rho(r_{1,2}) \\ \vdots \\ \rho(r_{2,1}) \\ \vdots \\ \rho(r_{N,K}) \end{pmatrix}. $$

Let $\psi \in \mathbb{C}^{(N+2)^2}$. We define the following norm

$$ \|\psi\|_{\Delta}^{2} = \delta_{x}^{2} \sum_{n=0}^{N+1} \sum_{k=0}^{N+1} |\psi_{n,k}|^{2}. \qquad (2.2) $$

**Proposition 6** *Using the matrices $A_k$, $A_r$, $V_W$ and the expressions obtained in (2.1)*
we rewrite the discretized energy in matrix form as follows

$$
\begin{aligned}
E_{\varepsilon, \delta}^{\Delta}(\psi^1, \psi^2) = & \frac{1}{2} \sum_{\ell=1}^{2} \sum_{\substack{n=0 \\ k=0}}^{N} (|p_{n+1,k}^{\ell} - p_{n,k}^{\ell}|^2 + |p_{n,k+1}^{\ell} - p_{n,k}^{\ell}|^2) \\
& + \frac{1}{2} \sum_{\ell=1}^{2} \sum_{\substack{n=0 \\ k=0}}^{N} (|q_{n+1,k}^{\ell} - q_{n,k}^{\ell}|^2 + |q_{n,k+1}^{\ell} - q_{n,k}^{\ell}|^2) \\
& - \Omega_{\varepsilon} \delta_x \sum_{\ell=1}^{2} \sum_{\substack{n=0 \\ k=0}}^{N} (y_k(q_{n,k}^{\ell} p_{n+1,k}^{\ell} - q_{n+1,k}^{\ell} p_{n,k}^{\ell}) + x_n(q_{n,k+1}^{\ell} p_{n,k}^{\ell} - q_{n,k}^{\ell} p_{n,k+1}^{\ell})) \\
& + \delta_x^2 \sum_{\substack{n=1 \\ k=1}}^{N} \left( \frac{1}{4\varepsilon^2} \left( \rho(r_{n,k}) - |\psi_{n,k}^1|^2 - |\psi_{n,k}^2|^2 \right)^2 + \frac{\delta-1}{2\varepsilon^2} |\psi_{n,k}^1|^2 |\psi_{n,k}^2|^2 \right) \\
= & \sum_{\ell=1}^{2} \left( \frac{1}{2}(P^{\ell})^t A_k P^{\ell} + \frac{1}{2}(Q^{\ell})^t A_k Q^{\ell} - \Omega_{\varepsilon} \delta_x (Q^{\ell})^t A_r P^{\ell} \right) \\
& + \frac{1}{4\varepsilon^2} \left\| V_W - \sum_{\ell=1}^{2} ((P^{\ell})^2 - (Q^{\ell})^2) \right\|_{\Delta}^2 \\
& + \frac{\delta-1}{2\varepsilon^2} \left\| ((P^1)^2 + (Q^1)^2) * ((P^2)^2 + (Q^2)^2) \right\|_{\Delta}
\end{aligned}
$$

**Remark 7** Following Remark 3, we can also rewrite the centrifugal energy in a matrix form

$$ (E_{cf})_{\varepsilon, \delta}^{\Delta}(\psi^1, \psi^2) = \frac{\delta_x \delta_y \Omega^2}{2} \|R_{cf} * (|\psi^1|^2 + |\psi^2|^2)\|_1, $$

with

$$ R_{cf} = \begin{pmatrix} x_1^2 + y_1^2 & x_1^2 + y_2^2 & \cdots & x_1^2 + y_N^2 \\ x_2^2 + y_1^2 & x_2^2 + y_2^2 & \cdots & x_2^2 + y_N^2 \\ \vdots & \vdots & \ddots & \vdots \\ x_N^2 + y_1^2 & x_N^2 + y_2^2 & \cdots & x_N^2 + y_N^2 \end{pmatrix}. $$

## 1.4 Computation of the gradient of the discrete energy in a matrix form

Still assuming $N = K$, we have the following proposition

**Proposition 8** The gradient in the matrix form is equal to

$$ \nabla_{P^1,Q^1} E_{\varepsilon,\delta}^{\Delta} = \begin{pmatrix} A_k P^1 - \Omega_{\varepsilon} \delta_x (Q^1)^t A_r - \frac{\delta_x^2}{\varepsilon^2} (V_W - \sum_{\ell=1}^2 ((P^\ell)^2 - (Q^\ell)^2)) * P^1 + \frac{\delta_x^2 (\delta-1)}{\varepsilon^2} P^1 * ((P^2)^2 + (Q^2)^2) \\ A_k Q^1 - \Omega_{\varepsilon} \delta_x (P^1)^t A_r - \frac{\delta_x^2}{\varepsilon^2} (V_W - \sum_{\ell=1}^2 ((P^\ell)^2 - (Q^\ell)^2)) * Q^1 + \frac{\delta_x^2 (\delta-1)}{\varepsilon^2} Q^1 * ((P^2)^2 + (Q^2)^2) \end{pmatrix}, $$

$$ \nabla_{P^2,Q^2} E_{\varepsilon,\delta}^{\Delta} = \begin{pmatrix} A_k P^2 - \Omega_{\varepsilon} \delta_x (Q^2)^t A_r - \frac{\delta_x^2}{\varepsilon^2} (V_W - \sum_{\ell=1}^2 ((P^\ell)^2 - (Q^\ell)^2)) * P^2 + \frac{\delta_x^2 (\delta-1)}{\varepsilon^2} P^2 * ((P^1)^2 + (Q^1)^2) \\ A_k Q^2 - \Omega_{\varepsilon} \delta_x (P^2)^t A_r - \frac{\delta_x^2}{\varepsilon^2} (V_W - \sum_{\ell=1}^2 ((P^\ell)^2 - (Q^\ell)^2)) * Q^2 + \frac{\delta_x^2 (\delta-1)}{\varepsilon^2} Q^2 * ((P^1)^2 + (Q^1)^2) \end{pmatrix}. $$
**Remark 9** Note that here the vectors $P^l$, $Q^l$ and $V_W$ are in $\mathbb{R}^{N^2}$, and the matrices $A_r$, $A_k$ are in $\mathbb{R}^{N^2 \times N^2}$ so that $\nabla_{P^l,Q^l} E_{\varepsilon,\delta}^{\Delta}$ is in $\mathbb{R}^{2N^2}$.

# 2 Fast Fourier Transform scheme

Using the finite difference scheme, we noticed the frequent appearance of a square shaped behaviour in the simulations when $\Omega$ is big enough (see Chapter 3, Section 2, the one component case). So we decided to use the discrete Fourier transform and its inverse for the discretization of the terms in the energy (1.13) that involve gradients (mainly the kinetic and rotational energies). This choice allows the use of Fast Fourier Transform algorithms for the computation of the gradient of the discrete energy (see Section 2.3). To do so, we set the following definitions, which are related to that of Python Numpy.

## 2.1 Fast Fourier transform discretization

Let $v \in \mathbb{C}^{(N+2) \times (K+2)}$. The discrete Fourier transforms of $v$ in the $x$- and $y$-direction are respectively given by:

$$ 
\begin{aligned}
\hat{v}_{n,k}^{x}(\xi_p) &= \frac{\delta_x}{\sqrt{2\pi}} \sum_{n=0}^{N+1} v_{n,k} e^{-ix_n \xi_p} \\
&= \frac{\delta_x}{\sqrt{2\pi}} e^{-\frac{i\pi(N+1)}{2}} e^{ip\pi} \sum_{n=0}^{N+1} v_{n,k} e^{-2\pi i \frac{np}{N+1}} e^{i\pi n},
\end{aligned}
$$

for $p \in \{0, \dots, N+1\}$ and

$$ 
\begin{aligned}
\hat{v}_{n,k}^{y}(\lambda_q) &= \frac{\delta_y}{\sqrt{2\pi}} \sum_{k=0}^{K+1} v_{n,k} e^{-iy_k \lambda_q} \\
&= \frac{\delta_y}{\sqrt{2\pi}} e^{-\frac{i\pi(K+1)}{2}} e^{iq\pi} \sum_{k=0}^{K+1} v_{n,k} e^{-2\pi i \frac{kq}{K+1}} e^{i\pi k},
\end{aligned}
$$

for $q \in \{0, \dots, K+1\}$, with

$$ \xi_p = -\frac{\pi(N+1)}{2L} + p\delta_{\xi}, \quad \lambda_q = -\frac{\pi(K+1)}{2L} + q\delta_{\lambda}, \quad \delta_{\xi} = \delta_{\lambda} = \frac{\pi}{L}. \quad (2.3) $$

The inverse discrete Fourier transforms in the $\xi$- and $\lambda$-direction are respectively given by:

$$ 
\begin{aligned}
\tilde{v}_{p,q}^{\xi}(x_n) &= \frac{\delta_{\xi}}{\sqrt{2\pi}} \sum_{p=0}^{N+1} v_{p,q} e^{ix_n \xi_p} \\
&= \frac{\delta_{\xi}}{\sqrt{2\pi}} e^{\frac{i\pi(N+1)}{2}} e^{-in\pi} \sum_{p=0}^{N+1} v_{p,q} e^{2\pi i \frac{np}{N+1}} e^{-i\pi p}
\end{aligned}
$$

$$ 
\begin{aligned}
\tilde{v}_{p,q}^{\lambda}(y_k) &= \frac{\delta_{\lambda}}{\sqrt{2\pi}} \sum_{q=0}^{K+1} v_{p,q} e^{iy_k \lambda_q} \\
&= \frac{\delta_{\lambda}}{\sqrt{2\pi}} e^{\frac{i\pi(K+1)}{2}} e^{-ik\pi} \sum_{q=0}^{K+1} v_{p,q} e^{2\pi i \frac{kq}{K+1}} e^{-i\pi q}.
\end{aligned}
$$
Notice that we have $\delta_x\delta_\xi = \frac{2\pi}{N+1}$ and $\delta_y\delta_\lambda = \frac{2\pi}{K+1}$ and, for all $v$, we have $\tilde{v}^x = \tilde{v}^y = \frac{N+2}{N+1}v$. From now on we consider the case where $N = K$, in particular $\delta_x = \delta_y$.

Since we always use discrete Fourier transform and its inverse to compute the quantities of interest, we decide to set the following definitions, which correspond to that of Python Numpy and allow for the computation via the Fast Fourier Transform algorithm:

$$ \bullet \ (\mathbf{fft}_x(v))_{p,k} = \sum_{m=0}^{N+1} v_{m,k} e^{i\pi m} e^{-2\pi i \frac{mp}{N+1}}, $$

$$ \bullet \ (\mathbf{ifft}_x(v))_{n,k} = \frac{1}{N+1} e^{-i\pi n} \sum_{p=0}^{N+1} v_{p,k} e^{2\pi i \frac{np}{N+1}}, $$

$$ \bullet \ (\mathbf{fft}_y(v))_{n,q} = \sum_{l=0}^{N+1} v_{n,l} e^{i\pi l} e^{-2\pi i \frac{lq}{N+1}}, $$

$$ \bullet \ (\mathbf{ifft}_y(v))_{n,k} = \frac{1}{N+1} e^{-i\pi k} \sum_{q=0}^{N+1} v_{n,q} e^{2\pi i \frac{kq}{N+1}}. $$

Notice that we also have $\mathbf{ifft}_x(\mathbf{fft}_x(v)) = \mathbf{ifft}_y(\mathbf{fft}_y(v)) = \frac{N+2}{N+1}v$.

Using these definitions, we can now discretize the energy $E_{\varepsilon,\delta}^\Omega$.

## 2.2 Discretization of the energy

In this section, we will consider the same approximation as Section 2. Here, the borders are included so that:

$$ \psi^\ell = \begin{pmatrix} \psi_{0,0}^\ell & \psi_{0,1}^\ell & \cdots & \psi_{0,N+1}^\ell \\ \psi_{1,0}^\ell & \psi_{1,1}^\ell & \cdots & \psi_{1,N+1}^\ell \\ \vdots & \vdots & \ddots & \vdots \\ \psi_{N+1,0}^\ell & \psi_{N+1,1}^\ell & \cdots & \psi_{N+1,N+1}^\ell \end{pmatrix} \in \mathbb{C}^{(N+2)^2}, $$

with the convention that $\psi_{n,k}^\ell = 0$ if one of these conditions is true: $n = 0$ or $k = 0$ or $n = N + 1$ or $k = N + 1$.

Let us define for all $n, k \in \{0, \cdots, N + 1\}, X, Y, \Xi, \Lambda \in \mathbb{R}^{(N+2)^2}$ by

$$ X_{n,k} = x_n, \quad Y_{n,k} = y_k, \quad \Xi_{n,k} = \xi_n, \quad \Lambda_{n,k} = \lambda_k. $$

These will be used for the discretization of the rotational and kinetic energy. We still denote by $*$ the element wise multiplication operator between 2 matrices (see Example 5).

We define $E_{\varepsilon,\delta}^\Delta(\psi^1, \psi^2)$ as a discrete counterpart of $E_{\varepsilon,\delta}^\Omega(u_1, u_2)$ by setting for $\psi^1, \psi^2 \in \mathbb{C}^{(N+2)^2}$

$$ E_{\varepsilon,\delta}^{\Delta}(\psi^1, \psi^2) = \sum_{\ell=1,2} ((E_{kin})_{\varepsilon}^{\Delta}(\psi^{\ell}) + (E_{r})_{\varepsilon}^{\Delta}(\psi^{\ell})) + (E_{W})_{\varepsilon,\delta}^{\Delta}(\psi^1, \psi^2), \quad (2.4) $$

where:

* $(E_{kin})_\varepsilon^\Delta$, corresponding to the discretization of the kinetic energy, is defined as

$$ (E_{kin})_{\varepsilon}^{\Delta}(\psi^{\ell}) = \frac{\delta_x^2}{2} \sum_{n,k=0}^{N+1} \left( \left| \mathbf{ifft}_x(i\Xi * \mathbf{fft}_x(\psi^{\ell})) \right|^2 + \left| \mathbf{ifft}_y(i\Lambda * \mathbf{fft}_y(\psi^{\ell})) \right|^2 \right)_{n,k} \quad (2.5) $$
• $(E_r)_{\varepsilon}^{\Delta}$, corresponding to the discretization of the rotational energy, is defined as

$$ (E_r)_{\varepsilon}^{\Delta}(\psi^{\ell}) = -\Omega_{\varepsilon}\delta_x^2 \sum_{n,k=0}^{N+1} \Re \left( -i\overline{\psi^{\ell}} * \left[ -Y * \mathbf{ifft}_x(i\Xi * \mathbf{fft}_x(\psi^{\ell})) + X * \mathbf{ifft}_y(i\Lambda * \mathbf{fft}_y(\psi^{\ell})) \right] \right)_{n,k} $$

• $(E_W)_{\varepsilon,\delta}^{\Delta}$, corresponding to the discretization of the confinement energy, is defined as

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

where $r_{n,k} = \sqrt{x_n^2 + y_k^2}$.

**Remark 10** Similarly to Remark 3 let $(E_{cf})_{\varepsilon,\delta}^{\Delta}$, corresponding to the discretization of the centrifugal energy, be defined as

$$ (E_{cf})_{\varepsilon,\delta}^{\Delta}(\psi^1, \psi^2) = \frac{\delta_x \delta_y \Omega^2}{2} \sum_{n,k=0}^{N+1} \left[ (X^2 + Y^2) * (|\psi^1|^2 + |\psi^2|^2) \right]_{n,k}. $$

In this case, $\mathcal{E}_{\varepsilon,\delta}^{\Delta}(\psi^1, \psi^2) = E_{\varepsilon,\delta}^{\Delta}(\psi^1, \psi^2) + (E_{cf})_{\varepsilon,\delta}^{\Delta}(\psi^1, \psi^2)$ is the discrete counter part of $\mathcal{E}_{\varepsilon,\delta}^{\Omega}(2)$.

## 2.3 Computation of the gradient of the discrete energy $E_{\varepsilon,\delta}^{\Delta}$

Each energy term in (2.4) is a function of the $4N^2$ real variables $(P_{n,k}^1)_{1\le n,k\le N}$, $(Q_{n,k}^1)_{1\le n,k\le N}$, $(P_{n,k}^2)_{1\le n,k\le N}$, $(Q_{n,k}^2)_{1\le n,k\le N}$. Let us compute the gradient of each energy term (still using the convention $\psi_{n,k} = 0$ if one of these conditions is true: $n = 0$ or $k = 0$ or $n = N + 1$ or $k = N + 1$) with respect to these variables. Let $\tilde{V}_W \in \mathbb{R}^{(N+2)^2}$ be defined as:

$$ \tilde{V}_W = \begin{pmatrix} \rho(r_{0,0}) & \rho(r_{0,1}) & \cdots & \rho(r_{0,N+1}) \\ \rho(r_{1,0}) & \rho(r_{1,1}) & \cdots & \rho(r_{1,N+1}) \\ \vdots & \vdots & \ddots & \vdots \\ \rho(r_{N+1,0}) & \rho(r_{N+1,1}) & \cdots & \rho(r_{N+1,N+1}) \end{pmatrix}. $$

**Proposition 11** For $\ell \in \{1, 2\}$ we have for the discrete kinetic energy

$$ \frac{\partial E_{\text{kin}}^{\Delta}}{\partial P^{\ell}}(\psi^{\ell}) = \delta_x^2 \Re \left( \left[ \mathbf{ifft}_x (\Xi^2 \mathbf{fft}_x(\psi^{\ell})) \right]_{n,k} + \left[ \mathbf{ifft}_y (\Lambda^2 \mathbf{fft}_y(\psi^{\ell})) \right]_{n,k} \right)_{1 \le n,k \le N}, $$

$$ \frac{\partial E_{\text{kin}}^{\Delta}}{\partial Q^{\ell}}(\psi^{\ell}) = \delta_x^2 \Re \left( -i \left[ \mathbf{ifft}_x (\Xi^2 \mathbf{fft}_x(\psi^{\ell})) \right]_{n,k} - i \left[ \mathbf{ifft}_y (\Lambda^2 \mathbf{fft}_y(\psi^{\ell})) \right]_{n,k} \right)_{1 \le n,k \le N}. $$
**Proposition 12** For $\ell \in \{1, 2\}$ we have for the discrete rotational energy

$$ \frac{\partial(E_r)_\varepsilon^\Delta}{\partial P^\ell}(\psi^\ell) = -2\Omega_\varepsilon \delta_x^2 \Re \left( [X * \mathbf{i}fft_y (\Lambda * \mathbf{fft}_y(\psi^\ell)) - Y * \mathbf{i}fft_x (\Xi * \mathbf{fft}_x(\psi^\ell))]_{n,k} \right)_{1 \le n, k \le N}, \\ \frac{\partial(E_r)_\varepsilon^\Delta}{\partial Q^\ell}(\psi^\ell) = -2\Omega_\varepsilon \delta_x^2 \Re \left( -i [X * \mathbf{i}fft_y (\Lambda * \mathbf{fft}_y(\psi^\ell)) - Y * \mathbf{i}fft_x (\Xi * \mathbf{fft}_x(\psi^\ell))]_{n,k} \right)_{1 \le n, k \le N}. $$

**Proposition 13** For $\ell \in \{1, 2\}$ we have for the discrete interaction energy

$$ \frac{\partial(E_W)_\varepsilon^\Delta}{\partial P^\ell}(\psi^1, \psi^2) = -\frac{\delta_x^2}{\varepsilon^2} \left( \left[ P^\ell * \left( \tilde{V}_W - \sum_{j=1}^2 |\psi^j|^2 \right) - (\delta - 1) P^\ell * |\psi^{3-\ell}|^2 \right]_{n,k} \right)_{1 \le n, k \le N}, \\ \frac{\partial(E_W)_\varepsilon^\Delta}{\partial Q^\ell}(\psi^1, \psi^2) = -\frac{\delta_x^2}{\varepsilon^2} \left( \left[ Q^\ell * \left( \tilde{V}_W - \sum_{j=1}^2 |\psi^j|^2 \right) - (\delta - 1) Q^\ell * |\psi^{3-\ell}|^2 \right]_{n,k} \right)_{1 \le n, k \le N}. $$

**Proposition 14** For $\ell \in \{1, 2\}$ we have for the discrete centrifugal energy

$$ \frac{\partial(E_{cf})_{\varepsilon,\delta}^{\Delta}}{\partial P^{\ell}}(\psi^1, \psi^2) = \delta_x^2 \Omega^2 \left( \left[ (X^2 + Y^2) * P^{\ell} \right]_{n,k} \right)_{1 \le n, k \le N}, \\ \frac{\partial(E_{cf})_{\varepsilon,\delta}^{\Delta}}{\partial Q^{\ell}}(\psi^1, \psi^2) = \delta_x^2 \Omega^2 \left( \left[ (X^2 + Y^2) * Q^{\ell} \right]_{n,k} \right)_{1 \le n, k \le N}. $$

## 3 A criterion for the minimization of $E_{\varepsilon,\delta}^{\Delta}$ under constraints

We replace the continuous constraints $\|u^\ell\|_2^2 = \int_D |u^\ell|^2 dx dy = MN_\ell$ ($\ell = 1, 2$) with the discrete analogues

$$ \|\psi^\ell\|_{\Delta}^2 = \delta_x \delta_y \sum_{i=1}^{N} \sum_{j=1}^{N} |\psi_{i,j}^\ell|^2 = MN_\ell, \quad (\ell = 1, 2). \qquad (2.6) $$

In order to derive a criterion for the minimizers of $E_{\varepsilon,\delta}^{\Delta}$ ((2.1) and (2.4)) under the constraints (2.6), we introduce the following function which evaluates the energy on the constraints manifold defined by (2.6) close to some $(\psi^{1*}, \psi^{2*})$:

$$ f^{\Delta}(t) = E_{\varepsilon,\delta}^{\Delta} \left( \frac{\psi^{1*} + t\psi^1}{\|\psi^{1*} + t\psi^1\|_{\Delta}} \sqrt{N_1 M}, \frac{\psi^{2*} + t\psi^2}{\|\psi^{2*} + t\psi^2\|_{\Delta}} \sqrt{N_2 M} \right), \quad (2.7) $$

where $t \in \mathbb{R}$ is a real variable that later on will tend to 0, where $\psi^{1*}, \psi^{2*} \in \mathbb{C}^{(N+2)^2}$ are minimizers of $E_{\varepsilon,\delta}^{\Delta}$ under the constraints (2.6) satisfying the homogeneous Dirichlet boundary conditions and where $\psi^1, \psi^2 \in \mathbb{C}^{(N+2)^2}$ are arbitrary and satisfy the homogeneous Dirichlet boundary conditions. In the following, we will denote $\psi^{\ell*} = P^{\ell*} + iQ^{\ell*}$ and $\psi^{\ell} = P^{\ell} + iQ^{\ell}$ for $\ell \in \{1, 2\}$.

We identify $\mathbb{C}^{2N^2}$ with the subspace of $\mathcal{M}_{N+2}(\mathbb{C})^2$ consisting in pairs of matrices with zeros first and last row and column. We have the following proposition:
**Proposition 15** If $(\psi^{1*}, \psi^{2*})$ is a minimizer of $E_{\varepsilon, \delta}^{\Delta}$ under the constraints (2.6), then $K^{\Delta} = 0$ with $K^{\Delta}$ given by

$$ K^{\Delta} = \sum_{\ell=1}^{2} \left\| \nabla E_{\varepsilon, \delta}^{\Delta} (\psi^{\ell*}) - \frac{\delta_x^2}{\|\psi^{\ell*}\|_{\Delta}^2} \sum_{n,k=0}^{N+1} \left( \frac{\partial E_{\varepsilon, \delta}^{\Delta}}{\partial P^{\ell}} (\psi^{\ell*}) * P^{\ell*} + \frac{\partial E_{\varepsilon, \delta}^{\Delta}}{\partial Q^{\ell}} (\psi^{\ell*}) * Q^{\ell*} \right) \right\|_{\Delta} \quad (2.8) $$

**Proof:** Let us denote the scalar product on $\mathbb{C}^{2N^2}$ by $\langle u, v \rangle_{\Delta} = \Re \left( \sum_{k,l=1}^{N} u_{k,l} \bar{v}_{k,l} \right)$ with the following relation $\delta_x^2 \langle u, u \rangle_{\Delta} = \|u\|_{\Delta}^2$. The fact that $(\psi^{1*}, \psi^{2*})$ are minimizers of $E_{\varepsilon, \delta}^{\Delta}$ under the constraints (2.6) implies that

$$ (f^{\Delta})'(0) = 0, \quad (2.9) $$

for all $(\psi^1, \psi^2)$. Using (2.7), we set for $\ell = 1, 2$

$$ f_{\ell,P}^{\Delta}(t) = \frac{P^{\ell*} + tP^{\ell}}{\|\psi^{\ell*} + t\psi^{\ell}\|_{\Delta}} \sqrt{N_{\ell}M}, \quad \text{and} \quad f_{\ell,Q}^{\Delta}(t) = \frac{Q^{\ell*} + tQ^{\ell}}{\|\psi^{\ell*} + t\psi^{\ell}\|_{\Delta}} \sqrt{N_{\ell}M}, $$

so that, using (2.6), we have

$$ \|f_{\ell,P}^{\Delta}(t) + if_{\ell,Q}^{\Delta}(t)\|_{\Delta}^2 = N_{\ell}M. $$

We also set, in order to compute the first derivative of $f^{\Delta}$ at $t = 0$,

$$ a_{\ell} = \frac{\delta_x^2}{\|\psi^{\ell*}\|_{\Delta}^2} \left( \sum_{n,k=1}^{N} (p_{n,k}^{\ell*} p_{n,k}^{\ell} + q_{n,k}^{\ell*} q_{n,k}^{\ell}) \right) = \frac{\delta_x^2}{\|\psi^{\ell*}\|_{\Delta}^2} \sum_{n,k=1}^{N} \left[ P^{\ell*} * P^{\ell} + Q^{\ell*} * Q^{\ell} \right]_{n,k}, $$

so that we have

$$ (f_{\ell,P}^{\Delta})'(0) = P^{\ell} - a_{\ell}P^{\ell*}, \quad \text{and} \quad (f_{\ell,Q}^{\Delta})'(0) = Q^{\ell} - a_{\ell}Q^{\ell*}, $$

Let us remind that for all $t \in \mathbb{R}$ small enough, we have

$$ E_{\varepsilon, \delta}^{\Delta}(\psi^{1*} + t\psi^1, \psi^{2*} + t\psi^2) = E_{\varepsilon, \delta}^{\Delta}(\psi^{1*}, \psi^{2*}) + t\langle \nabla E_{\varepsilon, \delta}^{\Delta}(\psi^{1*}, \psi^{2*}), (\psi^1, \psi^2) \rangle_{\Delta} + O(t^2). $$

The next step is the computation of $(f^{\Delta})'(0)$ in order to derive a criterion for the minimizer. We obtain, using the definition (2.7) of $f^{\Delta}$ and the chain rule,

$$ \begin{aligned} (f^{\Delta})'(0) &= \langle \nabla E_{\varepsilon, \delta}^{\Delta} \left( \begin{pmatrix} f_{1,P}^{\Delta}(0) \\ f_{1,Q}^{\Delta}(0) \end{pmatrix}, \begin{pmatrix} f_{2,P}^{\Delta}(0) \\ f_{2,Q}^{\Delta}(0) \end{pmatrix} \right), \left( \begin{pmatrix} (f_{1,P}^{\Delta})'(0) \\ (f_{1,Q}^{\Delta})'(0) \end{pmatrix}, \begin{pmatrix} (f_{2,P}^{\Delta})'(0) \\ (f_{2,Q}^{\Delta})'(0) \end{pmatrix} \right) \rangle_{\Delta} \\ &= \sum_{\ell=1,2} \sum_{n,k=1}^{N} \left( \frac{\partial E_{\varepsilon, \delta}^{\Delta}}{\partial P^{\ell}} (\psi^{\ell*}) * (P^{\ell} - a_{\ell}P^{\ell*}) + \frac{\partial E_{\varepsilon, \delta}^{\Delta}}{\partial Q^{\ell}} (\psi^{\ell*}) * (Q^{\ell} - a_{\ell}Q^{\ell*}) \right)_{n,k}. \end{aligned} $$

As expected, $(f^{\Delta})'(0)$ is a linear function of $(P^1, Q^1, P^2, Q^2)$ of the form:

$$ (f^{\Delta})'(0) = \sum_{\ell=1}^{2} \sum_{n,k=1}^{N} \left( [K_{\ell,P}^{\Delta} * P^{\ell}]_{n,k} + [K_{\ell,Q}^{\Delta} * Q^{\ell}]_{n,k} \right), \quad (2.10) $$
with $K_{\ell,P}^{\Delta}$ and $K_{\ell,Q}^{\Delta}$ matrices in $\mathbb{R}^{N \times N}$ for $\ell \in \{1, 2\}$ given by

$$K_{\ell,P}^{\Delta} = \frac{\partial E_{\varepsilon,\delta}^{\Delta}}{\partial P^{\ell}}(\psi^{\ell*}) - \frac{\delta_x^2}{\|\psi^{\ell*}\|_{\Delta}^2} \sum_{n,k=0}^{N+1} \left( \frac{\partial E_{\varepsilon,\delta}^{\Delta}}{\partial P^{\ell}}(\psi^{\ell*}) * P^{\ell*} + \frac{\partial E_{\varepsilon,\delta}^{\Delta}}{\partial Q^{\ell}}(\psi^{\ell*}) * Q^{\ell*} \right)_{n,k} P^{\ell*},$$

and

$$K_{\ell,Q}^{\Delta} = \frac{\partial E_{\varepsilon,\delta}^{\Delta}}{\partial Q^{\ell}}(\psi^{\ell*}) - \frac{\delta_x^2}{\|\psi^{\ell*}\|_{\Delta}^2} \sum_{n,k=0}^{N+1} \left( \frac{\partial E_{\varepsilon,\delta}^{\Delta}}{\partial P^{\ell}}(\psi^{\ell*}) * P^{\ell*} + \frac{\partial E_{\varepsilon,\delta}^{\Delta}}{\partial Q^{\ell}}(\psi^{\ell*}) * Q^{\ell*} \right)_{n,k} Q^{\ell*}.$$

Using (2.9) and the expression given by (2.10), if $(P_1^*, Q_1^*, P_2^*, Q_2^*)$ minimizes $E_{\varepsilon,\delta}^{\Delta}$ under the constraints (2.6), then for any $(\psi^1, \psi^2)$, $(f^{\Delta})'(0) = 0$ is equivalent to $K_{\ell,P}^{\Delta} = K_{\ell,Q}^{\Delta} = 0$ for $\ell \in \{1, 2\}$. $\square$

Therefore, we use the smallness of $(K_{\ell,P})_{\ell=1,2}$ and $(K_{\ell,Q})_{\ell=1,2}$ in the method described in Section 4 as a criterion for numerical convergence. To this end, we set

$$K^{\Delta} = \sum_{\ell=1}^{2} \|K_{\ell,P}^{\Delta}\|_{\Delta} + \|K_{\ell,Q}^{\Delta}\|_{\Delta}.$$

**Remark 16** In the same way, we can use an adapted criterion for the minimization of $\mathcal{E}_{\varepsilon,\delta}^{\Delta}$.

# 4 Gradient Method

In order to minimize the discrete energy $E_{\varepsilon,\delta}^{\Delta}$ (in the form (2.1) or (2.4)) under the constraints (2.6), we propose the following gradient method with projection and with adaptive step.

## Initialization:

1. Choose $\psi^{1,1}, \psi^{2,1} \in \mathbb{C}^{(N+2)^2}$ normalized as: $\|\psi^{\ell,1}\|_{\Delta} = \sqrt{N_{\ell}M}$, ($\ell = 1, 2$) and satisfying the homogeneous Dirichlet boundary conditions,

2. choose the step $h > 0$,

3. choose two tolerance parameters $h_0, K_0 > 0$ for the convergence test.

## Iteration:

1. Compute the gradient: $\nabla E_{\varepsilon,\delta}^{\Delta}(\psi^{1,m}, \psi^{2,m})$,

2. Compute the auxiliary step $(\tilde{\psi}^{1,m}, \tilde{\psi}^{2,m}) = (\psi^{1,m}, \psi^{2,m}) - h \nabla E_{\varepsilon,\delta}^{\Delta}(\psi^{1,m}, \psi^{2,m})$, and set the homogeneous Dirichlet boundary conditions on $\tilde{\psi}^{1,m}$ and $\tilde{\psi}^{2,m}$,

3. Normalize the auxiliary step to obtain an attempt for the next step:

$$\psi^{\ell,m+1} = \frac{\tilde{\psi}^{\ell,m}}{\|\tilde{\psi}^{\ell,m}\|_{\Delta}} \sqrt{N_{\ell}M}.$$
4. If $E_{\varepsilon, \delta}^{\Delta}(\psi^{1,m+1}, \psi^{2,m+1}) > E_{\varepsilon, \delta}^{\Delta}(\psi^{1,m}, \psi^{2,m})$ then we replace $h$ by $\frac{h}{2}$ (provided $\frac{h}{2} \ge h_0$, otherwise we stop without convergence) and we compute a new auxiliary step $\tilde{\psi}^{1,m}$, $\tilde{\psi}^{2,m}$ for the same $m$ by going back to step 2. Else we compute $K^{\Delta}$ at point $\psi^{1,m+1}$, $\psi^{2,m+1}$ using (2.8). If $K^{\Delta} \ge K_0$ then we replace $m$ by $m+1$ and we start at step 1. Else we take $\psi^{1,m+1}$, $\psi^{2,m+1}$ as an approximate minimizer and we stop with convergence.

# 5 Post processing algorithms

## 5.1 Indices computation of single vortices

In this subsection, we introduce an algorithm that we developed for the computation of the indices of the vortices in the minimizers of the discrete energy functional. This algorithm is used in the Chapter 3.

The algorithm relies on 4 numerical parameters $tol_1 > 0$, $tol_2 > 0$, $N_{min} \in \mathbb{N}^*$ and $N_{max} \in \mathbb{N}^*$ with $N_{min} \le N_{max}$. It follows the 4 steps below. The three first steps identify vortices' centers and the last step computes the vortices' indices. In this section, $\psi$ denotes the squared complex matrix with $N^2$ entries for one of the two species of the Bose-Einstein condensate.

First, we determine the potential centers of the vortices and establish a lot of candidates as the set of $(n, k) \in \{1, \dots, N\}^2$ such that $|\psi_{n,k}|^2 < tol_1$.

Second, we build a second list $\mathbb{P}$ based on the first list above using the following rule. For each potential center $(n, k)$ in the list established in step 1, we consider the values of $|\psi|^2$ on the squares

$$ \mathbb{S}_{\lambda}(n, k) = \left\{ (n \pm \lambda, k + j) \mid j \in \{-\lambda, \dots, \lambda\} \right\} \cup \left\{ (n + j, k \pm \lambda) \mid j \in \{-\lambda, \dots, \lambda\} \right\}, $$

of length $2\lambda\delta_x$, for $\lambda \in \{N_{min}, \dots, N_{max}\}$. If for one of these $\lambda = \lambda(n, k)$, the values of $|\psi|^2$ at all points of the square $\mathbb{S}_{\lambda}(n, k)$ are such that $|\psi|^2 - |\psi_{n,k}|^2 > tol_2$, then we add the center $(n, k)$ to the second list $\mathbb{P}$, and we set $\lambda(n, k)$ as the characteristic length of the potential vortex. In other words, we have determined a list $\mathbb{P}$ of couple of points $(n, k)$ satisfying the following conditions:

* $|\psi_{n,k}|^2 < tol_1$,

* $|\psi_{i,j}|^2 > |\psi_{n,k}|^2 + tol_2$, for all the couples $(i, j)$ in $\mathbb{S}_{\lambda(n,k)}(n, k)$.

Third step, we consider each center $(n, k)$ from the list $\mathbb{P}$ and we identify if another center is inside the square $\cup_{1 \le \lambda \le \lambda(n,k)} \mathbb{S}_{\lambda}(n, k)$. If this is the case, we eliminate the center with the biggest $|\psi|^2$ at the center from the list. We repeat this step until we are left with isolated centers. Let us denote by $\mathbb{T}$ all the couples $(n, k)$ corresponding to isolated centers. The last step consists in computing the indices. For each $(n, k) \in \mathbb{T}$, we start off with any couple $(i, j)$ in $\mathbb{S}_{\lambda(n,k)}(n, k)$. Then we calculate their associated angle $\theta_0 = \arg(\psi_{i,j})$. Note that there are $8\lambda(n, k)$ couples in $\mathbb{S}_{\lambda(n,k)}(n, k)$. After computing the first angle $\theta_0$, we proceed to compute the other angles $\theta_1, \theta_2, \dots, \theta_{8\lambda(n,k)}$ in the following way:

* After computing the angle $\theta_m$, the next angle $\tilde{\theta}_{m+1}$ is computed as an argument of the next value of $\psi$ on the square $\mathbb{S}_{\lambda(n,k)}(n, k)$ with anticlockwise orientation (see Figure 2.1).
<figure><img src="image_7.png" /><figcaption>Figure 2.1 – $S_λ(n, k)$ with angles $θ_0$ to $θ_{8λ}$.</figcaption></figure>

• We set $θ_{m+1} := \tilde{θ}_{m+1} + 2k\pi$ with $k = \operatorname{argmin}_{l \in \mathbb{Z}} |\tilde{θ}_{m+1} - θ_m + 2\pi l|$.

Finally the index of the vortex $(n, k)$ is set to $(θ_{8λ(n,k)} - θ_0)/2π$.

## 5.2 Indices computation of vortex sheet

In this subsection, we introduce an algorithm that we developed for the computation of the index along a vortex sheet of a minimizer of the discrete energy functional. Let us note that this algorithm *does* need the interference of a human being at one point. This algorithm is used in Chapter 3.

The algorithm relies on 4 main steps in total and it uses 3 numerical values $m \le M$ (close to 1/2) and a small $tol_3 > 0$. The first step consists in identifying the contours of the vortex sheets. The second step optimizes these contours. The third step sorts the optimized contours, in order to compute, in the last step, each vortex sheet index alongside its contour. We will use the example in Figure 2.2 to illustrate the different steps of the algorithm.

The contour of a vortex sheet in one of the components of the condensate consists in a region where the squared modulus of the minimizer moves fast from 0 to 1 = $\max ρ$ (or the opposite). The first step consists in finding regions on the discrete square $[-L, L]^2$ where this occurs. We identify, for each species $ψ^ℓ$, the coordinates on the grid for the values $|ψ^ℓ|^2$ between $m$ and $M$. Next, we add the grid coordinates close to the circle $\partial D$, of which we retain only the coordinates where $|ψ^ℓ|^2 \le tol_3$. These grid points constitute the union of potential contours of vortex sheets (see Figure 2.3a for an example). Let us sort these points, to ease the search for connex components in contours in the next step. Let us note $\mathbb{K}$ the set of coordinates we have found so far. We are looking for a union of *connex* curves that describe the
<figure><img src="image_8.png" /><figcaption>(a) The first component of a minimizer for the energy $E_{\varepsilon, \delta}^{\Delta}$ presenting vortex sheets (see Figure 3.36a for details on the parameters).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b) An example of the contours detected after the first part of the first step of the contour detection algorithm.</figcaption></figure>

Figure 2.2 – An example of the contours detected after the first step of the contour detection algorithm ($m = 0.4, M = 0.6$, and $tol_3 = 0.3$).

<figure><img src="image_10.png" /><figcaption>(a) An example of the contours detected after the first part of the first step of the contour detection algorithm.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) An example of the *connex* curves obtained at the end of the first step of the contour detection algorithm.</figcaption></figure>

Figure 2.3 – A particular sheet from the example displayed in Figure 2.2 after the first step of the contour detection algorithm ($m = 0.4, M = 0.6$, and $tol_3 = 0.3$).
borders of each vortex sheet. Each curve can be determined with as many points as we want, but we want to avoid taking too many points per curve. Therefore, we want to limit the number of neighbours to any point on the curve to 3. For each grid point of coordinates $(i, j)$, we denote by $S_1(i, j)$ the set of its neighbours. For all grid point $(i, j) \in \mathbb{K}$, we act as follows:

1. If $(i, j)$ has no neighbours in $\mathbb{K}$, we just remove it from the set $\mathbb{K}$.

2. If $(i, j)$ has 1 neighbour, we add to the set $\mathbb{K}$ a couple $(i', j') \notin \mathbb{K}$ verifying

$$ (i', j') = \argmin_{(n,k) \in S_1(i,j)} ||\psi_{i,j}|^2 - |\psi_{n,k}|^2|. $$

3. If $(i, j)$ has 2 or 3 neighbours in $\mathbb{K}$, we do nothing.

4. If $(i, j)$ has at least 4 neighbours in $\mathbb{K}$, we remove the couple $(i', j')$ from $\mathbb{K}$ verifying:

$$ (i', j') = \argmax_{(n,k) \in S_1(i,j) \cap \mathbb{K}} ||\psi_{i,j}|^2 - |\psi_{n,k}|^2|, $$

and we repeat this step until we are left with only 3 neighbours.

At this point, we have detected a union of different curves defining the borders of the vortex sheets. An example if displayed in Figure 2.3b. Our next step is to separate the connex components in $\mathbb{K}$, to categorize the different sheets (since each minimizer can have more than one vortex sheet). Therefore, step 2 consists in associating to

<figure><img src="image_12.png" /><figcaption>(a) After Step 2: We can use a different colour of each contour (the colour changes with $p$).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b) After Step 3: We can use colour shading to show the anticlockwise orientation of each contour.</figcaption></figure>

Figure 2.4 – An example of the contours detected and orientated by the contour detection algorithm after the second and third step.

each border an integer number $p \ge 1$. At the $p^{th}$ vortex sheet, we create a new list $K_p$. We start by adding to $K_p$ a grid point at random from the list $\mathbb{K}$ whilst removing it from $\mathbb{K}$. Then, we add all of its neighbours one by one whilst removing them from $\mathbb{K}$ until we are left with no neighbours in $\mathbb{K}$ to all the grid points in $K_p$. Of course, if there are still points in $\mathbb{K}$, we increase $p$ to $p+1$ and start this step again until $\mathbb{K}$ is empty. After dividing the contours of the vortex sheets into categories, we look at all the categories and let a human decide if a detected category $K_p$ makes sense or not, as well as whether two categories should be merged into one. This is
the only human interaction needed in the sheet detection algorithm. To merge two categories $K_q$ and $K_{q'}$ we follow these simple steps. First, we search for the closest two couples $(n, k)$ and $(n', k')$ from $K_q$ and $K_{q'}$ (using the usual distance in $\mathbb{R}^2$) and we compute their middle $(i', j') = (\lfloor \frac{n+n'}{2} \rfloor, \lfloor \frac{k+k'}{2} \rfloor)$. If this middle has neighbours in both categories, we merge them into a new set $K_q \cup K_{q'} \cup (i', j')$. Otherwise, we repeat once the same process between $K_q \cup \{(i', j')\}$ and $K_{q'} \cup \{(i', j')\}$. In the end, we are left with different categories for different vortex sheet's contours. An example of what we obtain after Step 2 is displayed in Figure 2.4a.

Step 3 consists in building, from each set $K_p$ a list of grid points corresponding to an anticlockwise path along the contour of the $p^{th}$ vortex sheet. We proceed in the following way:

1. For each list $K_p$, we compute its barycentre. We choose our starting couple from $K_p$ as one with a close $x$-axis coordinate to that of the barycentre and with the biggest $y$-axis coordinate. We add this to our new list $K_p^{sort}$.

2. While the length of $K_p^{sort}$ is less than 70% of the original length of $K_p$ and the last point and the first point of $K_p^{sort}$ are not neighbours:

(a) Out of all the neighbours in $K_p$ to the last grid point we added to $K_p^{sort}$, we choose one and prioritize the anticlockwise direction. Then, we remove this grid point from $K_p$.

(b) If we don't find any neighbour, we delete this grid point from $K_p^{sort}$ and from $K_p$. Next, we return to the last grid point we added to $K_p$ before that and we repeat the previous step.

An example of what we obtain after Step 3 is displayed in Figure 2.4b.

The last and fourth step consists in computing the indices of each vortex sheet contour detected $K_p^{sort}$. The algorithm is similar to the last step in Section 5.1 but we take the list $K_p^{sort}$ instead of $S_{\lambda(n,k)}(n, k)$.

**Remark 17** In the case of a circular giant hole, we know it is going to be a disk shaped hole centered in $(0,0)$. We compute the radius of the disk using a simple algorithm that goes as follows.

Let $\psi$ be a minimizer with a giant hole and let $R$ be a set of radius $r_l > 0$ with $l \in \{0, \dots, l_0\}$ and $\theta_m = \frac{2\pi m}{m_0}$ an angle with $m \in \{0, \dots, m_0\}$ for $m_0 \gg 1$. At step $l$, we compute the following set

$$ \beta_l = \{(n, k) = \underset{(n,k) \in \{1, N\}^2}{\operatorname{argmin}} (|x_n - r_l \cos(\theta_t)|, |y_k - r_l \sin(\theta_t)|) \text{ for all } t \in \{0, \dots, m_0\}\}. $$

Note that $\beta_l$ contains at most $m_0$ unique couples. If $\max_{n,k \in \beta_l} |\psi_{n,k}|^2 > 0.1$, then we stop the algorithm and we will use $r_l$ and $\beta_l$ for the computation of the index. Otherwise we repeat this step with $l = l + 1$. Finally, to compute the index of the giant hole, we use the last step of Section 5.1 with $\beta_l$ instead of $S_{\lambda(n,k)}(n, k)$.

**Remark 18** Following Remark 17, in the case of a square giant hole (in the case of high rotation single-component minimizer using DF scheme), we know it is going to be a square shaped hole tilted by angle of $45^\circ$ centered in $(0,0)$. We compute the
side $c$ of the square using the radius $r$ found in Remark 17 with this simple relation
$c = \sqrt{2}r$. We then compute the set defining the tilted square $C$. Finally, to compute
the index of the square hole, we use the last step of Section 5.1 with $C$ instead of
$\mathbb{S}_{\lambda(n,k)}(n, k)$.

In the next chapter, we will present the numerical results obtained using the gra-
dient method with different numerical and physical parameters in order to illustrate
numerically the theory presented in Chapter 1 Sections 4.2, 4.3 and 4.4.
# Chapter 3

## Numerical results

This chapter is devoted to the numerical results in 2D of the numerical minimizers of the Gross-Pitaevskii energy using the Explicit with Projection Gradient method (EPG method) developed in Chapter 2. We display results covering the regimes discussed in Chapter 1 Sections 4.2, 4.3 and 4.4.

This chapter is divided into three main parts. The first part is dedicated to displaying the results of the finite difference scheme for one component and two components BEC. The second part of this chapter focuses on the numerical minimizers using the Fourier transformation approach. Both parts also display the post processing algorithms developed in Section 4 of Chapter 2 to validate the existence of numerical vortices and vortex sheets. The last part is devoted to the comparison between EPG and GPELab (see [9, 11]). We display different tests of one component and two components numerical minimizers for low and moderate number of points in the discretization to show the efficiency of the EPG method.

### 1 Common parameters used in all of the simulations

For all the simulations in this chapter, we consider the following parameters. The confinement is defined by the function

$$ \rho(x, y) = \min[1, 10(R^2 - x^2 - y^2)]. \qquad (3.1) $$

The discretization parameters are the following: the square's length $L = 7$, the radius $R = 4$, the initial step $h = 0.1$ and the tolerance $h_0 = 10^{-12}$.

### 2 Numerical results for a discretization of $E_{\varepsilon, \delta}^{\Delta}$ using Finite Difference scheme

This first section is devoted to displaying the numerical minimizers of Gross-Pitaevskii energy using finite difference scheme (2.1). We cover all the regimes described in Chapter 1.

We first cover the one component cases with and without rotation. We then move on to the two components condensates results. We display both segregation and coexistence regimes with and without rotation.
### 2.1 One component condensate without rotation

#### Parameters of the numerical simulations

In this first subsection, we consider the following parameters. For the initial datum we choose $\psi^0(x, y) = \exp(-10x^2 - 10y^2)/5$. We choose then for the rotational velocity $\Omega = 0$, the values $N_1 = 1$, $N_2 = 0$, (since $N_2 = 0$ then $\delta$ plays no role here). For the rest of the parameters, we refer to Section 1 of this chapter.

#### Numerical results

In this simulation of one component condensate, we study the case of no rotation ($\Omega = 0$) and strong confinement ($\epsilon = 10^{-2}$). We set the number of points in the $x$-axis and $y$-axis to $N = K = 512$. In this simulation, the gradient descent algorithm stopped due to the stopping criterion $K^{\Delta} \le 10^{-2}$. The results are shown in Figure 3.1b. As we can see, there are no vortices present in the numerical minimizer. This is in accordance with the theory presented in Section 4.2 (first case $\Omega < \Omega_{\epsilon}^{1}$).

<figure><img src="image_14.png" /><figcaption>Figure 3.1 – Comparison between the positive part of $\rho$ and the squared modulus of a minimizer for the energy $E^{\Delta}$ with no rotation $\Omega = 0$.</figcaption></figure>

In Figure 3.2, we present alongside the axis $y = 0$, the function $\rho(x, y)$ (in blue), the positive part of $\rho(x, y)$ (in black), the squared modulus of a minimizer for the energy $E_{\epsilon}^{\Delta}$ when $\Omega = 0$, for $N = K = 512$, for $\epsilon = 0.1$ (in green mark x) and for $\epsilon = 0.01$ (in red mark x). As we can see, the Figure 3.2 numerically illustrates the convergence of the squared modulus of a minimizer for the energy $E_{\epsilon}^{\Delta}$ to $\rho$ when $\epsilon$ tends to 0 at fixed rotational speed $\Omega = 0$. We observe that, as $\epsilon$ gets smaller, the squared modulus of a minimizer of the energy converges to the positive part of $\rho$ and vanished outside of $[-R, R]$. This is in accordance with Theorem 31 in Chapter 4.
<figure><img src="image_15.png" /><figcaption>Figure 3.2 – The graph of ρ(x, 0) (in blue), the positive part of ρ(x, 0) (in black), the squared modulus of a minimizer for the energy $E_{\varepsilon}^{\Delta}$ for $\varepsilon = 0.1$ (in green mark x) and for $\varepsilon = 0.01$ (in red mark x) alongside the axis $y = 0$.</figcaption></figure>

## 2.2 One component condensate with rotation

We now consider the theoretical results described in Chapter 1 Section 4.2. The results identify four different regimes for the behavior of the one component Bose-Einstein condensate depending on how big $\Omega_{\varepsilon}$ is as $\varepsilon \to 0$. These four regimes are separated by three characteristic rotational speeds $\Omega_{\varepsilon, i}$, $i = 1, 2, 3$ (see (1.14)). We explain below how we identify the four different regimes numerically. We conclude with numerical simulations for small $\varepsilon$ in each one of the four regimes.

### Identification of the regimes

In order to identify the four regimes described in Chapter 1 Section 4.2, we proceed as follows. First, we take $\varepsilon = 10^{-1}$ and we use the gradient algorithm described in Section 4 of Chapter 2 to compute minimizers of $E_{\varepsilon, \delta}^{\Delta}$ for several rotational speeds. Then, we do the same for $\varepsilon = 5 \times 10^{-2}$. Based on these many simulations, and the expressions (1.14), we estimate the three critical values for the rotation speed which are shown in Figure 3.3. This provides us with an estimation of the three critical values of the rotational speed when $\varepsilon = 10^{-2}$. Choosing rotation speed between these 3 estimated critical values for $\varepsilon = 10^{-2}$, we manage to observe the four different regimes, as detailed below.

In view of Figure 3.3, for $\varepsilon = 10^{-2}$ we choose the following values for the rotational velocity: $\Omega = \{1, 30, 60, 85, 110\}$. We choose these values to exhibit one case of the first two zones and three cases for the third zone. The case $\Omega = 60$ (Figure 3.6) is close enough to zone 2 and allows to see the transition from zone 2 to zone 3. The case $\Omega = 85$ (Figure 3.7) exhibits an annulus with a hole looking like a square. This is due to a grid orientation effect and might disappear if we would change numerical parameters such as the size of the computational domain. The last case $\Omega = 110$ (Figure 3.8) also displays an annulus with vortices but now the square shaped hole is more clear.
<figure><img src="image_16.png" /><figcaption>Figure 3.3 – $\Omega$ as a function of $\epsilon$ and the associated zones.</figcaption></figure>

### Parameters of the numerical simulations

For all the one component simulations with rotation ($\Omega \neq 0$), we choose $\psi^0 = \frac{1}{5}\exp(-10x^2 - 10y^2)$ as initial datum, the values $N_1 = 1$ and $N_2 = 0$. For the rest of the parameters we refer to Section 1 of this chapter.

### Small rotational speed ($\Omega < \Omega_{\epsilon}^{1}$)

In a first simulation of one component condensate, we study the case of a very low rotational velocity ($\Omega = 1$) and strong confinement ($\epsilon = 10^{-2}$). We set the number of points in the x-axis and y-axis to $N = K = 256$. In this simulation, the gradient descent algorithm converged due to the stopping criterion $K^{\Delta}$ which has a value smaller than $5 \times 10^{-4}$ after 213300 iterations. The results are shown in Figure 3.4.

<figure><img src="image_17.png" /><figcaption>Figure 3.4 – The squared modulus of a minimizer of the energy $E_{\epsilon, \delta}^{\Delta}$ for $\Omega = 1$ and $\epsilon = 10^{-2}$.</figcaption></figure>
As we can see, there are no vortices present in the numerical minimizer. This is
due to the low rotational velocity and is in accordance with the theory presented in
Chapter 1 Section 4.2 (first case $\Omega < \Omega_{\varepsilon}^{1}$).

Moderate rotational speed ($\Omega_{\varepsilon}^{1} < \Omega < \Omega_{\varepsilon}^{2}$)

In a second simulation of one component condensate in rotation, we study the case
of a moderate rotational velocity ($\Omega = 30$) and strong confinement ($\varepsilon = 10^{-2}$). We
set the number of points in the *x*-axis and *y*-axis to $N = K = 512$. Then, we
interpolate its real and imaginary parts to a new grid with $N = K = 1024$. In this
simulation, the gradient descent algorithm converged due to the stopping criterion
$K^{\Delta}$ which has a value smaller than $5 \times 10^{-3}$ after 714900 iterations. The results are
shown in Figure 3.5a. As we can see, there are vortices in the numerical minimizer.
This is due to the moderate rotational velocity and is in accordance with the theory
presented in Chapter 1 Section 4.2 (second case $\Omega_{\varepsilon}^{1} < \Omega < \Omega_{\varepsilon}^{2}$).

<figure><img src="image_18.png" /><figcaption>(a) The squared modulus of a minimizer
of a one component condensate.</figcaption></figure>

<caption>(b) The vortices indices of a minimizer
of a one component condensate.</caption>

Figure 3.5 – The squared modulus of a minimizer (a) of the energy $E_{\varepsilon,\delta}^{\Delta}$ for $\Omega = 30$ and $\varepsilon = 10^{-2}$. In the right panel, the numerical computation of the indices of the vortices is carried out with $N_{min} = 1$, $N_{max} = 5$, $tol_{1} = 0.05$ and for $tol_{2} = 0.01$.

In Figure 3.5b, we compute the indices of the numerical vortices with $N_{min} = 1$, $N_{max} = 5$, $tol_{1} = 0.05$ and $tol_{2} = 0.01$ (see Chapter 2 Section 5.1 for a description of the method). As we can see, the index of all the numerical vortices is equal to one. This validates numerically that the zeros of the function are singly quantized vortices.

**Remark 19** Using the forward finite difference scheme, we encounter a square
shape displayed in some of the minimizers of the energy $E_{\varepsilon,\delta}^{\Delta}$ due to the orienta-
tion grid effect (see for example Figures 3.6a and 3.7a). This only occurs when $\Omega$
is big enough for the centrifugal force to come into play. We avoid this problem by
choosing the Fast Fourier transform approach instead (see simulations in Section 3
of this chapter).
### Big rotational speed ($\Omega_{\epsilon}^2 < \Omega < \Omega_{\epsilon}^3$)

In a third simulation of one component condensate, we study the case of a fairly high rotational velocity ($\Omega = 60$) and strong confinement ($\epsilon = 10^{-2}$). We set the number of points in the x-axis and y-axis to $N = K = 512$. In this simulation, the gradient descent algorithm converged due to the step $h$ being too small ($h < tol$) while the stopping criterion $K^\Delta$ has a value of 0.015 after 1164200 iterations. The results are shown in Figure 3.6a. As we can see, there is a giant hole in the center, surrounded by vortices on the annulus of the numerical minimizer. This is due to the centrifugal force coming into play with the fairly high rotational velocity and is in accordance with the theory presented in Chapter 1 Section 4.2 (third case $\Omega_{\epsilon}^2 < \Omega < \Omega_{\epsilon}^3$).

<figure><img src="image_19.png" /><figcaption>(a) The squared modulus of a minimizer of the energy $E_{\epsilon,\delta}^{\Delta}$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(b) The vortices indices for $N_{min} = 1$, $N_{max} = 5$, $tol_1 = 0.05$ and $tol_2 = 0.01$.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(c) The vortices indices of a minimizer of the energy with $r = 0.782$ (see Remark 17))</figcaption></figure>

**Figure 3.6** – The squared modulus of a minimizer (a) of the energy $E_{\epsilon,\delta}^{\Delta}$ for $\Omega = 60$ and $\epsilon = 10^{-2}$. The vortices indices of the minimizer are presented in Figures (b) and (c).

In Figure 3.6b, we compute the index of the numerical vortices of Figure 3.6a with $N_{min} = 1$, $N_{max} = 5$, $tol_1 = 0.05$ and $tol_2 = 0.01$ (see Chapter 2 Section 5.1 for a description of the method). As we expected, the index of all the numerical vortices in Figure 3.6b around the giant hole is equal to one. In Figure 3.6c, we compute the index of the giant hole using the algorithm described in Remark 17 with $r = 0.782$. The index of the giant hole is equal to 30. This is in accordance with the theory presented in Chapter 1 Section 4.2 (third case).

In a fourth simulation of one component condensate, we study the case of a high rotational velocity ($\Omega = 85$) and strong confinement ($\epsilon = 10^{-2}$). We set the number of points in the x-axis and y-axis to $N = K = 512$. In this simulation, the gradient descent algorithm converged due to the stopping criterion $K^\Delta$ that has a value smaller than $9 \times 10^{-3}$ after 1150200 iterations. The results are shown in Figure 3.7a. As we can see, the giant hole in the center is now bigger, surrounded by less vortices on the "annulus" of the numerical minimizer than before. In this simulation we can see the giant hole is actually becoming square shaped (see Remark 19). This is due to the centrifugal force coming into play with a very high rotational velocity.

In Figure 3.7b, we compute the index of the numerical vortices of Figure 3.7a with $N_{min} = 1$, $N_{max} = 5$, $tol_1 = 0.05$ and $tol_2 = 0.01$ (see Chapter 2 Section 5.1 for a description of the method). As we can see, the index of all the numerical vortices
<figure><img src="image_22.png" /><figcaption>(a) The squared modulus of a minimizer of the energy $E_{\varepsilon, \delta}^{\Delta}$.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(b) The vortices indices of a minimizer for $N_{min} = 1$, $N_{max} = 5$, $tol_1 = 0.05$ and for $tol_2 = 0.1$.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(c) The index of the giant vortex with $c=2.014$ (see Remark 18).</figcaption></figure>

Figure 3.7 – The squared modulus of a minimizer (a) of the energy $E_{\varepsilon,\delta}^{\Delta}$ for $\Omega = 85$ and $\varepsilon = 10^{-2}$. The vortices indices of the minimizer are presented in Figures (b) and (c).

in Figure 3.7b around the giant hole is equal to one. In Figure 3.7c, we compute the index of the giant hole using the algorithm described in Remark 18 with $c = 2.014$. The index of the giant vortex is now equal to 105. This is in accordance with the theory presented in Chapter 1 Section 4.2 (third case).

### Huge rotational speed ($\Omega > \Omega_{\varepsilon}^3$)

In this last simulation of one component condensate, we study the case of a huge rotational velocity ($\Omega = 110$) and strong confinement ($\varepsilon = 10^{-2}$). We set the number of points in the x-axis and y-axis to $N = K = 512$. In this simulation, the gradient descent algorithm converged due to the stopping criterion $K^{-\Delta}$ that has a value smaller than $10^{-2}$ after 1443214 iterations. The results are shown in Figure 3.8a. As we can see, the giant hole is becoming a central giant vortex, surrounded by less vortices on the "annulus" of the numerical minimizer than before. In this simulation we can also see the giant hole is actually a square shaped hole (see Remark 19). This is due to the centrifugal force coming into play with a rotational velocity and is in accordance with the theory presented in Chapter 1 Section 4.2 (fourth case $\Omega \approx \Omega_{\varepsilon}^3$).

In Figure 3.8b, we compute the index of the numerical vortices of Figure 3.8a with $N_{min} = 1$, $N_{max} = 5$, $tol_1 = 0.1$ and $tol_2 = 0.02$ (see Chapter 2 Section 5.1 for a description of the method). As we can see, the index of all the numerical vortices in Figure 3.7b around the giant hole is equal to one. In Figure 3.8c, we compute the index of the giant hole using the algorithm described in Remark 18 with $c = 2.54$. The index of the giant vortex is now equal to 161. This is in accordance with the theory presented in Chapter 1 Section 4.2 (last case).
<figure><img src="image_25.png" /><figcaption>(a) The squared modulus of a minimizer of the energy $E_{\varepsilon,\delta}^{\Delta}$.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(b) The vortices indices of a minimizer for $N_{min} = 1$, $N_{max} = 5$, $tol_1 = 0.1$ and for $tol_2 = 0.02$.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(c) The index of the giant vortex with c=2.54 (see Remark 18).</figcaption></figure>

Figure 3.8 – The squared modulus of a minimizer (a) of the energy $E_{\varepsilon,\delta}^{\Delta}$ for $\Omega = 110$ and $\varepsilon = 10^{-2}$. The vortices indices of the minimizer are presented in Figures (b) and (c).

### 2.3 Two components condensate without rotation ($\Omega = 0$)

We move now on the theoretical results described in the first part of Chapter 1 Sections 4.3 and 4.4.

#### Parameters used for the simulations

In this subsection, we consider the following parameters. For the initial datum we choose

$$ \psi^1(x, y) = \psi^2(x, y) = \exp(-10x^2 - 10y^2)/5. \qquad (3.2) $$

We use a rotational velocity $\Omega = 0$, a confinement $\varepsilon = 5 \times 10^{-2}$, and the values $N_1 = 0.55$, $N_2 = 0.45$. The stopping criterion value for $K^{\Delta}$ is set to $K_0 = 10^{-2}$ and we set the number of point in the x-axis and y-axis to $N = 256$. For the rest of the parameters, we refer to Section 1 of this chapter.

#### Strong confinement in a strong segregation regime ($\delta_{\varepsilon} \to \infty$ as $\varepsilon \to 0$)

<figure><img src="image_28.png" /><figcaption>(a) The squared modulus of the first component.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>(b) The squared modulus of the second component.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>(c) The sum of the squared modulus of both components.</figcaption></figure>

Figure 3.9 – The squared modulus of the first component (a) and the second component (b) of a minimizer of the energy $E_{\varepsilon,\delta}^{\Delta}$ in the case of two components condensate with no rotation and $\delta = 4000$. Figure (c) displays their sum.
In a first simulation of two components condensate without rotation, we study the case of strong segregation regime ($\delta = 4000$) and strong confinement ($\varepsilon = 5 \times 10^{-2}$). The results are shown in Figure 3.9. As we can see in Figures 3.9a and 3.9b, we are in the segregation regime since the two components tend not to overlap. In Figure 3.9c, we can see the sum of the squared modulus of the two components displays a red curve which separates the squared modulus of the two components with $\min_{(x,y)\in D} |\psi^{1*}|^2 + |\psi^{2*}|^2 \approx 0.243$. This is in accordance with the theory presented in Chapter 1 Section 4.4 (first case of no rotation with $\delta_{\varepsilon}\varepsilon^2 \to +\infty$).

Strong confinement in a moderate segregation regime ($\delta_{\varepsilon}$ fixed as $\varepsilon \to 0$)

<figure><img src="image_31.png" /><figcaption>(a) The squared modulus of the first component.</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>(b) The squared modulus of the second component.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>(c) The sum of the squared modulus of both components.</figcaption></figure>

Figure 3.10 – The squared modulus of (a) the first component and (b) the second component of a minimizer for a two components condensate and (c) their sum, in the case of no rotation ($\Omega = 0$) and $\delta = 1.5$.

In a second simulation of two components condensate without rotation, we study the case of moderate segregation regime ($\delta = 1.5$) and strong confinement ($\varepsilon = 5 \times 10^{-2}$). The results are shown in Figure 3.10. As we can see in Figures 3.10a and 3.10b, we are in the segregation regime ($\delta > 1$). In Figure 3.10c, we can see that the sum of the squared modulus of the two components displays a white curve which separates the squared modulus of the two components with $\min_{(x,y)\in D} |\psi^{1*}|^2 + |\psi^{2*}|^2 \approx 0.78$. This is in accordance with the theory presented in Chapter 1 Section 4.4 (second case with no rotation where $\delta > 1$ is fixed).

Strong confinement in a weak segregation regime ($\delta_{\varepsilon} \to 1$ as $\varepsilon \to 0$)

In a last simulation of two components condensate without rotation, we study the case of segregation regime ($\delta = 1.02$) and strong confinement ($\varepsilon = 5 \times 10^{-2}$). The results are shown in Figure 3.11. As we can see in Figures 3.11a and 3.11b, we are in the segregation regime ($\delta > 1$). In Figure 3.11c, we can see that the sum of the squared modulus of the two components does not display a separation area with $\min_{(x,y)\in D} |\psi^{1*}|^2 + |\psi^{2*}|^2 \approx 0.971$. With the notations of Chapter 1 Section 4.4 we have here $\tilde{\varepsilon} = 0.35$. This is in accordance with the theory presented in Chapter 1 Section 4.4 (second case with no rotation where $\delta_{\varepsilon} \to 1$ and $\tilde{\varepsilon} \to 0$).
<figure><img src="image_34.png" /><figcaption>Figure 3.11 – The squared modulus of the first component (a) and the second component (b) of a minimizer of the energy $E_{\epsilon,\delta}^{\Delta}$ in the case of two components condensate and their sum (c) in the case of no rotation and $\delta = 1.02$.</figcaption></figure>

### 2.4 Two components condensate with rotation ($\Omega \neq 0$) in a segregation regime ($\delta_{\epsilon} > 1$)

In this subsection, we consider the theoretical results described in Chapter 1 Section 4.4 for the two components rotating Bose-Einstein condensate in the segregation regime.

**Parameters used for a two components condensate with rotation**

For all the two components segregation regime with rotation, we consider the following parameters. For the initial datum we choose $\psi^1, \psi^2$ as defined in (3.2). We choose then the numerical values $N_1 = 0.55, N_2 = 0.45$. The stopping criterion value for $K^{\Delta}$ is set to $K_0 = 10^{-2}$ and we set the number of points in the x-axis and y-axis to $N = 256$. Then, we interpolate the real and imaginary parts of the minimizer to a grid with $N = 512$. We also set the physical parameters $\epsilon = 10^{-2}$ and $\delta = 1 + \epsilon = 1.01$ so that $\tilde{\epsilon} = \sqrt{\epsilon} = 0.1$. For the rest of the parameters, we refer to Section 1 of this chapter.

**Low rotation case : $\Omega = 1$**

We first consider the case of low rotation $\Omega = 1$ and strong confinement ($\epsilon = 10^{-2}$). The results are shown in Figure 3.12. The numerical experiment confirms that we are in a segregation regime (recall that $\delta > 1$) and the two components tend to not overlap. Moreover, the small rotational velocity $\Omega = 1$ is not big enough to produce vortices in the minimizer. This is in accordance with the theory presented in Chapter 1 Section 4.4 (first bullet point in the rotational case).
<figure><img src="image_35.png" /><figcaption>(a) The squared modulus of the first component.</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>(b) The squared modulus of the second component.</figcaption></figure>

Figure 3.12 – The squared modulus of the first component (a) and the second component (b) for a minimizer in the case of two components condensate with low rotation $\Omega = 1$ in the segregation regime.

### Moderate rotational case ($\Omega = 3$)

Next we consider the case of a moderate rotation $\Omega = 3$ and strong confinement ($\varepsilon = 10^{-2}$). The results are shown in Figure 3.13. The numerical experiment shows that the segregation holds and the moderate rotational velocity $\Omega = 3$ is now big enough to produce singly quantized vortices in one of the components of the numerical minimizer. This is in accordance with the theory presented in Chapter 1 Section 4.4 (second bullet point in the rotational case).

<figure><img src="image_37.png" /><figcaption>(a) The squared modulus of the first component.</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>(b) The squared modulus of the second component.</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>(c) The vortices indices of the first component for $N_{min} = 1, N_{max} = 5, tol_1 = 0.05$ and $tol_2 = 1 = 0.02$.</figcaption></figure>

Figure 3.13 – The squared modulus of the first component (a) and the second component (b) for a minimizer in the case of two components condensate with low rotation $\Omega = 3$ in the segregation regime. The vortices indices of the first component are displayed in Figure (c).

### High rotational speed ($\Omega = 6$)

Next we consider the case of high rotation $\Omega = 6$ and strong confinement ($\varepsilon = 10^{-2}$). The results are shown in Figure 3.14. The numerical experiment shows that the
segregation holds and the high rotational velocity $\Omega = 6$ is now big enough to produce singly quantized vortices in both components of the numerical minimizer. This is in accordance with the theory presented in Chapter 1 Section 4.4 (second bullet point in the rotational case).

<figure><img src="image_40.png" /><figcaption>(a) The squared modulus of the first component.</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>(b) The squared modulus of the second component.</figcaption></figure>

<figure><img src="image_42.png" /><figcaption>(c) The vortices indices of the first component.</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>(d) The vortices indices of the second component.</figcaption></figure>

Figure 3.14 – The squared modulus of the first component (a) and the second component (b) for a minimizer in the case of a two components condensate with high rotation $\Omega = 6$ in the segregation regime. The vortices indices of the first and second component are displayed in Figures (c) and (d) respectively with $N_{min} = 1$, $N_{max} = 5$, $tol_1 = 0.05$ and $tol_2 = 0.01$.

### Huge rotational speed ($\Omega = 15$)

In the last and final numerical simulation of the segregation case, we consider the case of huge rotation $\Omega = 15$ and strong confinement ($\epsilon = 10^{-2}$). The results are shown in Figure 3.15. The numerical experiment shows that the segregation holds since the support of the minimizers tend to *not* overlap. Moreover, the huge rotational velocity $\Omega = 15$ is big enough to produce vortex sheets in both components of the numerical minimizer. As we can see in Figures 3.15b and 3.15e, almost all the indices of the numerical vortices are equal to one. Some vortices appear with an index of 2. This is due to the vortex detection algorithm which is not able to
identify vortices close to each others. Moreover, in Figures 3.15c and 3.15f, we use the post processing algorithm we developed in Section 5.2 of Chapter 2 to detect and compute the index of vortex sheets. The existence of vortex sheets is in accordance with the conjecture presented in Chapter 1 Section 4.4 (last bullet point in the rotational case). Indeed, this numerical simulation supports the theoretical possibility of having vortex sheets.

<figure><img src="image_44.png" /><figcaption>(a) The squared modulus of the first component.</figcaption></figure>

<figure><img src="image_45.png" /><figcaption>(b) The vortices indices of the first component.</figcaption></figure>

<figure><img src="image_46.png" /><figcaption>(c) The vortex sheet indices of the first component.</figcaption></figure>

<figure><img src="image_47.png" /><figcaption>(d) The squared modulus of the second component.</figcaption></figure>

<figure><img src="image_48.png" /><figcaption>(e) The vortices indices of the second component.</figcaption></figure>

<figure><img src="image_49.png" /><figcaption>(f) The vortex sheet indices of the second component.</figcaption></figure>

Figure 3.15 – The squared modulus of the first component (a) and the second component (d) for a minimizer in the case of two components condensate with huge rotation $\Omega = 15$ in the segregation regime. The vortices indices of the first and second component are displayed in Figures (b) and (e) respectively with $N_{min} = 1$, $N_{max} = 2$, $tol_1 = 0.05$ and $tol_2 = 10^{-5}$. The vortex sheets of the first and second component are displayed in Figures (c) and (f) with $m = 0.4$, $M = 0.6$ and $tol_3 = 0.3$.

**Remark 20** In order to detect most of the vortices, it is necessary that $N_{max}$ be close to $N_{min} = 1$ and $tol_2$ small enough.

## 2.5 Two components condensate in a coexistence regime ($\delta_{\varepsilon} < 1$)

In this subsection, we study the numerical behaviour of two components Bose-Einstein condensate in the coexistence regime ($\delta < 1$) as $\varepsilon \to 0$. Also, in this subsection, due to the high rotational speed in some simulations, we look for the minimizers of the energy $E_{\varepsilon,\delta}^{\Delta}$ with the centrifugal force (see Remark 2 in Chapter 1) in order to find the different regimes.
As mentioned in Chapter 1 Section 4.3, depending on the value of $\delta$ and $\Omega$, we identify four different regimes.

### Parameters used for the coexistence regimes

For all the two components coexistence regime with rotation, we consider the following parameters. For the initial datum we choose, unless stated otherwise,

$$ \begin{aligned} \psi^1(x, y) &= \exp(-10(x - 0.5)^2 - 10(y + 0.3)^2)/5, \\ \psi^2(x, y) &= \exp(-10(x + 0.7)^2 - 10(y - 0.1)^2)/5. \end{aligned} $$

We choose then the numerical values $N_1 = 0.55$, $N_2 = 0.45$. The stopping criterion value for $K^{\Delta}$ is set to $K_0 = 10^{-2}$ and we set the number of points in the x-axis and y-axis to $N = 128$. Then, we interpolate the real and imaginary parts of the minimizer to a grid with $N = 256$. We also set the physical parameter $\epsilon = 5 \times 10^{-2}$. For the rest of the parameters, we refer to Section 1 of this chapter.

### No vortices

In a first simulation of two components condensate in coexistence regime ($\delta = 0.3$), we study the case of no rotation ($\Omega = 0$) and strong confinement ($\epsilon = 5 \times 10^{-2}$). We use the same initial datum as (3.2). The results are shown in Figure 3.16.

<figure><img src="image_50.png" /><figcaption>Figure 3.16 – The squared modulus of the first component (a), and the second component (b) of a minimizer for the energy $E_{\epsilon,\delta}^{\Delta}$ with $\Omega = 0$</figcaption></figure>

We can see in Figure 3.16 that we are in the coexistence regime since each component is disk-shaped and both components are overlapping. Also, notice that the maximum value of the squared modulus of each component depends on the values of $N_1 = 0.55$ and $N_2 = 0.45$. This is in accordance with the theory presented in Chapter 1 Section 4.3, for the first case.

### Triangular vortex lattices

In a second simulation of a two components condensate in coexistence regime, we study the case of moderate rotational speed ($\Omega = 5$), strong confinement ($\epsilon =$
$5 \times 10^{-2}$) and weak interaction strength ($\delta = 0.1$). The results are shown in Figure 3.17. We can see in Figure 3.17 that we are in the coexistence regime where each

<figure><img src="image_51.png" /><figcaption>(a) The squared modulus of the first component.</figcaption></figure>

<figure><img src="image_52.png" /><figcaption>(b) The squared modulus of the second component.</figcaption></figure>

<figure><img src="image_53.png" /><figcaption>(c) The vortex indices of the first component.</figcaption></figure>

<figure><img src="image_54.png" /><figcaption>(d) The vortex indices of the second component.</figcaption></figure>

Figure 3.17 – The squared modulus of the first component (a), and the second component (b) of a minimizer for the energy $E_{\varepsilon,\delta}^{\Delta}$ for $\Omega = 5$, and $\delta = 0.1$. The vortices indices of the first component (c), and the second component (d), followed by red lines highlighting the triangular lattice.

component is disk-shaped. Also in Figure 3.17, we compute the index of the vortices detected with $N_{min} = 1$, $N_{max} = 5$, $tol_1 = 5 \times 10^{-2}$, $tol_2 = 10^{-2}$. As expected, all the indices of the numerical vortices are equal to one which validates numerically that the zeros of the function have a singly quantized phase circulation. Finally, the red lines highlighting the vortex lattice in Figures 3.17c and 3.17d form mostly triangles which is in accordance with the theory presented in Chapter 1 Section 4.3.

**Remark 21** In Figures 3.17, 3.18, 3.38 and 3.39, we consider all the pairs of vortices detected by our algorithm, and draw a red line between the two vortices of the pair if their euclidean distance is below some threshold, which is taken close to the mean distance between 2 vortices.
### Square vortex lattices

In a third simulation of a two components condensate in coexistence regime ($\delta_{\epsilon} < 1$), we study the case of high rotational speed ($\Omega = 6$), strong confinement ($\delta_{\epsilon} = 5 \times 10^{-2}$) and strong interaction strength ($\delta = 0.8$). The results are shown in Figure 3.18. We can see in Figure 3.18 that we are in the coexistence regime where

<figure><img src="image_55.png" /><figcaption>(a) The squared modulus of the first component.</figcaption></figure>

<figure><img src="image_56.png" /><figcaption>(b) The squared modulus of the second component.</figcaption></figure>

<figure><img src="image_57.png" /><figcaption>(c) The vortex indices of the first component.</figcaption></figure>

<figure><img src="image_58.png" /><figcaption>(d) The vortex indices of the second component.</figcaption></figure>

Figure 3.18 – The squared modulus of the first component (a), and the second component (b) of a minimizer for the energy $E_{\epsilon, \delta}^{\Delta}$ with $\Omega = 6$ and $\delta = 0.8$. The vortices indices of the first component (c), and the second component (d), followed by red lines highlighting the square lattice.

each component is disk-shaped. Also in Figure 3.18, we compute the indices of the vortices detected with $N_{min} = 1$, $N_{max} = 3$, $tol_1 = 5 \times 10^{-2}$, $tol_2 = 2 \times 10^{-2}$. As expected, all the indices of the numerical vortices are equal to one which validates numerically that the zeros of the function are singly quantized vortices. Finally, the red lines highlighting the vortex lattice in Figures 3.18c and 3.18d form mostly squares which is in accordance with the theory presented in Chapter 1 Section 4.3.

### Double core and Stripe vortex lattice

In this last regime, for $\delta$ close to 1, we observe depending on the initial datum and the rotational speed, either a stripe vortex lattice or a double core vortex lattice.
For a first simulation, we consider the case of fairly high rotational speed ($\Omega = 4$) and strong confinement ($\epsilon = 5 \times 10^{-2}$) and strong interaction strength ($\delta = 0.99$). The results are shown in Figure 3.19.

<figure><img src="image_59.png" /><figcaption>(a) The squared modulus of the first component.</figcaption></figure>

<figure><img src="image_60.png" /><figcaption>(b) The squared modulus of the second component.</figcaption></figure>

<figure><img src="image_61.png" /><figcaption>(c) The vortex indices of the first component.</figcaption></figure>

<figure><img src="image_62.png" /><figcaption>(d) The vortex indices of the second component.</figcaption></figure>

Figure 3.19 – The squared modulus of the first component (a), and the second component (b) of a minimizer for the energy $E_{\epsilon,\delta}^{\Delta}$ with $\Omega = 4$, $\epsilon = 5 \times 10^{-2}$ and $\delta = 0.99$. The vortices indices of the first component (c) and the second component (d) for $N_{min} = 1$, $N_{max} = 5$, $tol_1 = 0.05$ and $tol_2 = 0.02$.

We can see in Figure 3.19 that we are in the coexistence regime where each component is disk-shaped (although the second component is small on the border of the disk, but comparing it to the segregation regime the case where $\delta \to 1$, we can see the disk shape). In Figures 3.19c and 3.19d, we compute the indices of the vortices detected. First, as expected, almost all the vortices are paired up 2 by 2. Second, all the indices of the numerical vortices are equal to one which validates numerically that the zeros of the function have a singly quantized phase circulation. This is in accordance with the theory presented in Chapter 1 Section 4.3, last case.

For the last simulation, we use the following initial datum:

$$ \psi_{n,k}^{1} = \psi_{n,k}^{2} = \sin(x_{n,k} + y_{n,k}), $$
still with the convention that $\psi^1$ and $\psi^2$ satisfies the homogeneous Dirichlet boundary conditions. We consider the case of fairly moderate rotational speed ($\Omega = 3$) and strong confinement ($\epsilon = 5 \times 10^{-2}$) and very strong interaction strength ($\delta = 0.99$). The results are shown in Figure 3.20.

<figure><img src="image_63.png" /><figcaption>(a) The squared modulus of the first component.</figcaption></figure>

<figure><img src="image_64.png" /><figcaption>(b) The squared modulus of the second component.</figcaption></figure>

<figure><img src="image_65.png" /><figcaption>(c) The vortex indices of the first component.</figcaption></figure>

<figure><img src="image_66.png" /><figcaption>(d) The vortex indices of the first component.</figcaption></figure>

Figure 3.20 – The squared modulus of the first component (a), and the second component (b) of a minimizer for the energy $E_{\epsilon, \delta}^{\Delta}$ with $\Omega = 3$ and $\delta = 0.99$. The vortices indices of the first component (c) and the second component (d) for $N_{min} = 1$, $N_{max} = 5$, $tol_1 = 0.05$ and $tol_2 = 0.01$.

We can see in Figure 3.20 that we are in the coexistence regime still where each component is disk-shaped. Also in Figure 3.20, we compute the index of the vortices detected. As expected, all the indices of the numerical vortices are equal to one which validates numerically that the zeros of the function have a singly quantized phase circulation. We can see a stripe pattern in both components (Figures 3.20a, 3.20b), this is in accordance with the theory presented in Chapter 1 Section 4.3, last phase.
# 3 Numerical results for a discretization of the energy using Fast Fourier Transform

This second section is devoted to displaying the numerical minimizers of GP energy using Fourier transform (2.4). We cover all the regimes described in Chapter 1.

Similarly to Section 2, we first cover the one component cases with and without rotation. We then move on to the two components condensates results. We display both segregation and coexistence regimes with and without rotation.

## 3.1 One component condensate without rotation

In this first subsection, we present numerical results for one component condensate with no rotation. They show the convergence of $|\psi|^2 \to \rho$ as $\varepsilon$ tends to 0. This result will be proved later in Chapter 4 Theorem 31.

We introduce below the parameters used for the one component simulations.

### Parameters used for the one component simulations

For all the one component simulations we consider the following parameters. For the initial datum we choose $\psi^0(x, y) = \exp(-10x^2 - 10y^2)/5$. We choose the confinement parameter $\varepsilon = 10^{-2}$, the values $N_1 = 1$, $N_2 = 0$, (since $N_2 = 0$ then $\delta$ plays no role here). For the rest of the parameters, we refer to Section 1 of this chapter.

### Numerical simulations

In this first simulation of a one component condensate without rotation ($\Omega = 0$) we consider a strong confinement ($\varepsilon = 10^{-2}$). We set the number of points in the x-axis and y-axis to $N = K = 512$. In this simulation, the gradient descent algorithm converged due to the stopping criterion $K^\Delta \le 10^{-2}$. The results are shown in Figure 3.21b. As we can see, there are no vortices in the numerical minimizer. This is in accordance with the theory presented in Chapter 1 Section 4.2 (first case $\Omega < \Omega_\varepsilon^1$). Moreover the squared modulus of a minimizer converges to the positive part of the function $\rho$ as we can see in Figure 3.21a.

In Figure 3.22, we present alongside the axis $y = 0$, the function $\rho(x, y)$ (in blue), the positive part of $\rho(x, y)$ (in black), the squared modulus of a minimizer for the energy $E_\varepsilon^\Delta$ for $\Omega = 0$, $N = K = 512$, $\varepsilon = 0.1$ (in green mark x) and $\varepsilon = 0.01$ (in red mark x). As we can see, the Figure 3.22 numerically illustrates the convergence of the squared modulus of a minimizer for the energy $E_\varepsilon^\Delta$ to $\rho$ when $\varepsilon$ tends to 0 at fixed rotational speed $\Omega = 0$. We observe that, as $\varepsilon$ gets smaller, the squared modulus of a minimizer of the energy converges to the positive part of $\rho$. This is in accordance with Theorem 31 in Chapter 4.
<figure><img src="image_67.png" /><figcaption>Figure 3.21 – Comparison between the positive part of $\rho$ and the squared modulus of a minimizer for the energy $E_{\varepsilon}^{\Delta}$ without rotation ($\Omega = 0$).</figcaption></figure>

<figure><img src="image_68.png" /><figcaption>Figure 3.22 – The graph of $\rho(x, 0)$ (in blue), the positive part of $\rho(x, 0)$ (in black), the squared modulus of a minimizer for the energy $E_{\varepsilon}^{\Delta}$ for $\varepsilon = 0.1$ (in green mark x) and for $\varepsilon = 0.01$ (in red mark x) alongside the axis $y = 0$.</figcaption></figure>

## 3.2 One component condensate with rotation

We consider the theoretical results described in Chapter 1 Section 4.2. The results identify four different regimes for the behaviour of the one component Bose-Einstein condensate depending on how big $\Omega_{\varepsilon}$ is as $\varepsilon$ tends to zero. These four regimes are separated by three characteristic rotational speeds $\Omega_{\varepsilon}^i$, $i = 1, 2, 3$ (see (1.14)). We refer to Section 3.1 of this chapter for the parameters we use for all the one component simulations. Then we explain how we identify the four different regimes numerically. We conclude with numerical simulations for small $\varepsilon$ in each one of the four regimes.
Identification of the four regimes

In order to identify the four regimes, we proceed as follows and as we did in Section 2. First we take $\epsilon = 10^{-1}$ and we use the gradient algorithm described in Section 4 of Chapter 1 to compute minimizers of $E_{\epsilon, \delta}^{\Delta}$ for several rotational speeds. Then we do the same for $\epsilon = 5 \times 10^{-2}$. Based on all these simulations, and the expressions (1.14), we identify the 3 critical values for the rotational speed which are shown in the Figure 3.23. This provides us with some critical values of the rotational speed when $\epsilon = 10^{-2}$. The four regimes numerically observed for $\epsilon = 10^{-2}$ that we present below are separated by these three critical values. In view of Figure 3.23, for $\epsilon = 10^{-2}$

<figure><img src="image_69.png" /><figcaption>Figure 3.23 – $\Omega$ as a function of $\epsilon$ and the associated zones using the discrete energy (2.4) with $N=K=128$.</figcaption></figure>

we choose the following values for the rotational velocity: $\Omega = \{1, 20, 30, 40, 50, 70\}$. We choose these values to exhibit at least one case of each zone.

We can see a difference between the estimation of the critical constants obtained with the finite difference discretization of the energy (Figure 3.3) and the Fourier transform discretization of the energy (Figure 3.23). This might be due to several things:

* The grid orientation effect caused by the discretization choice (see Remark 19). Each approach may require different discretization parameters to mimic the Dirichlet boundary conditions at the boundary of the disk.

* The Fourier transform discretization gives more importance to the rotational energy than the finite difference approach.

* The values $\epsilon = 10^{-1}$ and $\epsilon = 5 \times 10^{-2}$ are maybe not small enough to predict the critical constants for $\epsilon = 10^{-2}$ using both approaches.
#### Small rotational speed ($\Omega_{\varepsilon} < \Omega_{\varepsilon}^{1}$)

In a first simulation of a one component condensate, we consider the case of small rotational speed ($\Omega = 1$) and strong confinement ($\varepsilon = 10^{-2}$). This simulation is carried out with $N = 512$. In this simulation, the gradient descent algorithm converged due to the stopping criterion $K^{\Delta} \le 10^{-2}$. The results are shown in Figure 3.24. As we can see, there are no vortices in the numerical minimizer which indicates we are in the first zone. This is in accordance with the theory presented in Chapter 1 Section 4.2 (first case $\Omega < \Omega_{\varepsilon}^{1}$).

<figure><img src="image_70.png" /><figcaption>Figure 3.24 – The squared modulus of a minimizer for the energy $E_{\varepsilon}^{\Delta}$ with small rotational speed ($\Omega = 1$) and strong confinement ($\varepsilon = 10^{-2}$) computed with $N = 512$ and $K^{\Delta} = 10^{-2}$.</figcaption></figure>

#### Moderate rotational speed ($\Omega_{\varepsilon}^{1} < \Omega_{\varepsilon} < \Omega_{\varepsilon}^{2}$)

In a second simulation of a one component condensate, we consider the case of a moderate rotational velocity ($\Omega = 20$) and strong confinement ($\varepsilon = 10^{-2}$), according to Figure 3.3. First, we compute a minimizer with $N = 512$. Then we interpolate its real and imaginary parts and compute another minimizer with $N = 1024$. In this simulation, the gradient descent algorithm converged due to the stopping criterion $K^{\Delta}$ which has a value of $5 \times 10^{-2}$. We compute numerically the indices of the vortices using the algorithm described in Chapter 2 Section 5.1. The results are shown in Figure 3.25. As we can see, there are several vortices in the numerical minimizer and there is no sign of a giant hole in the center. This corresponds to the case of a moderate rotational speed and is in accordance with the second case of the theory presented in Chapter 1 Section 4.2.

In a third simulation of a one component condensate, we still consider the case of a moderate rotational velocity ($\Omega = 30$) and strong confinement ($\varepsilon = 10^{-2}$). Once again we use an interpolation from $N = 512$ to $N = 1024$. In this simulation, the gradient descent algorithm converged due to the stopping criterion $K^{\Delta}$ which has a value of $5 \times 10^{-2}$. We also compute the indices of the vortices numerically using the algorithm described in Chapter 2 Section 5.1. The results are shown in Figure 3.26. As we can see, there are vortices in the numerical minimizer and there is a dark disk in the center (not present in Figure 3.25) which indicates that we are
<figure><img src="image_71.png" /><figcaption>Figure 3.25 – The squared modulus of a minimizer for $\Omega = 20$ and $\varepsilon = 10^{-2}$. In the right panel, the numerical computation of the indices of the vortices is carried out with $N_{min} = 1$, $N_{max} = 3$, $tol_1 = 0.1$ and $tol_2 = 0.05$.</figcaption></figure>

close to the limit between zone 2 and zone 3 as shown in Figure 3.3. This is due to
the moderate rotational velocity and is in accordance with the theory presented in
Chapter 2 Section 4.2. Moreover, the numerical index of almost all the numerical
vortices is equal to one. This validates numerically that almost all the zeros of the
wave function have a singly quantized phase circulation.

<figure><img src="image_72.png" /><figcaption>Figure 3.26 – The squared modulus of a minimizer for $\Omega = 30$ and $\varepsilon = 10^{-2}$. In the right panel, the numerical computation of the indices of the vortices is carried out with $N_{min} = 1$, $N_{max} = 3$, $tol_1 = 0.05$ and $tol_2 = 0.02$.</figcaption></figure>

Big rotational speed ($\Omega_{\varepsilon}^2 < \Omega_{\varepsilon} < \Omega_{\varepsilon}^3$)

In a fourth simulation of a one component condensate, we consider the case of a fairly
high rotational velocity ($\Omega = 40$) and strong confinement ($\varepsilon = 10^{-2}$) in accordance
with the regimes obtained in Figure 3.23. Once again we use an interpolation from
$N = 512$ to $N = 1024$. The gradient descent algorithm converged due to the
stopping criterion $K^{\Delta} \le 5 \times 10^{-2}$. The results are shown in Figure 3.27a. As we
can see, there is a giant hole in the center, surrounded by vortices on an annulus in
the numerical minimizer. This is due to the centrifugal force coming into play with
the fairly high rotational velocity and is in accordance with the theory presented in
Chapter 2 Section 4.2.
<figure><img src="image_73.png" /><figcaption>(a) The squared modulus of a minimizer of the energy $E_{\varepsilon, \delta}^{\Delta}$.</figcaption></figure>

<figure><img src="image_74.png" /><figcaption>(b) The vortices' indices of a minimizer of the energy for $N_{min} = 1$, $N_{max} = 3$, $tol_1 = 0.05$ and $tol_2 = 0.02$.</figcaption></figure>

<figure><img src="image_75.png" /><figcaption>(c) The index of the giant hole of a minimizer of the energy for with $r = 1.523$ (see Remark 17).</figcaption></figure>

Figure 3.27 – The squared modulus of a minimizer for $\Omega = 40$ and $\varepsilon = 10^{-2}$. The vortices indices of the singly quantized vortices (a) and the giant hole (b).

In Figure 3.27, we compute the indices of the vortices in the numerical minimizer displayed in Figure 3.27a. First, in Figure 3.27b, we compute the indices of the vortices in the annulus. Then, in Figure 3.27c, we compute the index of the giant hole with $r = 1.523$ (see Remark 17). As we can see, all the indices of the numerical vortices in Figure 3.27b around the giant hole are equal to one which validates numerically that almost all the zeros of the wave function have a singly quantized phase circulation. In Figure 3.27c, the index of the giant hole is equal to 100. This is in accordance with the theory presented in Chapter 2 Section 4.2.

In a fifth simulation of a one component condensate, we study the case of a fairly high rotational velocity ($\Omega = 50$) and strong confinement ($\varepsilon = 10^{-2}$) in accordance with the regimes obtained in Figure 3.23. Once again we use an interpolation from $N = 512$ to $N = 1024$. In this simulation, the gradient descent algorithm converged due to the stopping criterion $K^{\Delta} \le 5 \times 10^{-2}$. The results are shown in Figure 3.28a. As we can see, the giant hole in the center is now bigger, surrounded by less vortices on the annulus of the numerical minimizer than before. This is due to the centrifugal force coming into play with a very high rotational velocity and is in accordance with the theory presented in Chapter 2 Section 4.2 (fourth case $\Omega_{\varepsilon}^3 < \Omega$).

In Figure 3.28, we compute the indices of the vortices of the numerical minimizer displayed in Figure 3.28a. Firstly, in Figure 3.28b, we compute the indices of the vortices on the annulus. Then, in Figure 3.28c, we compute the index of the giant hole. As we can see, all the indices of the numerical vortices in Figure 3.28b around the giant hole are equal to one which validates numerically that the zeros of the function have a singly quantized phase circulation. While in Figure 3.28c, the index of the giant hole is equal to 263. This is in accordance with the theory presented in Chapter 2 Section 4.2 (fourth case).
<figure><img src="image_76.png" /><figcaption>(a) The squared modulus of a minimizer of the energy $E_{\varepsilon, \delta}^{\Delta}$.</figcaption></figure>

<figure><img src="image_77.png" /><figcaption>(b) The vortices indices of a minimizer of the energy for $N_{min} = 1$, $N_{max} = 3$, $tol_1 = 0.05$ and for $tol_2 = 0.02$.</figcaption></figure>

<figure><img src="image_78.png" /><figcaption>(c) The vortices indices of a minimizer of the energy with $r = 2.188$ (see Remark 17))</figcaption></figure>

Figure 3.28 – The squared modulus of a minimizer (a) for $\Omega = 50$ and $\varepsilon = 10^{-2}$. The vortices indices of the singly quantized vortices (b) and the giant hole (c).

### Huge rotational speed ($\Omega_{\varepsilon}^3 < \Omega_{\varepsilon}$)

In a last simulation of one component condensate, we study the case of a huge rotational velocity ($\Omega = 70$) and strong confinement ($\varepsilon = 10^{-2}$). Once again, we use an interpolation from $N = 512$ to $N = 1024$. The gradient descent algorithm converged due to the stopping criterion $K^{\Delta} \le 7 \times 10^{-2}$. The results are shown in Figure 3.29a. As we can see, the giant hole in the center is now bigger, surrounded by less vortices on the annulus of the numerical minimizer than before. This is due to the centrifugal force coming into play with a huge rotational velocity and is in accordance with the theory presented in Chapter 2 Section 4.2 (fourth case $\Omega_{\varepsilon}^3 < \Omega$).

<figure><img src="image_79.png" /><figcaption>(a) The squared modulus of a minimizer of the energy $E_{\varepsilon, \delta}^{\Delta}$.</figcaption></figure>

<figure><img src="image_80.png" /><figcaption>(b) The vortices indices of a minimizer of the energy for $N_{min} = 1$, $N_{max} = 3$, $tol_1 = 0.05$ and for $tol_2 = 0.02$.</figcaption></figure>

<figure><img src="image_81.png" /><figcaption>(c) The vortices indices of a minimizer of the energy with $r = 2.85$ (see Remark 17))</figcaption></figure>

Figure 3.29 – The squared modulus of a minimizer (a) for $\Omega = 70$, $\varepsilon = 10^{-2}$. The vortices indices of the singly quantized vortices (b) and the giant hole (c).

In Figure 3.29, we also compute the indices of the vortices in the numerical minimizer displayed in Figure 3.29a. First, in Figure 3.29b, we compute the indices of the vortices in the annulus. Then, in Figure 3.29c, we compute the index of the giant hole. As we can see, all the indices of the numerical vortices in Figure 3.29b around the giant hole are equal to one which validates numerically that the zeros of the function have a singly quantized phase circulation. While in Figure 3.29c,
the index of the giant hole is equal to 560. This is in accordance with the theory presented in Chapter 2 Section 4.2 (fourth case).

## 3.3 Two components condensate without rotation ($\Omega = 0$)

We move on the theoretical results described in the first part of Chapter 1 Section 4.4.

### Parameters used for two components simulations without rotation

We consider the following parameters for all the simulations of this subsection. We choose the numerical values $N_1 = 0.55$, $N_2 = 0.45$. The stopping criterion value for $K^\Delta$ is set to $K_0 = 10^{-2}$ and we set the number of point in the x-axis and y-axis to $N = 256$. We then interpolate the real and imaginary parts of the minimizer to a grid with $N = 512$. We also set the physical parameters $\varepsilon = 5 \times 10^{-2}$ and the rotational velocity $\Omega = 0$. For the rest of the parameters, we refer to Section 1 of this chapter.

**Remark 22** For the second simulation, displayed in Figure 3.31, we choose the same initial datum as before (see (3.2)). As to the first and third simulations (Figures 3.30 and 3.32), we started with the numerical minimizer obtained in Figure 3.31 as initial datum without changing the other parameters.

### Strong confinement in a strong segregation regime ($\delta_\varepsilon \varepsilon^2 \to \infty$ as $\varepsilon \to 0$)

In a first simulation of two components condensate without rotation, we study the case of strong segregation regime ($\delta = 4000$) and strong confinement ($\varepsilon = 5 \times 10^{-2}$). This simulation is carried out with $N = 256$. In this simulation, the gradient descent algorithm converged due to the stopping criterion $K^\Delta \le 10^{-2}$. The results are shown in Figure 3.30. As we can see in Figures 3.30a and 3.30b, we are in the segregation regime ($\delta > 1$). In Figure 3.30c, we can see that the sum of the squared modulus of the two components present a red and white curve which separates the two components with $\min_{(x,y)\in D} |\psi^{1*}|^2 + |\psi^{2*}|^2 \approx 0.14$. This is in accordance with the theory presented in Chapter 1 Section 4.4 (first case with no rotation where $\delta_\varepsilon \varepsilon^2 \to +\infty$).
<figure><img src="image_82.png" /><figcaption>(a) The squared modulus of the first component.</figcaption></figure>

<caption>(b) The squared modulus of the second component.</caption>

<caption>(c) The sum of the squared modulus of both components.</caption>

Figure 3.30 – The squared modulus of the first component (a) and the second component (b) of a minimizer for a two components condensate for $\Omega = 0$, $\varepsilon = 5 \times 10^{-2}$ and $\delta = 4000$. Their sum is displayed in figure (c).

### Strong confinement in segregation regime ($\delta_{\varepsilon} > 1$ fixed as $\varepsilon \to 0$)

In a second simulation of two component condensate without rotation, we study the case of moderate segregation regime ($\delta = 1.5$) and strong confinement ($\varepsilon = 5 \times 10^{-2}$). This simulation is carried out with $N = 256$. In this simulation, the gradient descent algorithm converged due to the stopping criterion $K^{\Delta} \le 10^{-2}$. The results are shown in Figure 3.31. As we can see in Figures 3.31a and 3.31b, we are in the segregation regime ($\delta > 1$). In Figure 3.31c, we can see that the sum of the squared modulus of the two components present a white and blue curve which is the separation areas between the two components with $\min_{(x,y)\in D} |\psi^{1*}|^2 + |\psi^{2*}|^2 \approx 0.78$. This is in accordance with the theory presented in Chapter 1 Section 4.4 (second case with no rotation where $\delta > 1$ is fixed).

<figure><img src="image_83.png" /><figcaption>(a) The squared modulus of the first component.</figcaption></figure>

<caption>(b) The squared modulus of the second component.</caption>

<caption>(c) The sum of the squared modulus of both components of the minimizer.</caption>

Figure 3.31 – The squared modulus of the first component (a) and the second component (b) of a minimizer for two components condensate for $\Omega = 0$, $\varepsilon = 5 \times 10^{-2}$ and $\delta = 1.5$. Their sum is displayed in figure (c).

### Strong confinement in weak segregation regime ($\delta_{\varepsilon} \to 1$ as $\varepsilon \to 0$)

In a last simulation of two components condensate without rotation, we study the case of a weak segregation regime ($\delta = 1.02$) and strong confinement ($\varepsilon = 5 \times 10^{-2}$). This simulation is carried out with $N = 256$. In this simulation, the gradient descent algorithm converged due to the stopping criterion $K^{\Delta} \le 10^{-2}$. The results are shown
in Figure 3.32 and Figure 3.32c. As we can see in Figures 3.32a and 3.32b, we are in the segregation regime ($\delta > 1$). In Figure 3.32c, we can see that the sum of the squared modulus of the two components does not present a separation area of any kind and $\min_{(x,y)\in D} |\psi^{1*}|^2 + |\psi^{2*}|^2 \approx 0.976$. Here $\tilde{\varepsilon} = 0.35$. This is in accordance with the theory presented in Chapter 1 Section 4.4 (second case with no rotation where $\delta_{\varepsilon} \to 1$ and $\tilde{\varepsilon} \to 0$).

<figure><img src="image_84.png" /><figcaption>(a) The squared modulus of the first component.</figcaption></figure>

<figure><img src="image_85.png" /><figcaption>(b) The squared modulus of the second component.</figcaption></figure>

<figure><img src="image_86.png" /><figcaption>(c) The sum of the squared modulus of both components.</figcaption></figure>

Figure 3.32 – The squared modulus of the first component (a)) and the second component (b) of a minimizer for two components condensate for $\Omega = 0$, $\varepsilon = 5 \times 10^{-2}$ and $\delta = 1.02$. Their sum is displayed in figure (c).

## 3.4 Two components condensate with rotation ($\Omega \neq 0$) in a segregation regime ($\delta_{\varepsilon} > 1$)

We move now on the theoretical results described in the second part of Chapter 1 Section 4.4.

### Parameters used for the rotating two components simulations

For the next simulations, in the segregation regime $\delta_{\varepsilon} \to 1$ and $\varepsilon \ll \tilde{\varepsilon}$, we took the following parameters. We choose the same initial datum as before (3.2). We choose then the numerical values $N_1 = 0.55$, $N_2 = 0.45$. We also set the physical parameters $\varepsilon = 10^{-2}$, $\delta = 1 + \varepsilon = 1.01$ so that $\tilde{\varepsilon} = \sqrt{\varepsilon} = 0.1$. For the rest of the parameters, we refer to Section 1 of this Chapter.

### Small rotational speed

In a first simulation of a rotating two components condensate in the segregation regime, we study the case of low rotation ($\Omega = 1$) and strong confinement ($\varepsilon = 10^{-2}$). We use an interpolation from $N = 256$ to $N = 512$. In this simulation, the gradient descent algorithm converged due to the stopping criterion $K^{\Delta} \le 2 \times 10^{-2}$. The results are shown in Figure 3.33. As we can see, we are in the segregation regime with no numerical vortices presented in neither of the components of the numerical minimizer. This is in accordance with the theory presented in Chapter 1 Section 4.4 (the first bullet point in the rotational case).
<figure><img src="image_87.png" /><figcaption>(a) The squared modulus of the first component.</figcaption></figure>

<figure><img src="image_88.png" /><figcaption>(b) The squared modulus of the second component.</figcaption></figure>

Figure 3.33 – The squared modulus of the first component (a), and the second component (b) of a minimizer for a two components condensate for $\Omega = 1$.

## Moderate rotational speed

In a second simulation of a rotating two components condensate in the segregation regime, we study the case of moderate rotational velocity ($\Omega = 3$) and strong confinement ($\epsilon = 10^{-2}$). Once again, we use an interpolation from $N = 256$ to $N = 512$. In this simulation, the gradient descent algorithm converged due to the stopping criterion $K^\Delta \le 2 \times 10^{-2}$. The results are shown in Figure 3.34. As we can see, there are vortices in one of the components of the numerical minimizer 3.34b. This is in accordance with the theory presented in Chapter 1 Section 4.4 (second bullet point of the rotational case).

<figure><img src="image_89.png" /><figcaption>(a) The squared modulus of the first component.</figcaption></figure>

<figure><img src="image_90.png" /><figcaption>(b) The squared modulus of the second component.</figcaption></figure>

<figure><img src="image_91.png" /><figcaption>(c) The vortices indices of the second component for $N_{min} = 1$, $N_{max} = 5$, $tol_1 = 0.05$ and $tol_2 = 0.02$.</figcaption></figure>

Figure 3.34 – The squared modulus of the first component (a), and the second component (b) of a minimizer for a two components condensate with $\Omega = 3$. The vortices indices are displayed in Figure (c).

In Figure 3.34c, we compute the indices of the vortices of the second component of the minimizer. As we can see, all the indices of the numerical vortices are equal to one which validates numerically that the zeros of the function have a singly quantized phase circulation.
#### Big rotational speed

In a third simulation of a rotating two components condensate in the segregation regime, we study the case of moderate rotational velocity ($\Omega = 6$) and strong confinement ($\epsilon = 10^{-2}$). Once again we use an interpolation from $N = 256$ to $N = 512$. In this simulation, the gradient descent algorithm converged due to the stopping criterion $K^\Delta \le 2 \times 10^{-2}$. The results are shown in Figure 3.35. As we can see, there are vortices in both components of the numerical minimizer (Figures 3.35c and 3.35d). This is in accordance with the theory presented in Chapter 1 Section 4.4 (second bullet point of the rotational case).

<figure><img src="image_92.png" /><figcaption>(a) The squared modulus of the first component.</figcaption></figure>

<figure><img src="image_93.png" /><figcaption>(b) The squared modulus of the second component.</figcaption></figure>

<figure><img src="image_94.png" /><figcaption>(c) The vortices indices of the first component.</figcaption></figure>

<figure><img src="image_95.png" /><figcaption>(d) The vortices indices of the second component.</figcaption></figure>

Figure 3.35 – The squared modulus of the first component (a), and the second component (b) of a minimizer for a two components condensate with $\Omega = 6$. The vortices indices of the first component (c) and the second component (d) with $N_{min} = 1$, $N_{max} = 5$, $tol_1 = 0.05$ and $tol_2 = 0.02$.

Also in Figure 3.35, we compute the vortices indices of the numerical minimizer. First, in Figure 3.35c, we compute the indices of the vortices of the first component with $N_{min} = 1$, $N_{max} = 5$, $tol_1 = 0.05$ and $tol_2 = 0.02$. Second, in Figure 3.35b,
we compute the index of the second component with the same parameters. As we can see, all the indices of the numerical vortices are equal to one which validates numerically that the zeros of the function have a singly quantized phase circulation.

## Huge rotational speed

In a last simulation of a rotating two components condensate in the segregation regime, we study the case of high rotational velocity ($\Omega = 15$) and strong confinement ($\epsilon = 10^{-2}$). Once again we use an interpolation from $N = 256$ to $N = 512$. In this simulation, the gradient descent algorithm converged due to the stopping criterion $K^\Delta \le 2 \times 10^{-2}$. The results are shown in Figure 3.36. As we can see, we are in the segregation regime and there are vortex sheet patterns in each component of Figure 3.36. This is in accordance with the theory presented in Chapter 1 Section 4.4 (last bullet point of the rotational case).

<figure><img src="image_96.png" /><figcaption>(a) The squared modulus of the first component.</figcaption></figure>

<figure><img src="image_97.png" /><figcaption>(b) The vortex indices of the first component.</figcaption></figure>

<figure><img src="image_98.png" /><figcaption>(c) The vortex sheet indices of the first component.</figcaption></figure>

<figure><img src="image_99.png" /><figcaption>(d) The squared modulus of the second component.</figcaption></figure>

<figure><img src="image_100.png" /><figcaption>(e) The vortex indices of the second component.</figcaption></figure>

<figure><img src="image_101.png" /><figcaption>(f) The vortex sheet indices of the second component.</figcaption></figure>

Figure 3.36 – The squared modulus of the first component (a), and the second component (b) of a minimizer for a two components condensate with $\Omega = 15$. The vortex indices of the first component (c) and for the second component (d) with $N_{min} = 1$, $N_{max} = 2$, $tol_1 = 0.1$ and $tol_2 = 10^{-5}$. The vortex sheet indices of the first component (e) and for the second component (f) for $m = 0.4$, $M = 0.6$ and for $tol_3 = 0.3$.

Also in Figure 3.36, we compute the vortices indices of the numerical minimizer. In Figure 3.36b (Figure 3.36e), we compute the index of the vortices detected by our algorithm presented in Section 5.1 of Chapter 2 of the first (respectively second)
component with $N_{min} = 1$, $N_{max} = 2$, for $tol_1 = 0.1$ and for $tol_2 = 10^{-5}$. As we can see, almost all the indices of the numerical minimizer are equal to one which validates numerically that these numerical vortices have a singly quantized phase circulation. Since these vortices are condensed, in the meaning that they are so close to each others, it validates numerically the existence of vortex sheets when $\delta_{\epsilon} \to 1$ and $\tilde{\epsilon} \to 0$.

In order to detect all of the vortices and their indices, we developed a vortex sheet detection algorithm. Given a certain minimizer with vortex sheets, it identifies the contours of each vortex sheet, optimize these contours to create a *connex* path, reorganize the path into an anticlockwise path then compute the index of each vortex sheet detected (see Chapter 2 Section 5.2 for more details on the algorithm).

In Figures 3.36c and 3.36f, we compute the index of numerical vortex sheets. In Figure 3.36c (Figure 3.36f), we compute the index of the vortices sheets detected by our algorithm of the first (respectively second) component with $m = 0.4$, $M = 0.6$ and $tol_3 = 0.3$ (see Chapter 2 Section 5.2 for a description of the method). As we can see, the indices of the different numerical vortex sheets detected are strictly positive which validates numerically that these vortex sheets have a singly quantized phase circulation. This validates numerically the existence of vortex sheets when $\delta_{\epsilon} \to 1$ and $\tilde{\epsilon} \to 0$.

**Remark 23** Notice that the index of any vortex sheet in Figures 3.36c and 3.36f is greater or equal than the sum of the single vortices detected in the same vortex sheet in Figures 3.36b and 3.36e. For example let us compare Figures 3.36f and 3.36e. Consider the green vortex sheet in the center. The first algorithm detected 18 singly quantized vortices and the second algorithm detected one vortex sheet with index 18. In contrast, for the red vortex sheet, the first algorithm detected 19 singly quantized vortices, and the second algorithm detected one vortex sheet with index 21. This indicates that the single vortex detection algorithm fails to detect all the single vortices with these parameters.

## 3.5 Two components condensate in the coexistence regime

$$ \delta_{\epsilon} < 1 $$

In this subsection, we observe the numerical behaviour of a two components Bose-Einstein condensate in the coexistence regime ($\delta < 1$) as $\epsilon \to 0$. Also, in this subsection, due to the high rotational speed in some simulations, we look for the minimizers of the energy $E_{\epsilon,\delta}^{\Delta}$ with the centrifugal force (see Remark 2 in Chapter 1) in order to find the different regimes.

As mentioned in Chapter 1 Section 4.4, depending on the value of $\delta$ and $\Omega$, we identify four different regimes.

### Parameters used for the coexistence regimes

For all the two component coexistence regime, we consider the following parameters. Unless stated otherwise, we choose for the initial datum the same as (3.2). We also choose the numerical values $N_1 = 0.55$, $N_2 = 0.45$. The stopping criterion value for $K^{\Delta}$ is set to $K_0 = 10^{-2}$. We also set the physical parameter $\epsilon = 5 \times 10^{-2}$. For the rest of the parameters, we refer to Section 1 of this chapter.
## No vortices

In a first simulation of a two components condensate in coexistence regime ($\delta_{\varepsilon} < 1$), we study the case of no rotational velocity ($\Omega = 0$) and strong confinement ($\varepsilon = 5 \times 10^{-2}$). This simulation is carried out with $N = 128$ which we then interpolated to $N = 256$. In this simulation, the gradient descent algorithm converged due to the stopping criterion $K^{\Delta} \le 10^{-2}$. The results are shown in Figure 3.37.

<figure><img src="image_102.png" /><figcaption>(a) The squared modulus of the first component.</figcaption></figure>

<figure><img src="image_103.png" /><figcaption>(b) The squared modulus of the second component.</figcaption></figure>

Figure 3.37 – The squared modulus of the first component (a), and the second component (b) of a minimizer for the energy $E_{\varepsilon,\delta}^{\Delta}$ with $\Omega = 0$ and $\delta = 0.9$.

We can see in Figure 3.37 that we are in the coexistence regime since each component is disk-shaped depending on the values of $N_1 = 0.55$ and $N_2 = 0.45$. This is in accordance with the theory presented in Chapter 1 Section 4.3, the first case.

According to Theorem 41 of Chapter 4, starting with an initial datum that is symmetrical leads to a symmetrical minimizer. For the next 3 simulations (triangular/square/stripe vortex lattice), we will use the following initial data:

$$ \psi^1 = \frac{1}{5} \exp(-10(x - 0.5)^2 - 10(y + 0.3)^2), $$

$$ \psi^2 = \frac{1}{5} \exp(-10(x + 0.7)^2 - 10(y - 0.1)^2). $$

## Triangular vortex lattices

In a second simulation of a two components condensate in coexistence regime, we study the case of high rotational speed ($\Omega = 7$), strong confinement ($\varepsilon = 5 \times 10^{-2}$) and weak interaction strength ($\delta = 0.2$). This simulation is carried out with $N = 128$ which we then interpolated to $N = 256$. In this simulation, the gradient descent algorithm converged due to the stopping criterion $K^{\Delta} \le 10^{-2}$. The results are shown in Figure 3.38. We can see in Figure 3.38 that we are in the coexistence regime where each component is disk-shaped. Also in Figure 3.38, we compute the index of the vortices detected with $N_{min} = 1$, $N_{max} = 3$, $tol_1 = 5 \times 10^{-2}$, $tol_2 = 2 \times 10^{-2}$. As expected, all the indices of the numerical vortices are equal to one which validates numerically that the zeros of the function have a singly quantized phase circulation.
<figure><img src="image_104.png" /><figcaption>(a) The squared modulus of the first component.</figcaption></figure>

<figure><img src="image_105.png" /><figcaption>(b) The squared modulus of the second component.</figcaption></figure>

<figure><img src="image_106.png" /><figcaption>(c) The vortex indices of the first component.</figcaption></figure>

<figure><img src="image_107.png" /><figcaption>(d) The vortex indices of the second component.</figcaption></figure>

Figure 3.38 – The squared modulus of the first component (a), and the second component (b) of a minimizer for the energy with centrifugal force $\mathcal{E}_{\varepsilon,\delta}^{\Delta}$ for $\Omega = 7$, $\varepsilon = 5 \times 10^{-2}$ and $\delta = 0.2$. The vortices indices of the first component (c), and the second component (d), followed by red lines highlighting the triangular lattice.

Finally, the red lines highlighting the vortex lattice (see Remark 21) in Figures 3.38c and 3.38d form mostly triangles which is in accordance with the theory presented in Chapter 1 Section 4.3 second case.

#### Square vortex lattices

In a third simulation of a two components condensate in coexistence regime ($\delta_{\varepsilon} < 1$), we study the case of high rotational speed ($\Omega = 8$), strong confinement ($\delta_{\varepsilon} = 5 \times 10^{-2}$) and strong interaction strength ($\delta = 0.8$). This simulation is carried out with $N = 128$ which we then interpolated to $N = 256$. In this simulation, the gradient descent algorithm converged due to the stopping criterion $K^{\Delta} \le 10^{-2}$. The results are shown in Figure 3.39. We can see in Figure 3.39 that we are in the coexistence regime where each component is disk-shaped. Also in Figure 3.39, we compute the index of the vortices detected with $N_{min} = 1$, $N_{max} = 3$, $tol_1 = 5 \times 10^{-2}$, $tol_2 = 2 \times 10^{-2}$. As expected, all the indices of the numerical vortices are equal to one which validates numerically that the zeros of the function have a singly quantized phase circulation. Finally, the red lines highlighting the vortex lattice (see Remark
<figure><img src="image_108.png" /><figcaption>(a) The squared modulus of the first component.</figcaption></figure>

<figure><img src="image_109.png" /><figcaption>(b) The squared modulus of the second component.</figcaption></figure>

<figure><img src="image_110.png" /><figcaption>(c) The vortex indices of the first component.</figcaption></figure>

<figure><img src="image_111.png" /><figcaption>(d) The vortex indices of the second component.</figcaption></figure>

Figure 3.39 – The squared modulus of the first component (a), and the second component (b) of a minimizer for the energy with centrifugal force $\mathcal{E}_{\varepsilon,\delta}^{\Delta}$ for $\Omega = 8$, $\varepsilon = 5 \times 10^{-2}$ and $\delta = 0.8$. The vortices indices of the first component (c), and the second component (d), followed by red lines highlighting the square lattice.

21) in Figures 3.39c and 3.39d form mostly squares which is in accordance with the theory presented in Chapter 1 Section 4.3 third case.

Double core and Stripe vortex lattice

In this last regime, for $\delta$ close to 1, we observe depending on the initial datum and the rotational speed, either stripe vortex lattice or double core vortex lattice.

For a first simulation, we use the same initial datum $\psi^1$ and $\psi^2$ as before. We consider the case of fairly high rotational speed ($\Omega = 4$) and strong confinement ($\varepsilon = 5 \times 10^{-2}$). This simulation is carried out with $N = 128$ which we then interpolated to $N = 256$. In this simulation, the gradient descent algorithm converged due to the stopping criterion $K^{\Delta} \le 10^{-2}$. The results are shown in Figures 3.40a and 3.40d.

Due to the high rotation, the centrifugal force comes into play. Indeed, a dark area appears in the middle of each component of the minimizer. Therefore, in Figures 3.40b and 3.40e, in addition to minimizing the energy $E_{\varepsilon,\delta}^{\Omega}$, we compute a minimizer of the energy $\mathcal{E}_{\varepsilon,\delta}^{\Delta}$ with centrifugal force for the same parameters, taking
the minimizer found in Figures 3.40a and 3.40d as initial data. We can see in Figure

<figure><img src="image_112.png" /><figcaption>(a) The squared modulus of the first component.</figcaption></figure>

<figure><img src="image_113.png" /><figcaption>(b) The squared modulus of the first component.</figcaption></figure>

<figure><img src="image_114.png" /><figcaption>(c) The vortex indices of the first component $N_{min} = 3$, $N_{max} = 4$.</figcaption></figure>

<figure><img src="image_115.png" /><figcaption>(d) The squared modulus of the second component.</figcaption></figure>

<figure><img src="image_116.png" /><figcaption>(e) The squared modulus of the second component.</figcaption></figure>

<figure><img src="image_117.png" /><figcaption>(f) The vortex indices of the second component for $N_{min} = 1$, $N_{max} = 4$.</figcaption></figure>

Figure 3.40 – The squared modulus of the first component (a) (resp. (b)), and the second component (d) (resp. (e)) of a minimizer for the energy $E_{\epsilon, \delta}^{\Delta}$ (resp. the energy with centrifugal force $\mathcal{E}_{\epsilon, \delta}^{\Delta}$) with $\Omega = 4$ and $\delta = 0.995$. The vortices indices of the first component (e) and the second component (f) for $tol_1 = 0.05$ and $tol_2 = 0.01$.

3.40 that we are in the coexistence regime where each component is disk-shaped (although the second component is small on the border of the disk, but comparing it to the segregation regime the case where $\delta \to 1$, we can see the disk shape). In Figures 3.40c and 3.40f, we compute the indices of the vortices detected. First, as expected, almost all the vortices are paired up 2 by 2. Second, almost all the indices of the numerical vortices are equal to one which validates numerically that the zeros of the function have a singly quantized phase circulation except for one vortex in the first component 3.40c which has an index of 2. This also presents a double core vortex that are too close to each other which cannot be detected as two different vortices by our index detecting algorithm. This is in accordance with the theory presented in Chapter 1 Section 4.3, last case.

For the second and last simulation, we use the following initial data:

$$ \psi_{n,k}^{1} = \psi_{n,k}^{2} = \sin(x_{n,k} + y_{n,k}), $$

still with the convention that $\psi^1$ and $\psi^2$ satisfies the homogeneous Dirichlet boundary conditions. We consider the case of fairly high rotational speed ($\Omega = 5$) and strong confinement ($\epsilon = 5 \times 10^{-2}$) and very strong interaction strength ($\delta = 0.98$). This simulation is carried out with $N = 128$, then we interpolate to $N = 256$. In this
simulation, the gradient descent algorithm converged due to the stopping criterion $K^\Delta \le 10^{-2}$. The results are shown in Figures 3.41a and 3.41d.

Due to the high rotation, the centrifugal force comes into play again and a dark area appears in the middle of each component of the minimizer. Therefore, in Figures 3.41b and 3.41e, we computed a minimizer of the energy $\mathcal{E}_{\varepsilon, \delta}^\Delta$ with centrifugal force for the same parameters, taking the minimizer found in Figures 3.41a and 3.41d as initial data.

<figure><img src="image_118.png" /><figcaption>(a) The squared modulus of the first component.</figcaption></figure>

<figure><img src="image_119.png" /><figcaption>(b) The squared modulus of the first component.</figcaption></figure>

<figure><img src="image_120.png" /><figcaption>(c) The vortex indices of the first component for $N_{min} = 1, N_{max} = 3$.</figcaption></figure>

<figure><img src="image_121.png" /><figcaption>(d) The squared modulus of the second component.</figcaption></figure>

<figure><img src="image_122.png" /><figcaption>(e) The squared modulus of the second component.</figcaption></figure>

<figure><img src="image_123.png" /><figcaption>(f) The vortex indices of the first component for $N_{min} = 1, N_{max} = 4$.</figcaption></figure>

Figure 3.41 – The squared modulus of the first component (a) (resp. (b)), and the second component (d) (resp. (e)) of a minimizer for the energy $E_{\varepsilon, \delta}^\Delta$ (resp. the energy with centrifugal force $\mathcal{E}_{\varepsilon, \delta}^\Delta$) with $\Omega = 5$ and $\delta = 0.998$. On the right panel, we display the vortices indices of the first component (e) and the second component (f) for $\text{tol}_1 = 0.05$ and $\text{tol}_2 = 0.01$.

We can see in Figure 3.41 that we are in the coexistence regime still where each component is disk-shaped. Also in Figure 3.41, we compute the index of the vortices detected. As expected, almost all the indices of the numerical vortices are equal to one which validates numerically that the zeros of the function have a singly quantized phase circulation (except for one vortex in the second component 3.41f which has an index of 2). We can see a stripe pattern in both components (Figures 3.41b, 3.41e), which is in accordance with the theory presented in Chapter 1 Section 4.3, last case.
# 4 Comparison between finite difference and Fourier transform approaches

This section is a brief comparison between the finite difference approach (Section 1 of Chapter 2) and the fast Fourier transform approach (Section 2 of Chapter 2). Although the finite difference approach seemed to be more time efficient, it also has some drawbacks.

In terms of computational time, the finite difference approach is way more efficient when the number of points in the discretization is big ($N \ge 512$). Using predefined sparse matrices, the finite difference approach allows for a faster computation of the minimizer (gradients are faster to compute) than the fast Fourier transformation approach since the number of operations is much less. An example of the time efficiency of the finite difference approach over the FFT is displayed in Figure 3.42a. In this figure, we can see that the computation time of a minimizer of the energy $E_{\epsilon,\delta}^{\Delta}$ using finite difference approach is almost 1.5 times faster with $N = 256$. Notice that the energy $E_{\epsilon,\delta}^{\Delta}$ in Figure 3.42b decreases faster using finite difference approach but then stagnates above the energy level of the FFT approach. This is mainly because the energy using these two approaches is not the same functional. Finally we see in Figure 3.42c that the value of $K$ mainly decreases over the number of iterations to finally reach the stopping criterion ($K^{\Delta} = 0.05$). In this figure we see as well that the gradient method using the finite difference approach needs less iterations to converge than using FFT approach.

In terms of qualitative results, the Fourier transformation approach is more relevant. In particular, it shows better numerical resemblance to the theoretical regimes described in Chapter 1. For instance, comparing Figures 3.17 and 3.38, they both exhibit the coexistence regime and a triangular lattice. The difference is that the triangular lattice using the Fourier transform approach is more distinctive. In addition to that, in the case of a one component condensate with high rotation, we can see a grid orientation effect affecting the results using the finite difference scheme. Indeed, the shape of the hole is different in this regime: it is a square shaped hole in Figures 3.7 and 3.8 (finite difference approach with $\Omega = 85$ and $\Omega = 110$) and it is a disk in Figures 3.28 and 3.29 (fast Fourier transform approach with $\Omega = 50$ and $\Omega = 70$).

<figure><img src="image_124.png" /><figcaption>(a) The computation time (in seconds) as a function of $-\log(K^{\Delta})$.</figcaption></figure>

<figure><img src="image_125.png" /><figcaption>(b) The evolution of energy as a function of the number of iterations.</figcaption></figure>

<figure><img src="image_126.png" /><figcaption>(c) The evolution of the value of $K^{\Delta}$ as a function of the number of iterations.</figcaption></figure>

Figure 3.42 – Data extracted from a simulation of a minimizer of the energy $E_{\epsilon,\delta}^{\Delta}$ for one component condensate using both finite difference and FFT approaches with $N = 256$, $\epsilon = 0.05$ and $\Omega = 5$.
## 5 GPELab comparison

In this section, we study the efficiency of EPG used in this thesis (using Fourier transformation) compared to GPELab (see [8, 9, 11]) for minimizing the energy $E_{\varepsilon,\delta}^{\Delta}$. At first sight, the EPG method is explicit whereas the GPELab method solves a linear system at each time step $\Delta_t$ (here $\Delta_t$ denotes the step of GPELab). We display for each comparison test three convergence criteria (see below) as well as the number of iterations, the energy (as defined in (1.13)) and the execution time. Let $\psi_n$ denote the $n^{th}$ iteration using both algorithms. We compute the following three criteria

1. the difference between two successive iterations $\|\psi_n - \psi_{n-1}\|_{\infty}^2$,

2. the criterion we developed (see equation (2.8)): $K^{\Delta}$,

3. the energy evolution $|E_{\varepsilon,\delta}^{\Delta}(\psi_n) - E_{\varepsilon,\delta}^{\Delta}(\psi_{n-1})|$.

For the simulations with GPELab, we use the first criterion as a convergence test (which is a native criterion in the code). For the simulations with EPG method, we use the first two criteria (in contrast to the previous numerical simulations of Sections 2 and 3 where the only criterion with $K^{\Delta}$ is used). For each test case, we compute the final value of the three criteria and we compare them.

**Remark 24** The energy computed by GPELab's method is different from (1.13). The continuous GPELab's energy goes as follows (see [9, 11] for more details):

$$ \mathcal{E}_{GPE}(u_1, u_2) = \int_D \sum_{\ell=1}^{2} \left[ \frac{1}{2} |\nabla u_{\ell}|^2 + V_{\ell} |u_{\ell}|^2 - \Omega u_{\ell}^* L_z u_{\ell} \right] + \int_D u^* \frac{\beta}{2} \mathbf{F} u, $$

where $V_{\ell} = -\frac{\rho}{2\varepsilon^2}$ is the confinement function associated to component number $\ell$, $\beta = \frac{1}{2\varepsilon^2}$ is the intra-component interaction, $u^*$ is defined as $(\bar{u}_1, \bar{u}_2)$, the matrix $\mathbf{F}$ is given by

$$ \mathbf{F} = \begin{pmatrix} \beta_{1,1} |u_1|^2 + \beta_{1,2} |u_2|^2 & 0 \\ 0 & \beta_{2,1} |u_1|^2 + \beta_{2,2} |u_2|^2 \end{pmatrix}, $$

where $\begin{pmatrix} \beta_{1,1} & \beta_{1,2} \\ \beta_{2,1} & \beta_{2,2} \end{pmatrix} = \begin{pmatrix} 1 & \delta \\ \delta & 1 \end{pmatrix}$ and the rotation operator is defined by $L_z = -i(x\partial_x - y\partial_y)$. The relation between the discrete analogue to $\mathcal{E}_{GPE}$ and $E_{\varepsilon,\delta}^{\Delta}$ (see (2.4)) is the following:

$$ \mathcal{E}_{GPE}^{\Delta} = E_{\varepsilon,\delta}^{\Delta} - \frac{\delta_x \delta_y}{4\varepsilon^2} \sum_{n=0}^{N} \sum_{k=0}^{K} (\rho_{n,k}^2)_{\rho>0}. $$

Observe that the difference between the two energies does not depend on $u = (u_1, u_2)$. In the case of one component condensate (by taking $u_2 = 0$), we take $\beta_{1,1} = 1$ and

$$ \beta_{1,2} = \beta_{2,1} = \beta_{2,2} = 0, \text{ so that the last term of } \mathcal{E}_{GPE} \text{ reads } \int_D u^* \frac{\beta}{2} \mathbf{F} u = \int_D \frac{\beta}{2} |u_1|^4. $$

### 5.1 One component condensate comparison

For this first series of comparisons, we consider a one component Bose-Einstein condensate. We compare both algorithms on the same computer machine, each test done separately.
**Remark on implementation**

We compare in this section our implementation of the EPG method in a ~ 500 lines Python code to that of the implicit gradient method with projection of GPELab All the simulations of this section are numerically simulated on a laptop with an i5 8th generation CPU with 15 Gb of RAM.

**Parameters taken for EPG method**

We take the following parameters for all the one component tests we did. For the initial datum, we choose $\psi^1 = \frac{1}{5}\exp(-10x^2 - 10y^2)$. We choose the confinement parameter $\epsilon = 5 \times 10^{-2}$, the values $N_1 = 1$, $N_2 = 0$ and $K_0 = 10^{-2}$. We refer to Section 1 of this chapter for the rest of the parameters.

**Parameters taken for GPELab**

The equivalent of these parameters in GPELab are the following. We choose Ncomponents= 1, Type='BESP', Delta= 0.5 (the coefficient in front of the kinetic energy), Beta= 200 (equivalent to $\frac{1}{2\epsilon^2}$), the confinement function $V(x, y) = \frac{1}{2\epsilon^2}\min[1, 10(R^2 - x^2 - y^2)]$ (we had to add our own function to the algorithm), $x_{min} = y_{min} = -7$ and $x_{max} = y_{max} = 7$. As for the initial datum, we have to add $\psi^1 = \frac{1}{5}\exp(-10x^2 - 10y^2)$. We also have to modify the normalization step after each time step $\Delta_t$.

**Parameters taken for detecting and computing vortex indices**

For all the vortex indices detection and computations of this section, we use the post processing algorithm for singly quantized vortices (see Section 5.1 of Chapter 2). We choose $N_{min} = 1$, $N_{max} = 5$, $tol_1 = 0.1$ and $tol_2 = 0.02$ to detect most of the vortices.

**Convergence test**

The GPELab algorithm converges whenever the infinity norm of the difference between two successive iteration is less than a certain value $G_0 = \Delta_t \times \text{Stop\_crit}$. It stops without convergence if the number of iterations done exceeds $10^6$.

The EPG method converges if either one of the following statements is true: at an iteration step $n$, if $\|\psi_n^1 - \psi_{n-1}^1\|_\infty < G_0$ or $K^\Delta < K_0$.

**Remark 25** GPELab's minimization method uses a fixed time step $\Delta_t$. The minimizer provided by the method may depend on the value of $\Delta_t$. We make sure we choose $\Delta_t$ small enough so that we converge to a comparable minimizer to that obtained by EPG method from the same initial datum. Another point is that the discretization slightly differs in GPELab from the EPG method so the two methods minimize two different discrete energies. However these two energies are consistent with (1.13).
Low number of points in the discretization $N = 128$

For a first comparison, we set the rotational speed to $\Omega = 1$. We also set $\Delta_t = 10^{-3}$ and $\text{Stop\_crit} = 10^{-2}$ so that $G_0 = 10^{-5}$. We display the results in Figure 3.43.

<figure><img src="image_127.png" /><figcaption>(a) A minimizer of the energy $E_{\varepsilon}^{\Delta}$ using EPG.</figcaption></figure>

<figure><img src="image_128.png" /><figcaption>(b) A minimizer of the energy $E_{GPE}^{\Delta}$ using GPELab.</figcaption></figure>

<table><caption>Figure 3.43 – Comparison between the results of GPELab (right) and EPG (left) for $\Omega = 1$.</caption><thead><tr><th></th><th><i>E</i><sub>ε,δ</sub><sup>Δ</sup></th><th><i>G</i><sub>0</sub></th><th><i>K</i><sup>Δ</sup></th><th>Iteration #</th><th><i>E</i><sub>ε,δ</sub><sup>Δ</sup>(<i>u</i><sub>n</sub>) − <i>E</i><sub>ε,δ</sub><sup>Δ</sup>(<i>u</i><sub>n−1</sub>)</th><th>Time(s)</th></tr></thead><tbody><tr><td>GPELab</td><td>69</td><td>10<sup>−5</sup></td><td>0.189</td><td>16 880</td><td>−10<sup>−4</sup></td><td>2 749</td></tr><tr><td>EPG</td><td>39</td><td>10<sup>−5</sup></td><td>0.045</td><td>116 830</td><td>−1.7 × 10<sup>−4</sup></td><td>2 538</td></tr></tbody></table>

In this first comparison (Figure 3.43), both methods converged with respect to the first criterion ($G_0$) sharing almost the same value for the third criterion (energy difference). Even though the execution time of both methods is approximately the same, the criterion we developed in this thesis (see Section 3 of Chapter 1) is 4 times smaller for the EPG method than for GPELab method.

Notice that the number of vortices in Figure 3.43 varies depending on the method used since each method will converge to a potentially different local minimizer (see also Remark 25). This explains the energy difference as well. Notice also that there is a significant difference between the number of iterations of both methods. It is justified by the fact that GPELab method is implicit and therefore each iteration takes more computational time. In addition to that, GPELab method works best with an odd number of points, this is why it has $2^n + 1$ points per direction (see Section 7 in [9] for reference) in contrast to the EPG method which has $2^n + 2$ points per direction (including the boundary condition). This holds for all the simulations of this section as well.

**Remark 26** Note that we display the third criterion without the absolute value to show that the energy difference between two successive iterations of both methods is negative. This is consistent with the fact that we are using a gradient method.

For a second comparison, we set the rotational velocity to $\Omega = 3$. We also set $\Delta_t = 2 \times 10^{-3}$ and $\text{Stop\_crit} = 25 \times 10^{-4}$ so that $G_0 = 5 \times 10^{-6}$. The results are displayed in Figure 3.44

In this second comparison (Figure 3.44), both methods converge with respect to the first criterion ($G_0$) sharing almost the same value for the third criterion (energy
<figure><img src="image_129.png" /><figcaption>(a) A minimizer of the energy $E_{\varepsilon}^{\Delta}$ using EPG.</figcaption></figure>

<figure><img src="image_130.png" /><figcaption>(b) A minimizer of the energy $E_{\varepsilon}^{\Delta}$ using GPELab.</figcaption></figure>

<table><caption>Figure 3.44 – Comparison between the results of GPELab (right) and EPG (left) for $\Omega = 3$.</caption><thead><tr><th></th><th><i>E</i><sub>ε,δ</sub><sup>Δ</sup></th><th><i>G</i><sub>0</sub></th><th><i>K</i><sup>Δ</sup></th><th>Iteration #</th><th><i>E</i><sub>ε,δ</sub><sup>Δ</sup>(<i>u</i><sub>n</sub>)−<i>E</i><sub>ε,δ</sub><sup>Δ</sup>(<i>u</i><sub>n−1</sub>)</th><th>Time(s)</th></tr></thead><tbody><tr><td>GPELab</td><td>−1427</td><td>5.10<sup>−6</sup></td><td>0.242</td><td>22 591</td><td>−2.3 × 10<sup>−5</sup></td><td>6 742</td></tr><tr><td>EPG</td><td>−1398</td><td>5.10<sup>−6</sup></td><td>0.028</td><td>217 000</td><td>−1.7 × 10<sup>−5</sup></td><td>5 330</td></tr></tbody></table>

difference). We notice that the execution time of GPELab is a bit bigger than that of EPG. On the other hand, the second criterion ($K^{\Delta}$) is 8 times smaller for EPG than for GPELab. We see that the energies of both methods are closer than in the previous experiment, and the numerical minimizers share similarities in the number of vortices and in their positions.

## Moderate number of points in the discretization $N = 256$

For a first comparison, we set the rotational speed to $\Omega = 1$. We also set $\Delta_t = 2 \times 10^{-3}$ and Stop\_crit= $10^{-2}$ so that $G_0 = 2 \times 10^{-5}$. We display the results in Figure 3.45.

In this first comparison with $N = 256$ (Figure 3.45), both methods converge with respect to the first criterion ($G_0$) with GPELab having 3 times less energy difference (third criterion) than EPG. We notice that the execution time of GPELab is almost 4 times greater than EPG and yet the second criterion ($K^{\Delta}$) is 3 times smaller for EPG than for GPELab.

We see more vortices in the numerical minimizers of EPG than in that of GPELab. These vortices are more positioned in the center for GPELab's minimizer and more on the outer disk for EPG's minimizer. All numerical vortices detected for both methods (see Section 5.1 of this chapter for the parameters taken for the index detection algorithm of Section 5.1 of Chapter 2) are singly quantized vortices.
<figure><img src="image_131.png" /><figcaption>(a) A minimizer of the energy $E_{\varepsilon}^{\Delta}$ using EPG.</figcaption></figure>

<figure><img src="image_132.png" /><figcaption>(b) A minimizer of the energy $E_{\varepsilon}^{\Delta}$ using GPELab.</figcaption></figure>

<figure><img src="image_133.png" /><figcaption>(c) The vortex indices of the minimizer obtained by EPG.</figcaption></figure>

<figure><img src="image_134.png" /><figcaption>(d) The vortex indices of the minimizer obtained by 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>38</td><td>2.10<sup>-5</sup></td><td>0.037</td><td>129 569</td><td>-1.4 × 10<sup>-5</sup></td><td>95 189</td></tr><tr><td>EPG</td><td>83</td><td>2.10<sup>-5</sup></td><td>0.011</td><td>271 041</td><td>-4.5 × 10<sup>-5</sup></td><td>24 688</td></tr></tbody></table>

Figure 3.45 – Comparison between the results of GPELab (right) and EPG (left) for $\Omega = 1$ and $N = 256$.

For a second comparison, we set the rotational speed to $\Omega = 3$. We also set $\Delta_t = 2 \times 10^{-3}$ and $\text{Stop\_crit} = 10^{-2}$ so that $G_0 = 2 \times 10^{-5}$. We display the results in Figure 3.46.

In this second comparison with $N = 256$ (Figure 3.46), GPELab converges with respect to the first criterion ($G_0$) while EPG converges with respect to the second criterion ($K_{\Delta}$). GPELab has 3 times less energy difference (third criterion) than EPG on the last iteration. We notice that the execution time of GPELab is almost 4 times greater than that of EPG. However, the second criterion ($K^{\Delta}$) is 5 times smaller for EPG than for GPELab. We see slightly more vortices in the numerical minimizer of GPELab than in that of EPG. The two minimizers share similarities in the vortices positions except in the center. Indeed for EPG, the vortices are aligned on a square lattice around the center and this is not the case for GPELab. However, all numerical vortices detected in the minimizers of both methods are singly quantized.
<figure><img src="image_135.png" /><figcaption>(a) A minimizer of the energy $E_{\epsilon}^{\Delta}$ using EPG.</figcaption></figure>

<figure><img src="image_136.png" /><figcaption>(b) A minimizer of the energy $E_{\epsilon}^{\Delta}$ using GPELab.</figcaption></figure>

<figure><img src="image_137.png" /><figcaption>(c) The vortex indices of the minimizer obtained by EPG.</figcaption></figure>

<figure><img src="image_138.png" /><figcaption>(d) The vortex indices of the minimizer obtained by GPELab.</figcaption></figure>

<table><thead><tr><th></th><th>E<sub>ε,δ</sub><sup>Δ</sup></th><th>10<sup>5</sup>·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>–1426</td><td>2.0</td><td>0.05</td><td>56 507</td><td>–1.3 × 10<sup>–5</sup></td><td>60 209</td></tr><tr><td>EPG</td><td>–1428</td><td>2.4</td><td>0.01</td><td>182 346</td><td>–3.5 × 10<sup>–5</sup></td><td>16 609</td></tr></tbody></table>

Figure 3.46 – Comparison between the results of GPELab (right) and EPG (left) for $\Omega = 3$.
For a final comparison, we set the rotational velocity to $\Omega = 6$. We also set $\Delta_t = 10^{-3}$ and $\text{Stop\_crit} = 10^{-2}$ so that $G_0 = 10^{-5}$. The results are displayed in Figure 3.47.

<figure><img src="image_139.png" /><figcaption>(a) A minimizer of the energy $E_{\varepsilon}^{\Delta}$ using EPG.</figcaption></figure>

<figure><img src="image_140.png" /><figcaption>(b) A minimizer of the energy $E_{\varepsilon}^{\Delta}$ using GPELab.</figcaption></figure>

<figure><img src="image_141.png" /><figcaption>(c) The vortex indices of the minimizer obtained by EPG.</figcaption></figure>

<figure><img src="image_142.png" /><figcaption>(d) The vortex indices of the minimizer obtained by GPELab.</figcaption></figure>

<table><caption>Figure 3.47 – Comparison between the results of GPELab (right) and EPG (left) for $\Omega = 6$.</caption><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>-7334</td><td>10<sup>-5</sup></td><td>0.07</td><td>224 473</td><td>-4.3 × 10<sup>-4</sup></td><td>199 321</td></tr><tr><td>EPG</td><td>-7338</td><td>5.10<sup>-6</sup></td><td>0.006</td><td>707 141</td><td>-1.1 × 10<sup>-5</sup></td><td>64 942</td></tr></tbody></table>

We started this last comparison with $N = 256$ and we noticed that EPG converged fast but with an important energy difference at the last iteration when we used $G_0 = 10^{-5}$. This is why we modify, for this comparison only, the values of the first and third criterion for EPG. The new values for EPG are $G_0 = 5 \times 10^{-6}$ and $K_0 = 5 \times 10^{-3}$. With these parameters, we compare both methods. The results are displayed in Figure 3.47. Both methods converge with respect to the first criterion ($G_0$), with EPG having 40 times less energy difference (third criterion) than that of GPELab. Even though EPG's first criterion of convergence is twice less than GPElabs, we notice that the execution time of GPELab is more than 3 times greater than the one of EPG. However, the second criterion ($K^\Delta$) is 12 times smaller for EPG than for GPELab.
We see almost the same number of vortices in both numerical minimizers of GPELab and EPG. The result of GPELab shows a circular pattern of vortices while the result of EPG is more disrupted. All numerical vortices detected for the minimizers of both methods are singly quantized.

## 5.2 Two components condensate comparison

For the second series of comparison, we consider a two components Bose-Einstein condensate. We also compare both algorithms on the same computer machine, each test done separately.

### Parameters taken for EPG method

We take the following parameters for all the two components tests we did. For the initial data, we choose the same Gaussian function as before (3.2). We choose the confinement parameter $\epsilon = 5 \times 10^{-2}$, the values $N_1 = 0.55$, $N_2 = 0.45$ and $K_0 = 10^{-2}$. We refer to Section 1 of this chapter for the rest of the parameters.

### Parameters taken for GPELab

The equivalent of these parameters in GPELab are the following. We choose Ncomponents= 2, Type='BESP', Delta= 0.5 (the coefficient in front of the kinetic energy), Beta= 200 (equivalent to $\frac{1}{2\epsilon^2}$), Beta_coupled= $\begin{pmatrix} \beta_{1,1} & \beta_{1,2} \\ \beta_{2,1} & \beta_{2,2} \end{pmatrix} = \begin{pmatrix} 1 & \delta \\ \delta & 1 \end{pmatrix}$, the confinement function $V(x, y) = \frac{1}{2\epsilon^2} \min[1, 10(R^2 - x^2 - y^2)]$ (the same function we added before), $x_{min} = y_{min} = -7$ and $x_{max} = y_{max} = 7$. As for the initial data, we have to add $\psi^1 = \psi^2 = \frac{1}{5} \exp(-10x^2 - 10y^2)$. We also have to modify the normalization step after each time step $\Delta_t$.

### Parameters taken for detecting and computing vortex indices

For all the vortex indices detection and computations of this section, we use the post processing algorithm for singly quantized vortices (see Section 5.1 of Chapter 2). We choose $N_{min} = 1$, $N_{max} = 5$, tol$_1$ = 0.05 and tol$_2$ = 0.02 to detect most of the vortices. If large vortices exists with small vortices (like in Figure 3.49), we use a combination of both algorithms (Sections 5.1 and 5.2 of Chapter 2) using $m = 0.4$, $M = 0.6$ and tol$_3$ = 0.3 as parameters for the vortex sheet detection algorithm.

### Numerical results

For a first comparison, we set the rotational speed to $\Omega = 3$ and the interaction strength to $\delta = 0.7$. We also set $\Delta_t = 5 \times 10^{-4}$ and Stop_crit= $2 \times 10^{-2}$ so that $G_0 = 10^{-5}$. We display the results in Figure 3.48. In this comparison of two components condensate in a coexistence regime, both methods converge with respect to the first criterion ($G_0$) with GPELab having 10 times less energy difference (third criterion) than EPG. We notice that the execution time of GPELab is almost 3 times greater than EPG and the second criterion ($K^\Delta$) is 7 times smaller for EPG than for GPELab.
<figure><img src="image_143.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_144.png" /><figcaption>(b) The first and second component of a minimizer of the energy $E_{\varepsilon}^{\Delta}$ using GPELab.</figcaption></figure>

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

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

<table><thead><tr><th></th><th>E<sup>Δ</sup><sub>ε,δ</sub></th><th>G<sub>0</sub></th><th>K<sup>Δ</sup></th><th>Iteration #</th><th>E<sup>Δ</sup><sub>ε,δ</sub>(u<sub>n</sub>) - E<sup>Δ</sup><sub>ε,δ</sub>(u<sub>n-1</sub>)</th><th>Time(s)</th></tr></thead><tbody><tr><td>GPELab</td><td>-2234</td><td>10<sup>-5</sup></td><td>0.07</td><td>206 833</td><td>-2.3 × 10<sup>-6</sup></td><td>163 690</td></tr><tr><td>EPG</td><td>-2227</td><td>10<sup>-5</sup></td><td>0.01</td><td>259 383</td><td>-2 × 10<sup>-5</sup></td><td>50 527</td></tr></tbody></table>

Figure 3.48 – Comparison between the results of GPELab (right) and EPG (left) method for $\Omega = 3$ and $\delta = 0.7$.

We see almost the same number of vortices in both numerical minimizers of EPG and GPELab. These minimizers both display lattices of vortices made of triangles and squares, which is in accordance with the theory presented in Section 4.3 of Chapter 1. All the numerical vortices detected for both methods are singly quantized.

For the second and last comparison test, we set the rotational speed to $\Omega = 3$ and the interaction strength to $\delta = 1.5$. We also set $\Delta_t = 5 \times 10^{-4}$ and $\text{Stop\_crit} = 2 \times 10^{-2}$ so that $G_0 = 10^{-5}$. Due to a symmetrical problem (see Chapter 4 Theorem 41), we display in Figure 3.49 3 results. The first numerical minimizer is computed with EPG using the same symmetrical initial data as (3.2). The second numerical minimizer is computed with EPG using the following non symmetrical initial data

$$ \psi^1 = \frac{1}{5} \exp(-10(x-0.5)^2 - 10(y+0.2)^2), \quad \psi^2 = \frac{1}{5} \exp(-10(x+0.5)^2 - 10(y-0.3)^2). $$

The third numerical minimizer is computed using GPELab with the symmetrical initial data (3.2).

In this comparison of a two components condensate in a segregation regime (Figure 3.49), we compare EPG's symmetrical results (denoted with *) to GPELab's results. Although EPG's non-symmetrical results (denoted with ***) are visually more comparable to that of GPELab, we only compare the efficiency of the two methods starting with the same symmetric initial data.

Both methods (GPELab and EPG *) converge with respect to the first criterion ($G_0$) with GPELab having 37 times less energy difference (third criterion) than EPG.