Guillaume Dujardin

Ingrid Lacroix-Violet

Anthony Nahas

A NUMERICAL STUDY OF VORTEX NUCLEATION IN 2D ROTATING BOSE-EINSTEIN CONDENSATES

Keywords: AMS Classification. 35Q40, 65N35, 65Z05 Bose-Einstein condensation, Gross-Pitaevkii energy, segregation and coexistence regimes, minimization under constraints, gradient algorithm, vortex detection, vortex quantization, numerical experiments

This article introduces a new numerical method for the minimization under constraints of a discrete energy modeling multicomponents rotating Bose-Einstein condensates in the regime of strong confinement and with rotation. Moreover, we consider both segregation and coexistence regimes between the components. The method includes a discretization of a continuous energy in space dimension 2 and a gradient algorithm with adaptive time step and projection for the minimization. It is well known that, depending on the regime, the minimizers may display different structures, sometimes with vorticity (from singly quantized vortices, to vortex sheets and giant holes). In order to study numerically the structures of the minimizers, we introduce in this paper a numerical algorithm for the computation of the indices of the vortices, as well as an algorithm for the computation of the indices of vortex sheets. Several computations are carried out, to illustrate the efficiency of the method, to cover different physical cases, to validate recent theoretical results as well as to support conjectures. Moreover, we compare this method with an alternative method from the literature.

Introduction

Bose-Einstein condensation was predicted by Satyendra Nath Bose [START_REF] Bose | Plancks gesetz und lichtquantenhypothese[END_REF] and Albert Einstein [START_REF] Einstein | Quantentheorie des einatomigen idealen gases[END_REF] in 1924 and 1925. It describes a state of matter in which separate atoms or subatomic particles, cooled to near absolute zero (a few µK), coalesce into a single quantum mechanical entity-that is, one that can be described by a single wave function-on a near-macroscopic scale. Bose-Einstein condensates were first realized experimentally in 1995 [START_REF] Davis | Bose-einstein condensation in a gas of sodium atoms[END_REF][START_REF] Anderson | Observation of Bose-Einstein condensation in a dilute atomic vapor[END_REF]. When the so-called Bose-Einstein condensates (BEC) are set to rotation, topological defects often manifest themselves as vortices that correspond to zeros of the wave function with phase circulation. This phenomenon was first observed in two components BEC [START_REF] Matthews | Vortices in a Bose-Einstein condensate[END_REF]. 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 [START_REF] Ho | Bose-Einstein condensates with large number of vortices[END_REF][START_REF] Mason | Classification of the ground states and topological defects in a rotating twocomponent Bose-Einstein condensate[END_REF]. In the case of 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 [START_REF] Mizushima | Coreless and singular vortex lattices in rotating spinor Bose-Einstein condensates[END_REF], or vortex sheets [START_REF] Mason | Classification of the ground states and topological defects in a rotating twocomponent Bose-Einstein condensate[END_REF][START_REF] Kasamatsu | Vortices in multicomponent Bose-Einstein condensates[END_REF].

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 [START_REF] Ignat | The critical velocity for vortex existence in a two-dimensional rotating bose-einstein condensate[END_REF] (see also references therein). In [START_REF] Aftalion | Vortex patterns and sheets in segregated two component Bose-Einstein condensates[END_REF][START_REF] Aftalion | Thomas-Fermi approximation for coexisting two component Bose-Einstein condensates and nonexistence of vortices for small rotation[END_REF], 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 [START_REF] Bao | Ground-state solution of bose-einstein condensate by directly minimizing the energy functional[END_REF][START_REF] Bao | Ground states and dynamics of multicomponent bose-einstein condensates[END_REF][START_REF] Bao | Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow[END_REF][START_REF] Bao | An efficient and spectrally accurate numerical method for computing dynamics of rotating Bose-Einstein condensates[END_REF][START_REF] Bao | Efficient and spectrally accurate numerical methods for computing ground and first excited states in Bose-Einstein condensates[END_REF][START_REF] Bao | Efficient numerical methods for computing ground states and dynamics of dipolar bose-einstein condensates[END_REF][START_REF] Bao | Ground states of two-component bose-einstein condensates with an internal atomic josephson junction[END_REF][START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF][START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates[END_REF][START_REF] Bao | Computing the ground state and dynamics of the nonlinear schrödinger equation with nonlocal interactions via the nonuniform fft[END_REF], the authors develop numerical methods in several physical contexts and with several space discretizations, which require solving a linear system at each time-step. A theoretical analysis of the convergence of some of these methods has been carried out in several simple contexts (see for example [START_REF] Faou | Convergence of a normalized gradient algorithm for computing ground states[END_REF] for the classical Gross-Pitaevskii energy in one space dimension, for the one component case, without rotation nor confinement, in a neighborhood of the ground state). Another option is to use Sobolev gradients, as opposed to L 2 gradients, as developed in [START_REF] Danaila | A new Sobolev gradient method for direct minimization of the Gross-Pitaevskii energy with rotation[END_REF].

In this article, we study numerically the behaviour of a rotating Bose-Einstein condensate (BEC) in two dimensions in a strong confinement regime based on the numerical minimization of the Gross-Pitaevskii (GP) energy (2.1). We consider the both cases of one component and two components condensates. For two components condensates, we consider both segregation and coexistence regimes between the two components. We introduce a new discretization for the GP energy using the Fast Fourier transformation (FFT) scheme. We use an explicit L 2 gradient method with adaptive step and projection for the minimization of our discrete energy under constraints. We derive a stopping criterion based on the evaluation of the gradient of the energy on the constrained manifold. Moreover we provide a post processing algorithm for computing the indices of the vortices and of the vortex sheets of the minimizers. Finally, we compare the efficiency of our explicit gradient method with projection, named EPG, to the method used in GPELab ( [START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations I: computation of stationary solutions[END_REF][START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations II: Dynamics and stochastic simulations[END_REF]). The numerical results of the test cases in this paper will illustrate how an explicit projected gradient method together with an energy discretization allowing for the use of FFT in the computation of its gradient makes it possible to outperform (linearly) implicit methods such as that of GPELab.

The outline of this paper is as follows. We introduce the continuous model for rotating Bose-Einstein condensates with two components in Section 2 and we recall the different regimes for one component and two components condensates. In order to numerically study these regimes, we discretize in Section 3 the continuous Gross-Pitaevskii energy. Our discretization uses the discrete Plancherel formula and allows for using FFT in the computation of the gradient. We introduce in this section a gradient method for the minimization of the energy with 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 develop in the same section. We develop in Section 4 two post processing algorithms for the computation of indices. The first deals with vortices and the second deals with vortex sheets. We present numerical results in Section 5 for the regimes described in Section 2. In particular, we validate numerically recent theoretical results and we support some conjectures as for example the existence of vortex sheets in a segregation regime. Section 6 is devoted to the comparison of the efficiency of the EPG method with that of GPELab [START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations I: computation of stationary solutions[END_REF][START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations II: Dynamics and stochastic simulations[END_REF].

2. The model and the regimes 2.1. The model. We consider the following model for the energy of a two components rotating Bose-Einstein condensate in dimension 2 in the limit of strong confinement and strong rotation studied in [START_REF] Aftalion | Vortex patterns and sheets in segregated two component Bose-Einstein condensates[END_REF]. This model reads

E Ω ε,δ pu 1 , u 2 q " 2 ÿ "1 1 2 ż D › › ∇u ´iΩu x K › › 2 dxdy `Wε,δ pu 1 , u 2 q,
where D Ă R 2 is the bounded physical domain of interest, Ω P R is the rotation speed, ε and δ are positive constants, x K " p´y, xq, u 1 , u 2 P H 1 pD, Cq are the wave functions related to each component of the condensate. Moreover, W ε,δ is the confining part of the energy defined as

W ε,δ pu 1 , u 2 q " 1 4ε 2 ż D pρprq ´|u 1 | 2 q 2 dxdy `1 4ε 2 ż D pρprq ´|u 2 | 2 q 2 dxdy `δ 2ε 2 ż D |u 1 | 2 |u 2 | 2 dxdy ´1 4ε 2 ż D ρ 2 prqdxdy,
where ρ is a function of r " a x 2 `y2 to be defined later (ρ " 1 in [START_REF] Aftalion | Vortex patterns and sheets in segregated two component Bose-Einstein condensates[END_REF]), δ ą 0 measures the stength of the interaction between the two components, 1{ε 2 measures the strength of the internal interaction in each component of the condensate. In this paper, we take D as the disk of radius R ą 0 centered at the origin of R 2 .

Our aim is to compute, in the regime ε Ñ 0, minimizers of the energy E Ω ε,δ in H 1 0 pD, Cq with the constraints (1)

ż D |u 1 | 2 dxdy " M N 1 , and
ż D |u 2 | 2 dxdy " M N 2 ,
where M " ş D ρprqdxdy, and N 1 , N 2 ě 0 with N 1 `N2 " 1. Note that, in the sequel, δ ą 0 and Ω P R may or may not depend on ε when considering a sequence of minimizing problems with smaller and smaller ε.

Note that one has also

(2) W ε,δ pu 1 , u 2 q " 1 4ε 2 ż D pρprq ´|u 1 | 2 ´|u 2 | 2 q 2 dxdy `δ ´1 2ε 2 ż D |u 1 | 2 |u 2 | 2 dxdy.
In view of this expression of W ε,δ , we refer to the case δ ą 1 as the segregation regime since minimizing the energy tends to split the two components in this case, and to the case δ ă 1 as the coexistence regime since minimizing the energy tends to mix the two components of the condensate in this other case.

Remark 1. A straightfoward computation yields for all P t1, 2u,

1 2 ż D › › ∇u ´iΩu x K › › 2 dxdy " 1 2 ż D }∇u } 2 dxdy ´Ω ż D `iu x K ¨∇u ˘dxdy `Ω2 2 ż D r 2 |u | 2 dxdy. (3)
Following Remark 1 of [START_REF] Aftalion | Vortex patterns and sheets in segregated two component Bose-Einstein condensates[END_REF] and the one component analysis carried out in [START_REF] Correggi | Critical rotational speeds in the Gross-Pitaevskii theory on a disc with Dirichlet boundary conditions[END_REF], in the regime |Ω| ăă 1{ε, the contribution of the third term in (3) plays no role in the asymptotics ε Ñ 0. Therefore, from now on and unless stated otherwise, we shall consider in the regime |Ω| ăă 1{ε the minimization of the energy

E Ω ε,δ pu 1 , u 2 q " 2 ÿ "1 " 1 2 ż D }∇u } 2 dxdy ´Ω ż D `iu x K ¨∇u ˘dxdy  `Wε,δ pu 1 , u 2 q (4) " E Ω ε,δ pu 1 , u 2 q ´Ω2 2 ż D r 2 p|u 1 | 2 `|u 2 | 2 qdxdy.
2.2. The different regimes. We recall in this section the theoretical results obtained in [START_REF] Correggi | Critical rotational speeds in the Gross-Pitaevskii theory on a disc with Dirichlet boundary conditions[END_REF][START_REF] Aftalion | Vortex patterns and sheets in segregated two component Bose-Einstein condensates[END_REF][START_REF] Goldman | Phase segregation for binary mixtures of Bose-Einstein condensates[END_REF][START_REF] Kasamatsu | Vortex phase diagram in rotating two-component Bose-Einstein condensates[END_REF][START_REF] Mason | Classification of the ground states and topological defects in a rotating twocomponent Bose-Einstein condensate[END_REF]. These results describe the asymptotics of the minimizers of (4) in different regimes. These regimes depend on the behavior of the rotation speed Ω in the strong confinement limit (ε Ñ 0) in one component condensates (see 2.2.1). They also depend on the value of the segregation/coexistence parameter (δ ě 1 or δ ă 1) in two components condensates (see 2.2.2 and 2.2.3).

2.2.1. One component condensates. According to [START_REF] Correggi | Critical rotational speeds in the Gross-Pitaevskii theory on a disc with Dirichlet boundary conditions[END_REF], for one component condensates (N 1 " 1, N 2 " 0), the behaviour of the structure of the vortices of the minimizers of (4) with constraints (1) in the regime ε Ñ 0 changes with the dependency of the rotation speed Ω " Ω ε with respect to ε. Namely, there exists three critical rotation speeds Ω

1 ε ăă Ω 2 ε ăă Ω 3 ε such that ' if Ω ε ă Ω ε,1 , then minimizers u ε 1 have no vortices, ' if Ω 1 ε ă Ω ε ă Ω 2 ε
, then the vortices of minimizers u ε 1 appear on a hexagonal lattice and are singly quantized,

' if Ω 2 ε ă Ω ε ă Ω 3
ε , the centrifugal force comes into play: a hole appears in minimizers u ε 1 and the condensate looks like an annulus, with vortices located on a lattice on the annulus, ' if Ω 3 ε ă Ω ε , 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

(5) Ω 1 ε " logp1{εq, Ω 2 ε " 1 ε , and Ω 3 ε " 1 ε 2 logp1{εq
.

Note that δ plays no role in this case since u 2 " 0 because N 2 " 0. According to Remark 1, the regimes Ω ε ă Ω 2 ε are similar for minimizers of E Ωε ε,δ and E Ωε ε,δ . In the regimes Ω 2 ε ă Ω ε , the annulus behaviour is due to the centrifugal force in E Ωε ε,δ . One of the goals of this paper is to confirm numerically the four regimes above (described in [START_REF] Correggi | Critical rotational speeds in the Gross-Pitaevskii theory on a disc with Dirichlet boundary conditions[END_REF]). To do so, we introduce a numerical algorithm in Section 3 and we use it to provide numerical simulations in Section 5.1.

2.2.2.

Two components condensates in the segregation regime (δ ą 1). In the recent paper [START_REF] Aftalion | Vortex patterns and sheets in segregated two component Bose-Einstein condensates[END_REF], the authors consider two components condensates in the segregation case δ ą 1. For N 1 P p0, 1q (recall that N 2 " 1 ´N1 ), they introduce the minimizing perimeter

N 1 " min ωĂD |ω|"N 1 perpωq,
and they perform an analysis of the minimizers of E Ωε ε,δε in the regime ε Ñ 0 depending on whether Ω ε " 0 or not.

First, they address the minimization of (4) with constraints (1) when Ω " 0. In this case, the squared moduli of the minimizers u ε 1 and u ε 2 tend to 1 in two separate regions of D and the authors of [START_REF] Aftalion | Vortex patterns and sheets in segregated two component Bose-Einstein condensates[END_REF] prove that their sum

v 2 ε " |u ε 1 | 2 `|u ε 2 | 2 and the normalized energy εE 0 ε,δε
have the following behaviour ' in the regime δ ε ε 2 Ñ `8 (strong segregation regime [START_REF] Aftalion | A minimal interface problem arising from a two component Bose-Einstein condensate via Γ-convergence[END_REF]), then inf D v 2 ε tends to 0 and εE 0 ε,δε tends to some constant times N 1 ; ' in the regime ε Ñ 0 with δ ą 1 fixed, then inf D v 2 ε tends to some number between 0 and 1 and εE 0 ε,δε tends to some constant (which depends on δ) times N 1 [START_REF] Goldman | Sharp interface limit for two components Bose-Einstein condensates[END_REF]; ' in the regime δ ε Ñ 1 with ε " ε{ ? δ ε ´1 Ñ 0, the rescaled energy εE 0 ε,δε tends to N 1 {2 and it is expected that inf D v ε tends to 1 [START_REF] Goldman | Phase segregation for binary mixtures of Bose-Einstein condensates[END_REF]. Second, they address the minimization of (4) with constraints (1) when Ω " Ω ε Ñ `8 as ε tends to 0. With ε defined as above, they consider a regime where δ " δ ε Ñ 1 and ε Ñ 0 as ε tends to 0 and they prove that ' there exists two constants C 1 , C 2 ą 0 such that, for all i, if Ω ε ą C i logp1{εq, then the infimum of the limiting density |u ε i | 2 vanishes as ε tends to 0; ' in the regime of moderate rotational speed logp1{εq ăă Ω ε ăă 1{ε, 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{pε logp1{εqq ăă Ω ε ăă 1{ε 2 , the leading order in the energy E Ω ε ε,δε is the vortex energy, and the authors of [START_REF] Aftalion | Vortex patterns and sheets in segregated two component Bose-Einstein condensates[END_REF] conjecture the possibility of observing vortex sheets.

2.2.3.

Two components condensates in the coexistence regime (δ ă 1). In the regime δ ă 1, the two components still have mass in the region where ρ ą 0, but the supports of the two components tend to overlap (see [START_REF] Kasamatsu | Vortex phase diagram in rotating two-component Bose-Einstein condensates[END_REF][START_REF] Mason | Classification of the ground states and topological defects in a rotating twocomponent Bose-Einstein condensate[END_REF]). In this regime (δ ă 1), depending on the rotational speed Ω, we should observe four different qualitative behaviours (see [START_REF] Kasamatsu | Vortex phase diagram in rotating two-component Bose-Einstein condensates[END_REF][START_REF] Mason | Classification of the ground states and topological defects in a rotating twocomponent Bose-Einstein condensate[END_REF]) for the minimizers when ε tends to 0:

' The first behaviour (δ P r0, 1q), with a very low velocity Ω, is when the two components coexist and there are no vortices: each component has mass in the disc where ρ ą 0 and the profile of the components depends on N 1 and N 2 . ' The second behaviour, when δ " 0 (which means there is no interactions between the two components), we should observe, depending on the (high) rotation speed Ω, the existence of a triangular vortex lattice. ' As δ P p0, 1q 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 this case, after a certain value of δ, the vortices in each component form a square lattice. ' The last behaviour corresponds to δ Ñ 1 (when ε Ñ 0). In that case, we should observe either stripe or double-core vortex lattices in the minimizers.

The discretisation of the problem and the minimization method

Minimizing the energy E Ω ε,δ defined in (4) under the constraints (1) is a continuous minimization problem over an infinite dimensional manifold that we replace in this section by a discrete minimization problem over a finite dimensional manifold. We consider a gradient descent method with projection over the constrained manifold, which is a first-order iterative optimization algorithm for finding a local minimum of a differentiable function. We name this method EPG. By calculating the gradient of the energy E Ω ε,δ , we are able to approach more and more towards the solution of the discrete problem after each iteration of the method. In Section 3.1, we describe the discretization of the energy E Ω ε,δ defined in [START_REF] Anderson | Observation of Bose-Einstein condensation in a dilute atomic vapor[END_REF]. Then, we compute its gradient in Section 3.2. We establish a criterion for the minimization of the energy on the finite dimensional manifold in Section 3.3. We propose a full description of the minimization method in Section 3.4 as a conclusion.

3.1. Discretization of the energy. For the discretization of the energy defined in (4), we decide to include the disk D of radius R centered at the origin of R 2 in a larger square of side length 2L centered at the origin. In addition, we extend the function ρ used to define the energy in (4) to the square in such a way that it takes (very) negative values outside D. By imposing homogeneous Dirichlet conditions on the boundary of the square, such an extension of ρ should, in the regime ε Ñ 0, lead to minimizers with very small squared modulus outside the disk, the restriction to the disk of which can be considered as approximations to the continuous minimizers of E Ω ε,δ defined in (4). This is due in particular to the first term in the definition of W ε,δ in [START_REF] Aftalion | A minimal interface problem arising from a two component Bose-Einstein condensate via Γ-convergence[END_REF].

Therefore, we choose L ą R so that the disc D of radius R centered at the origin of R 2 is included in the square r´L, Ls 2 . Then, we discretize r´L, Ls 2 with N `2 equidistant points with respect to the x-axis and N `2 equidistant points with respect to the y-axis for some integer N . We do this by setting δ x " δ y " 2L N `1 and x n " ´L `nδ x , and y k " ´L `kδ y ,

for n, k P t0, ..., N `1u.

We use the letter ψ to denote the discrete counterpart to the continuous wave functions denoted by u in Section 2. As explained above, we think of ψ n,k as an approximation of u px n , y k q for " 1, 2. In matrix form, we use the notation:

ψ " ¨ψ 0,0 ψ 0,1 ... ψ 0,N `1 ψ 1,0 ψ 1,1 ... ψ 1,N `1 . . . . . . . . . . . . ψ N `1,0 ψ N `1,1 ... ψ N `1,N `1‹ ‹ ‹ ' P C pN `2q 2 ,
with the convention that ψ n,k " 0 if either n " 0 or k " 0 or n " N `1 or k " N `1. In order to separate real and imaginary parts of the unknowns, we set ψ n,k " p n,k `iq n,k for " 1, 2, and set accordingly ψ " P `iQ . This allows us to write functions of the complex-valued variables ψ as functions of the real-valued variables P and Q . We decide to use the discrete Fourier transform and its inverse for the discretization of the terms in the energy E Ω ε,δ defined in (4) that involve gradients. This choice allows the use of Fast Fourier Transform algorithms for the computation of the gradient of the discrete energy (see section 3.2). To do so, we set the following definitions, which are related to that of Python Numpy. For v P C pN `2q 2 and for n, p, k, q P t0, . . . , N `1u, the discrete Fourier transforms of v in the x-and y-direction and their inverses are given by:

' `fft x pvq ˘p,k " ř N `1 m"0 v m,k e iπm e ´2πi mp N `1 , ' `ifft x pvq ˘n,k " 1 N `1 e ´iπn ř N `1 p"0 v p,k e 2πi np N `1 ,
' `fft y pvq ˘n,q " ř N `1 l"0 v n,l e iπl e ´2πi lq N `1 , ' `ifft y pvq ˘n,k " 1 N `1 e ´iπk ř N `1 q"0 v n,q e 2πi kq N `1 .

Remark 2. Note that, with the definitions above, the discrete Fourier transforms are not inverses by pairs. Indeed, we have for all v, ifft x pfft x qpvq " v `Opδxq and ifft y pfft y qpvq " v `Opδyq. We will not need any such inversion formula in this paper however.

For the discretization of the rotational energy, we define X, Y, Ξ, Λ P R pN `2q 2 for all n, k P t0, . . . , N `1u by where δ ξ " δ λ " π L . Observe that all the columns of the matrix Ξ are equal to the vector ξ " pξ k q 0ďkďN `1 and all the rows of the matrix Λ are equal to the vector λ " pλ n q 0ďnďN `1. Moreover, we denote by ˚the term by term multiplication operator between 2 matrices.

X n,k " x n , Y n,k " y k , Ξ n,k " ξ n ,
We define E ∆ ε,δ pψ 1 , ψ 2 q as a discrete counterpart of E Ω ε,δ pu 1 , u 2 q by setting for

ψ 1 , ψ 2 P C pN `2q 2 (6) E ∆ ε,δ pψ 1 , ψ 2 q " ÿ "1,2 ´pE cin q ∆ pψ q `pE rot q ∆ ε pψ q ¯`pE W q ∆ ε,δ pψ 1 , ψ 2 q,
where the superscript ∆ indicates that these energies depend on discrete variables. Let us describe the discretized energy terms which appear above:

' pE cin q ∆ corresponds to the kinetic energy and is defined by

pE cin q ∆ pψ q " δ 2 x 2 N `1 ÿ n,k"0
ˆˇˇi fft x `iΞ ˚fft x pψ q ˘ˇˇ2 `ˇˇi fft y `iΛ ˚fft y pψ q ˘ˇˇ2 ˙n,k , ' pE rot q ∆ ε corresponds to the rotational energy and is defined by

pE rot q ∆ ε pψ q " ´Ωε δ 2 x N `1 ÿ n,k"0
´´iψ ˚" ´Y ˚ifft x `iΞ ˚fft x pψ q ˘`X ˚ifft y `iΛ ˚fft y pψ q ˘ı¯n ,k , ' pE W q ∆ ε,δ corresponds to the confinement energy and is defined by

pE W q ∆ ε,δ pψ 1 , ψ 2 q " δ 2 x 4ε 2 N `1 ÿ n,k"0 ˆρpr n,k q ´|ψ 1 n,k | 2 ´|ψ 2 n,k | 2 ˙2 `δ2 x pδ ´1q 2ε 2 N `1 ÿ n,k"0 |ψ 1 n,k | 2 |ψ 2 n,k | 2 , where r n,k " b x 2 n `y2 k .
3.2. Computation of the gradient of the discrete energy [START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates[END_REF]. Each energy term in ( 6) is a function of the 4N 2 real variables pP 1 n,k q 1ďn,kďN , pP 2 n,k q 1ďn,kďN , pQ 1 n,k q 1ďn,kďN and pQ 2 n,k q 1ďn,kďN using the definitions from Section 3.1. Let us compute the gradient of each of these energy terms (still using the convention p n,k " q n,k " 0 if n " 0 or k " 0 or n " N `1 or k " N `1) with respect to these variables for the usual scalar product on R 4N 2 . Let V conf P R pN `2q 2 be defined as:

V conf " ¨ρpr 0,0 q ρpr 0,1 q . . . ρpr 0,N `1q ρpr 1,0 q ρpr 1,1 q . . . ρpr 1,N `1q . . .

ρpr N `1,0 q ρpr N `1,1 q . . . ρpr N `1,N `1q ‹ ‹ ‹ ' .
For P t1, 2u we have for the discrete kinetic energy:

BpE cin q ∆ BP pψ q " δ 2 x ˜"ifft x `Ξ˚2 ˚fft x pψ q ˘ın,k `"ifft y `Λ˚2 ˚fft y pψ q ˘ın,k ¸1ďn,kďN ,
BpE cin q ∆ BQ pψ q " δ 2 x ˜"ifft x `´iΞ ˚2 ˚fft x pψ q ˘ın,k `"ifft y `´iΛ ˚2 ˚fft y pψ q ˘ın,k ¸1ďn,kďN .
For P t1, 2u we have for the discrete rotational energy:

BpE rot q ∆ ε BP pψ q " ´2Ω ε δ 2 x ˜"X ˚ifft y `Λ ˚fft y pψ q ˘ın,k ´"Y ˚ifft x `Ξ ˚fft x pψ q ˘ın,k ¸1ďn,kďN , BpE rot q ∆ ε BQ pψ q " ´2Ω ε δ 2 x ˜´i " X ˚ifft y `Λ ˚fft y pψ q ˘ın,k `i " Y ˚ifft x `Ξ ˚fft x pψ q ˘ın,k ¸1ďn,kďN .
For P t1, 2u we have for the discrete interaction energy:

BpE W q ∆ ε,δ BP pψ 1 , ψ 2 q " ´δ2 x ε 2 ˆ"P ˚´V conf ´ˇψ 1 ˇˇ2 ´ˇψ 2 ˇˇ2
¯ın,k `p1 ´δq

" P ˚|ψ 3´ | 2 ı n,k ˙1ďn,kďN , BpE W q ∆ ε,δ BQ pψ 1 , ψ 2 q " ´δ2 x ε 2 ˆ"Q ˚´V conf ´ˇψ 1 ˇˇ2 ´ˇψ 2 ˇˇ2 ¯ın,k `p1 ´δq " Q ˚|ψ 3´ | 2 ı n,k ˙1ďn,kďN .

3.3.

A criterion for the minimization of E ∆ ε,δ under constraints. We identify C 2N 2 with the subspace of M N `2pCq 2 consisting in pairs of matrices with zero first and last row and column. Let us denote the usual scalar product on

C 2N 2 by xu, vy " ¨N ÿ n,k"1 u n,k Ě v n,k ',
and we set }v} 2 ∆ " δ 2 x xv, vy. We replace the continuous constraints (1) with the discrete analogues

(7) }ψ } 2 ∆ " δ 2 x N `1 ÿ k"0 N `1 ÿ n"0 |ψ k,n | 2 " M N p " 1, 2q.
In order to derive a stopping criterion for the minimization of E ∆ ε,δ defined in (6) under the constraints [START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations II: Dynamics and stochastic simulations[END_REF], we prove the following proposition. Proposition 3. If pψ 1˚, ψ 2˚q P C pN `2q 2 are minimizers of E ∆ ε,δ under the constraints (7) satisfying the homogeneous Dirichlet boundary conditions, then the four N ˆN real-valued matrices

K ∆ 1,P , K ∆ 2,P , K ∆ 1,Q , K ∆ 2,Q defined by 1 K ∆ ,P " BE ∆ ε,δ BP pψ 1˚, ψ 2˚q ´δ2 x }ψ ˚}2 ∆ N ÿ n,k"1 ˆBE ∆ ε,δ BP pψ 1˚, ψ 2˚q ˚P ˚`BE ∆ ε,δ BQ pψ 1˚, ψ 2˚q ˚Q ˚˙n,k P ˚,
and

K ∆ ,Q " BE ∆ ε,δ BQ pψ 1˚, ψ 2˚q ´δ2 x }ψ ˚}2 ∆ N ÿ n,k"1 ˆBE ∆ ε,δ BP pψ 1˚, ψ 2˚q ˚P ˚`BE ∆ ε,δ BQ pψ 1˚, ψ 2˚q ˚Q ˚˙n,k Q ˚,
for P t1, 2u, vanish.

Proof. Assume that pψ 1˚, ψ 2˚q P C pN `2q 2 are minimizers of E ∆ ε,δ under the constraints (7) satisfying the homogeneous Dirichlet boundary conditions. For all ψ 1 , ψ 2 P C pN `2q 2 satisfying the homogeneous Dirichlet boundary conditions too and for all t P R small enough, we define ( 8)

f ∆ ptq " E ∆ ε,δ ˆψ1˚`t ψ 1 }ψ 1˚`t ψ 1 } ∆ a N 1 M , ψ 2˚`t ψ 2 }ψ 2˚`t ψ 2 } ∆ a N 2 M ˙.
The fact that pψ 1˚, ψ 2˚q are minimizers of E ∆ ε,δ under the constraints [START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations II: Dynamics and stochastic simulations[END_REF] implies that (9) pf ∆ q 1 p0q " 0, for all pψ 1 , ψ 2 q. In order to rewrite ( 9) as a criterion depending on pψ 1˚, ψ 2˚q , we set ψ ˚"

P ˚`iQ ˚and ψ " P `iQ with " 1, 2 using the definitions from Section 3.1. Let us denote for all " 1, 2 and small t P R: g ∆ ,P ptq "

P ˚`tP }ψ ˚`tψ } ∆ a N M , and g ∆ ,Q ptq " Q ˚`tQ }ψ ˚`tψ } ∆ a N M ,
so that, using [START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations II: Dynamics and stochastic simulations[END_REF], we have

› › g ∆ ,P ptq `ig ∆ ,Q ptq › › 2 ∆ " N M. Introducing the notation, a " δ 2 x }ψ ˚}2 ∆ ¨N ÿ n,k"1 pp ln ,k p l n,k `qln ,k q l n,k q '" δ 2 x }ψ ˚}2 ∆ pxP ˚ , P y `xQ ˚ , Q yq ,
we infer pg ∆ ,P q 1 p0q " P ´a P ˚, and

pg ∆ ,Q q 1 p0q " Q ´a Q ˚.
Using the definition (8) of f ∆ and the chain rule, we obtain

pf ∆ q 1 p0q " 2 ÿ "1 ˜C BE ∆ ε,δ BP pψ 1˚, ψ 2˚q , pg ∆ ,P q 1 p0q G `C BE ∆ ε,δ BQ pψ 1˚, ψ 2˚q , pg ∆ ,Q q 1 p0q G" 2 ÿ "1 N ÿ n,k"1 ˆBE ∆ ε,δ BP pψ 1˚, ψ 2˚q ˚`P ´a P ˚˘`B E ∆ ε,δ BQ pψ 1˚, ψ 2˚q ˚`Q ´a Q ˚˘˙n ,k
, 1 Recall that ˚is the term-by-term product defined in Section 3.1.

As expected, pf ∆ q 1 p0q is an R-linear form depending on pP 1 , Q 1 , P 2 , Q 2 q. The expression above implies that (10)

pf ∆ q 1 p0q " N ÿ n,k"1 `K∆ 1,P ˚P 1 `K∆ 1,Q ˚Q1 `K∆ 2,P ˚P 2 `K∆ 2,Q ˚Q2 ˘n,k ,
with K ∆ ,P and K ∆ ,Q the matrices defined in the statement of Proposition 3. This concludes the proof.

Therefore, in the method described in Section 3.4, we will use the fact that

(11) K ∆ " 2 ÿ "1 }K ∆ ,P } ∆ `2 ÿ "1 }K ∆ ,Q } ∆ ,
is below some threshold as a criterion for numerical convergence.

3.4. Gradient Method. In order to minimize the discretized energy defined in ( 6) under the constraints ( 7), we propose the following gradient method with projection and adaptive step.

' Initialization:

-Choose ψ 1,1 , ψ 2,1 P C pN `2q 2 , normalized as: }ψ ,1 } ∆ " ? N M , p " 1, 2q and satisfying the homogeneous Dirichlet boundary conditions, choose a step h ą 0, choose two tolerance parameters h 0 , K 0 ą 0 for the convergence test. ' Iteration:

(1) Compute the gradient:

∇E ∆ ε,δ `ψ1,m , ψ 2,m ˘.
(2) Compute the auxiliary step

ˆψ 1,m ψ2,m
˙" ˆψ1,m ψ 2,m ˙´h∇E ∆ ε,δ pψ 1,m , ψ 2,m q, and set the homogeneous Dirichlet boundary conditions on ψ1,m and ψ2,m . (3) Normalize the auxiliary step to obtain an attempt for the next step:

ψ ,m`1 " ψ ,m } ψ ,m } ∆ a N M , " 1, 2.
(

) If E ∆ ε,δ pψ 1,m`1 , ψ 2,m`1 q ą E ∆ ε,δ pψ 1,m , ψ 2, 4
m q, then we replace h by h{2 (provided h{2 ě h 0 , otherwise we stop without convergence) and we compute a new auxiliary step p ψ1,m , ψ2,m q for the same m by going back to step 2. Otherwise, we compute K ∆ at point pψ 1,m`1 , ψ 2,m`1 q using (11). If K ∆ ě K 0 , then we replace m by m `1 and we restart at step 1. Otherwise, we take pψ 1,m`1 , ψ 2,m`1 q as an approximate minimizer and we stop with convergence.

Post-processing algorithms

4.1. Indices of the Vortices. In this section, we introduce an algorithm that we developed for the computation of the indices of the vortices in the minimizers of the discrete energy functional defined in [START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates[END_REF] under the constraints [START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations II: Dynamics and stochastic simulations[END_REF], obtained using the minimization method described in Section 3. This algorithm is used in Section 5.

The algorithm relies on 4 numerical parameters tol 1 ą 0, tol 2 ą 0, N min P N ˚and N max P N ˚with N min ď N max . It follows the 4 steps below. The three first steps identify the vortices' centers and the last step computes the vortices' indices. In this section, ψ denotes the square complex matrix with N 2 entries for either of the two components of the Bose-Einstein condensate.

Step 1: We determine the potential centers of the vortices and establish a list of candidates as the set of pn, kq P t1, . . . , N u 2 such that |ψ n,k | 2 ă tol 1 .

Step 2: We build a second list P based on the first list above using the following rule. For each potential center pn, kq in the list established in Step In other words, we have determined a list P of couple of points pn, kq satisfying the following conditions:

' |ψ n,k | 2 ă tol 1 , ' |ψ i,j | 2 ą |ψ n,k | 2 `tol 2
, for all the couples pi, jq in S λpn,kq pn, kq.

Step 3: We consider each center pn, kq from the list P and we identify if another center is inside the square Ť 1ďλďλpn,kq S λ pn, kq. If this is the case, we eliminate the center with the biggest |ψ| 2 at the center from the list. We repeat this step until we are left with isolated centers. Let us denote by T the list of all the couples pn, kq corresponding to isolated centers.

Step 4: We compute the indices using the following rules. For each pn, kq P T, we start off with any couple pi, jq in S λpn,kq pn, kq. Then we compute their associated angle θ 0 " argpψ i,j q. Note that there are 8λpn, kq couples in S λpn,kq pn, kq (see Figure 1). After computing the first angle θ 0 , we proceed to compute the other angles θ 1 , θ 2 , ¨¨¨, θ 8λpn,kq in the following way:

' After computing the angle θ m , the next angle θm`1 is computed as an argument of the next value of ψ on the square S λpn,kq pn, kq with anticlockwise orientation (see figure 1). ' We set θ m`1 :" θm`1 `2kπ with k " argmin lPZ | θm`1 ´θm `2πl|. Eventually, index of the vortex pn, kq is equal to pθ 8λpn,kq ´θ0 q{2π. 4.2. Indices of the vortex sheets. In this section, 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 defined in (6) under the constraints (7) using the minimization method described in Section 3. Let us note that this algorithm requires some input from a human being at some point. This algorithm is used in Section 5.

The algorithm relies on 5 main steps. It uses 3 numerical values: m ď 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 splits these contours into discrete connex curves. The third step requires a human being to decide whether each contour should be discarded or merged with another one. The fourth step sorts the optimized contours, and provides them with an orientation. The fifth and last step computes each vortex sheet's index alongside its contour. Step 1: The contour of a vortex sheet in either 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 other way around). The first step consists in finding regions on the discrete square r´L, Ls 2 where this occurs. We identify, for each component ψ , the coordinates on the grid for the values |ψ | 2 between m and M . Next, we add the grid coordinates close to the circle BD, of which we retain only the coordinates where |ψ | 2 ď tol 3 . These grid points constitute the union of potential contours of vortex sheets (see Figure 2A for an example). Let us sort these points, to ease the search for connex components in contours in the next step. Let us note K the set of coordinates we have found so far. We are looking for a union of connex curves that describe the 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 limit the number of neighbours to any point on the curve to 3. For each grid point of coordinates pi, jq, we denote by S 1 pi, jq the set of its neighbours. For all grid point pi, jq P K, we act as follows:

(1) If pi, jq has no neighbours, we just remove it from the set K.

(2) If pi, jq has 1 neighbour, we add to the set K a couple pi 1 , j 1 q R K verifying pi 1 , j 1 q " argmin pn,kqPS 1 pi,jq ˇˇ|ψ i,j | 2 ´|ψ n,k | 2 ˇˇ.

(3) If pi, jq has 2 or 3 neighbours, we do nothing. (4) If pi, jq has at least 4 neighbours, we remove the couple pi 1 , j 1 q from K verifying: Step 2: We partition K into the union of a finite number of contours K p corresponding to connex curves. We start with p " 1. At the p th vortex sheet, we create dynamically a new list K p " H. We start by adding to K p a grid point at random from the list K and we remove it from K. Then, we add all of its neighbours one by one to K p and we remove them from K. We repeat this until we are left with no neighbours in K to all the grid points in K p . Then, as long as there are points in K, we increase p to p `1 and restart this step again until K is empty.

pi 1 , j 1 q " argmax pn,kqPS 1 pi,jqXK ˇˇ|ψ i,j | 2 ´|ψ n,k | 2 ˇˇ,
Step 3: After dividing the contours of the vortex sheets into discrete connex curves K p , we allow for human intervention. First, a human being decides if each K p should be considered or ignored. Second, a human being decides whether couples of two discrete connex curves K p and K p 1 should be merged or not. To merge two discrete connex curves K p and K p 1 , we follow these steps. First, we search for the closest two couples pn, kq and pn 1 , k 1 q from K q and K q 1 (using the usual distance in R 2 ) and we compute their middle pi 1 , j 1 q " pt n`n 1 2 u, t k`k 1 2 uq. If this middle has neighbours in both categories, we merge them into a new set K q Y K q 1 Y pi 1 , j 1 q. Otherwise, we repeat once the same process between K q Y tpi 1 , j 1 qu and K q 1 Y tpi 1 , j 1 qu. In the end, we are left with different categories for different vortex sheet's contours. An example of what we obtain after Step 3 is displayed in Fig 3A .  Step 4: We build, from each set K p a list of grid points corresponding to an anticlockwise path along the countour of the p th vortex sheet. We proceed in the following way:

(1) For each list K p , we compute it's barycenter. We choose our starting couple from K p as one with a close x-axis coordinate to that of the barycenter and with the biggest y-axis coordinate. We add this to our new list K sort p . (2) While the last point and the first point of K sort p are not neighbours: (a) Out of all the neighbours in K p to the last grid point we added to K sort p , 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 to the last grid point we added do K sort p , we delete this grid point from K sort p 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. Remark 4. In the third point of step 4, we choose as criterion 70% arbitrarily. Other choices can be made as long as it is strictly bigger than 50%. If we choose a too high value, then the algorithm could be so constrained that it might not give a proper result.

Step 5: We compute the indices of each vortex sheet contour detected K sort p . The algorithm is similar to the last step in Section 4.1 with the list K sort p instead of S λpn,kq pn, kq.

Remark 5 (Computation of the index of a giant hole). For an estimation of the index of a giant hole, we first estimate its radius numerically. The algorithm goes as follows. We are given integers l 0 ě 1 and m 0 ąą 1 and a tolerance parameter tol 4 ą 0. We choose tol 4 " 0.1 in all the numerical experiments. For l P t0, ¨¨¨, l 0 u, we choose r l ą 0 a possible radius for the giant hole in such a way that l Þ Ñ r l is increasing, r 0 is too small a radius and r l 0 is too big a radius. For m P t0, ¨¨¨, m 0 u, we set θ m " 2πm m 0 . Let ψ be a minimizer with a giant hole. Starting from l " 0 and as long as l ď l 0 , we compute successively

β l " m 0 ´1 ď m"0 ˜˜argmin nPt1,¨¨¨,N u |x n ´rl cospθ m q| ¸ˆ˜a rgmin kPt1,¨¨¨,N u |y k ´rl cospθ m q| ¸¸. If max pn,kqPβ l |ψ n,k | 2 ą tol 4
, then we stop the algorithm and we will use r l and β 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 4.1 with β l instead of S λpn,kq pn, kq.

Numerical results

This section is devoted to the numerical experiments carried out using the EPG minimization method introduced in Section 3 and the post-processing algorithms described in Section 4 for the detection of structures and the computation of the indices. In Section 5.1, we consider a one component Bose-Einstein condensate in strong confinement and rotation, in connection with the theoretical results of Section 2.2.1. In Section 5.2, we consider a two components Bose-Einstein condensate in a strong confinement regime in the segregation case, without rotation. This allows to connect with the theoretical results of Section 2.2.2. In Section 5.3, we consider a two components Bose-Einstein condensate with strong confinement and rotation in the segregation case. This allows to connect with the theoretical results of Section 2.2.2. Last, in Section 5.4, we consider a two components Bose-Einstein condensate with strong confinement and rotation in the coexistence case. This allows for connection with theoretical conjectures of Section 2.2.3. 5.1. One component condensate with rotation. The goal of this section is to illustrate numerically the theoretical results described in Section 2.2.1. These results identify four different regimes for the behavior of the one component Bose-Einstein condensate depending on how big Ω ε is as ε tends to zero. These four regimes are separated by three characteristic rotational speeds Ω i ε , i " 1, 2, 3 (see ( 5)). We introduce below 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 ε in each of the four regimes.

5.1.1.

Parameters used for the one component simulations. For all the one component simulations, we consider the following parameters. The confinement is defined by the function [START_REF] Bao | Ground states of two-component bose-einstein condensates with an internal atomic josephson junction[END_REF] ρpx, yq " minr1, 10pR 2 ´x2 ´y2 qs.

For the initial datum we choose, ψ 0 px, yq " expp´10x 2 ´10y 2 q{5. The square of computation is of length 2L " 14 and the radius of the disc D is R " 4. We initialize the step size to h " 0.1 and we set h 0 " 10 ´12 . We also compute M " } ? ρ ą0 } 2 ∆ " 50.084.

Identification of the four regimes.

In order to identify the four regimes, we proceed as follows. First, we take ε " 10 ´1 and we use the EPG algorithm described in Section 3.4 to compute minimizers of E ∆ ε,δ for several rotational speeds. Then, we do the same for ε " 5 ˆ10 ´2. Based on these many simulations, and the expressions (5), we estimate the three critical values for the rotation speed which are shown in Figure 4. This provides us with an estimation of the three critical values of the rotational speed when ε " 10 ´2. Choosing rotation speed between these 3 estimated critical values for ε " 10 ´2, we manage to observe the four different regimes, as detailed below. ). We consider the case of no rotation (Ω " 0) and strong confinement (ε " 10 ´2). We set the number of points in the x-axis and y-axis to N " 512. The gradient descent algorithm converged due to the stopping criterion K ∆ (see [START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF]) which has a value of 10 ´2. The results are shown in Figure 5. As we can see, there are no vortices in the numerical minimizer. This is in accordance with the theory presented in section 2.2.1 (first case Ω ă Ω 1 ε ). Moreover the squared modulus of a minimizer converges to the non-negative part of the function ρ when ε tends to 0, as illustrated in Figure 6. ). We consider the case of a moderate rotation speed (Ω " 20) and strong confinement (ε " 10 ´2), according to Figure 4. First, we compute a minimizer with N " 512. Then, we interpolate its real and imaginary parts and compute another minimizer with N " 1024. The EPG method converged due to the stopping criterion K ∆ which has a value of 5 ˆ10 ´2. We compute numerically the indices of the vortices using the algorithm described in Section 4.1. The results are shown in Figure 7. As we can see, there are several vortices in the numerical minimizer and there is no sign of a giant hole in the center. The numerical index of most of the numerical vortices is equal to one. This validates numerically that most zeros of the wave function have a phase circulation. This is in accordance with the second case of the theory presented in Section 2.2.1.

Figure 7. The squared modulus of a minimizer for Ω " 20, ε " 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.

In another simulation, we consider the case of a moderate yet higher rotation speed (Ω " 30) and strong confinement (ε " 10 ´2). Once again, we use an interpolation from N " 512 to N " 1024. In this simulation, the EPG algorithm converged due to the stopping criterion K ∆ which has a value of 5 ˆ10 ´2. We also compute the indices of the vortices numerically using the algorithm described in Section 4.1. The results are shown in Figure 8. As we can see, there are vortices in the numerical minimizer and there is a dark disk in the center (not present in Figure 7) which indicates that we are close to the limit between zone 2 and zone 3 as predicted in Figure 4. The numerical index of most of the numerical vortices is equal to one. This validates numerically that most zeros of the wave function have a phase circulation. This is in accordance with the theory presented in Section 2.2.1.

Figure 8. The squared modulus of a minimizer for Ω " 30, ε " 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.

Big rotational speed (Ω 2

ε ă Ω ε ă Ω 3 ε ). We consider the case of a fairly high rotation speed (Ω " 40) and strong confinement (ε " 10 ´2) in accordance with the regimes obtained in Figure 4. This simulation is carried out with N " 1024. The EPG algorithm converged due to the stopping criterion K ∆ ď 5 ˆ10 ´2. The results are shown in Figure 9A. 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 rotation speed and is in accordance with the theory presented in section 2.2.1.

In Figures 9B and9C, we compute the indices of the vortices of the numerical minimizer displayed in Figure 9A. First, in Figure 9B, we compute the indices of the vortices in the annulus using the algorithm described in Section 4.1. Then, in Figure 9C, we compute the index of the giant hole (see Remark 5). The numerical parameters are indicated in the captions. As we can see, all the indices of the numerical vortices in Figure 9B around the giant hole are equal to one which validates numerically that almost all the zeros of the wave function have a phase circulation. In Figure 9C, the index of the giant hole is equal to 100. This is in accordance with the theory presented in section 2.2.1. 5.1.6. Huge rotational speed (Ω 3 ε ă Ω ε ). We consider the case of a huge rotation speed (Ω " 70) and strong confinement (ε " 10 ´2). Once again, we use an interpolation from N " 512 to N " 1024. The EPG algorithm converged due to the stopping criterion K ∆ ď 5 ˆ10 ´2. The results are shown in Figure 10A. 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 rotation speed and is in accordance with the theory presented in section 2.2.1 (fourth case Ω 3 ε ă Ω). In Figures 10B and10C, we compute the indices of the vortices of the numerical minimizer displayed in Figure 10A. First, in Figure 10B, we compute the indices of the vortices in the annulus using the algorithm of Section 4.1. Then, in Figure 10C, we compute the index of the giant hole (see Remark 5). As we can see, most of the indices of the numerical vortices in Figure 10B around the giant hole are equal to one, which validates numerically that most of the zeros of the function have a phase circulation. The index of the giant hole is equal to 522 (see Figure 10C). This is in accordance with the theory presented in section 2.2.1.

5.2.

Two components condensate without rotation (Ω " 0). In this section, we move on the theoretical results described in the first part of Section 2.2.2 that deal with two components Bose-Einstein condensates, without rotation, in different (strong, moderate and weak) segregation (δ ą 1) regimes.

Common numerical parameters.

In this section, we use the following parameters. The function ρ is the same as in [START_REF] Bao | Ground states of two-component bose-einstein condensates with an internal atomic josephson junction[END_REF]. The physical parameters are ε " 5 ˆ10 ´2, N 1 " 0.55 and N 2 " 0.45. The initial data are given by ψ 1 " ψ 2 " expp´10x 2 ´10y 2 q{5. The discretization parameters are L " 7, R " 4, N " 256, h " 0.1 and h 0 " 10 ´12 . The stopping criterion value for K ∆ is set to 10 ´2.

Note that for the simulations of Figures 11 and13), we started with the minimizers obtained in Figure 12 as initial data.

5.2.2.

The strong segregation regime: δ ε ˆε2 Ñ `8. We consider the case of strong segregation regime (δ ε " 4000) and strong confinement so that δ ε ˆε2 " 10. The results are shown in Figure 11. As we can see in Figure 11A and11B, the numerical support of the two components tend to not overlap, thereby confirming that we are in the segregation regime. In Figure 11C, we can see that the sum of the squared modulus of the two components has a minimum inside the disc (away from the boundary of the disc), with min px,yqPD |ψ 1˚|2 `|ψ 2˚|2 « 0.14. This is in accordance with the theory presented in section 2.2.2.

5.2.3.

The moderate segregation regime: δ ą 1 fixed and small ε. We consider the case of moderate segregation regime (δ " 1.5) and strong confinement, still without rotation (Ω " 0). The results are shown in Figure 12. The results displayed in Figures 12A and12B confirm that we are in the segregation regime (δ ą 1). Figure 12C presents the sum of the two squared moduli of the two components of the minimizer. We can see a white curve corresponding to the separation between the two components with min px,yqPD |ψ 1˚|2 `|ψ 2˚|2 « 0.78. This is in accordance with the theory presented in section 2.2.2 (second case with no rotation and fixed δ ą 1).

5.2.4.

The weak segregation regime: δ ε Ñ 1 with ε{ ? δ ε ´1 Ñ 0. We consider the case of segregation regime (δ ε " 1.02) and strong confinement, still without rotation (Ω " 0). The results are shown in Figure 13. Figures 13A and13B confirm that we are still in the segregation regime. In contrast to the two previous cases, we can see in figure 13C that the sum of the two squared moduli does not present a seperation area between the components. Indeed, apart from close to the boundary of the disc, the sum of the squared moduli is almost constant, with an approximate value of 0.976. With the notations of Section 2.2.2, we have ε " 0.35 in this third simulation. The results are in accordance with the theory (third case with no rotation, δ ε close to 1 and ε close to 0 in Section 2.2.2).

Two components condensate with rotation (Ω ‰ 0) in the segregation regime.

In this section, we consider physical parameters related to the theoretical results described in the second part of Section 2.2.2, that deal with two components Bose-Einstein condensates, in the segregation regime (δ ą 1), with rotation.

Common parameters.

In this section, we use the following parameters. The function ρ is the same as in [START_REF] Bao | Ground states of two-component bose-einstein condensates with an internal atomic josephson junction[END_REF]. The physical parameters are ε " 10 ´2, δ ε " 1 `ε " 1.01 (so that ε " ? ε " 0.1), N 1 " 0.55 and N 2 " 0.45. The initial data are given by ψ 1 " ψ 2 " expp´10x 2 ´10y 2 q{5. The discretization parameters are L " 7, R " 4, h " 0.1 and h 0 " 10 ´12 . All the numerical results in this Section 5.3 are obtained by first minimizing the discrete energy with N " 256, then interpolating the real and imaginary parts of each component to a grid of N " 512 points in each direction, then minimizing the corresponding discrete energy. For all the experiments in this section, the EPG algorithm of Section 3.2 converged because of the stopping criterion on K ∆ , with a value less or equal to 2.0 ˆ10 ´2.

These parameters correspond to the segregation regime (δ greater yet close to 1) and ε ăă ε. We consider the cases Ω P t1, 3, 6, 15u.

Low rotation case:

Ω " 1. We first consider the case of low rotation and strong confinement. The results are shown in Figure 14. The numerical experiment confirms that we are in a segregation regime (since δ ą 1) and the two components of the minimizer tend to not overlap. Moreover, the small rotation speed (Ω " 1) is not big enough to produce vortices in the minimizer. This is in accordance with the theory presented in section 2.2.2 (first bullet point in the rotational case). 5.3.3. Moderate rotation case: Ω " 3. We consider the case of moderate rotation speed and strong confinement. The results are shown in Figure 15. The numerical experiment shows that segregation holds, and the moderate rotation speed is big enough to produce singly quantized vortices in one of the components of the numerical minimizer (see Fig. 15C). This is in accordance with the theory presented in section 2.2.2 (second bullet point of the rotational case). 5.3.4. Average rotation case: Ω " 6. We consider the case of average rotation speed and strong confinement. The results are shown in Figures 16. Numerically, the experiment yields vortices in both components of the numerical minimizer, since the rotation is more important than in the two previous experiments. Moreover, the computation of the indices of the vortices in both components of the minimizer, using the algorithm described in Section 4.1, shows that the vortices are singly quantized (Fig. 16C and16D). This is in accordance with the theory presented in section 2.2.2 (second bullet point of the rotational case). 5.3.5. High rotation case: Ω " 15. We consider the case of high rotation velocity and strong confinement. The results are shown in Figures 17. The numerical experiment shows that we still are in a segregation regime, since the supports of the minimizers tend to not overlap. Moreover, the rotation speed is sufficiently big to observe numerically the formation of vortex sheets in each component of the minimizer. This is in accordance with the theory presented in section 2.2.2 (last bullet point of the rotational case). The computation of the indices of the vortex sheets is carried out using the algorithm described in Section 4.2, with m " 0.4, M " 0.6 and tol 3 " 0.3 in each component of the numerical minimizer. The computed indices of the vortex sheets are positive, which validates numerically the existence of vortex sheets with a phase circulation when δ ε Ñ 1 and ε Ñ 0, which was conjectured in [START_REF] Aftalion | Vortex patterns and sheets in segregated two component Bose-Einstein condensates[END_REF] (see also Section 2.2.2). 5.4. Two components condensate in the coexistence regime δ ă 1. In this section, we consider the numerical behavior of two components Bose-Einstein condensates in the coexistence regime (δ ă 1) as ε Ñ 0. As mentioned in 2.2.2, depending on the value of δ and Ω, we expect four different qualitative behaviours for the minimizers. 5.4.1. Common parameters. In this section, we use the following parameters. The function ρ is the same as in [START_REF] Bao | Ground states of two-component bose-einstein condensates with an internal atomic josephson junction[END_REF]. The physical parameters are ε " 5 ˆ10 ´2, N 1 " 0.55 and N 2 " 0.45. The discretization parameters are L " 7, R " 4, h " 0.1 and h 0 " 10 ´12 . All the numerical results in this Section 5.3 are obtained by first minimizing the discrete energy with N " 128, then interpolating the real and imaginary parts of each component to a grid of N " 256 points in each direction, then minimizing the corresponding discrete energy. For all the experiments in this section, the EPG algorithm of Section 3.2 converged because of the stopping criterion on K ∆ , with a value less or equal to 10 ´2. These parameters correspond to the coexistence regime (δ ă 1). Remark 6. When minimizing E ∆ ε,δ with high rotation speed Ω, giant holes are created in some of the next simulations. To circumvent this, from now on, we will refer to a discrete analogue of E Ω ε,δ (the energy with the centrifugal force) as E ∆ ε,δ . Of course, E ∆ ε,δ still denotes the discrete energy without the centrifugal force. See Remark 1.

5.4.2.

No vortices. We first consider the case of no rotation (Ω " 0), strong confinement and weak interaction strength δ " 0.9 (i.e. 1 ´δ " 0.1).

We start with a symmetric numerical initial datum corresponding to that of Section 5.2.1. The results are shown in Figure 18. The numerical experiment confirms that we are in a coexistence regime since the two components of the minimizer are disk-shaped, with almost constant moduli away from the boundary, with similar (hence overlapping) supports. Of course, the values of the constant moduli for ψ 1 and ψ 2 depend on the values of N 1 and N 2 . This is in accordance with the theory presented in Section 2.2.3. 5.4.3. Triangular vortex lattices. In this section, we consider the case of strong interaction strength (δ " 0.2, i.e. 1 ´δ " 0.8) and high rotation pΩ " 7q. In order to break off the ε,δ , in the case of no rotation (Ω " 0) and weak interaction strength (δ " 0.9). symmetry, we choose a decentered initial datum, given by ψ 1 px, yq " 1 5 expp´10px ´0.5q 2 ´10py `0.3q 2 q, ψ 2 px, yq " 1 5 expp´10px `0.7q 2 ´10py ´0.1q 2 q. (13)

The results are shown in Figure 19. The numerical minimizer is consistent with the coexistence regime, and the high rotation is big enough to produce singly quantized vortices in both components forming a triangular vortex lattice. This is in accordance with the second behaviour in the theory presented in Section 2.2.3. (Ω " 7) and strong interaction strength (δ " 0.2) for the energy E ∆ ε,δ followed by red lines highlighting the triangular lattices. 5.4.4. Square vortex lattices. In this section, we consider the case of weak interaction strength (δ " 0.8, i.e. 1 ´δ " 0.2) and high rotation (Ω " 8). We use the same non-symmetric initial data as before (see Equation 13). The results are shown in Figure 20. Numerically, the minimizer confirms that we are in a coexistence regime, and the weak interaction strength with the high rotation produce a square vortex lattice with singly quantized vortices in both components of the minimizer. This is in accordance with the third behavior in the theory presented in Section 2.2.3. In a first simulation, we use [START_REF] Bao | Efficient numerical methods for computing ground states and dynamics of dipolar bose-einstein condensates[END_REF] as an initial datum as before. We consider the case of fairly high rotation speed (Ω " 4) and weak interaction strength (δ " 0.998). The numerical results obtained are displayed in Figure 21. They confirm that we are in a coexistence regime where each component is disk-shaped. Indeed, in both components, the points inside the disk of radius R where the modulus of the minimizer is very small are isolated, and there is at least a bit of mass in each component in all the rest of the disk. This contrasts for example with all the numerical experiments carried out in the segregation regime (see Section 5.3).

We compute the index of the detected vortices, and we display the results in Figure 21 as well. As expected, almost all the vortices are paired up 2 by 2, forming double-core vortices. Moreover, almost all the indices of the numerical vortices are equal to one, which validates numerically that the zeros of the function have a non-trivial and singly quantized phase circulation. The only exception is one vortex in the first component (see Figure 21A) which has an index of 2. An explanation for this is that it consists in two vortices that are too close to each other, and hence cannot be detected as two different vortices by our vortex detection algorithm (see Section 4.1). This is in accordance with the last behaviour predicted by the theory presented in Section 5. [START_REF] Anderson | Observation of Bose-Einstein condensation in a dilute atomic vapor[END_REF].

In a second simulation, we use the following initial datum ψ 1 n,k " ψ 2 n,k " sinpx n,k `yn,k q, still with the convention that ψ 1 and ψ 2 satisfy the Dirichlet boundary conditions. We consider The numerical results confirm that we are in the coexistence regime, where each component is disk-shaped. In Figure 22, we also compute and display the index of the vortices detected. As expected, almost all the indices of the numerical vortices are equal to 1, which validates numerically that the zeros of the function have a non-trivial and singly-quantized phase circulation. An explanation for the two exceptions (doubly-quantized vortices in the second component displayed in Figure 22B) is once again that our algorithm simply fails to detect two separate vortices. Moreover, we can see a stripe pattern in both components (Figures 22A and22B). This is in accordance with the last behaviour in the theory presented in 5.4.

Comparison with GPELab

In this section, we study the efficiency of the EPG method introduced in Section 3 compared to that of GPELab's method (see [START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations I: computation of stationary solutions[END_REF][START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations II: Dynamics and stochastic simulations[END_REF]) for minimizing the energy E ∆ ε,δ . At first sight, the EPG method is explicit whereas GPELab's method solves a linear system at each time step. The time step is denoted by ∆ t in GPELab, while we denote it by h in this paper. We display for each comparison test three convergence criteria (see below) as well as the number of iterations, the energy and the execution time. Let ψ n denote the n th iteration using either EPG or GPELab's algorithm. We consider and compute the following three convergence criteria:

(1) The difference between two successive iterations }ψ n ´ψn´1 } 2 8 , (2) The criterion K ∆ that we developed (see Section 3.3), (3) The energy evolution |E ∆ ε,δ pψ n q ´E∆ ε,δ pψ n´1 q|. Even though we compute all these quantities in all the numerical experiments in this Section, we use as a stopping criterion the first criterion above (which is native in the code) for GPELab's method, and the first two criteria for the EPG method (in contrast to the previous numerical simulations of Section 5 where only the criterion on K ∆ is used). For each test case, we compute the final value of the three criteria and we compare them.

Remark 7. The energy computed by GPELab's method is different from the one given in (4). The continuous GPELab's energy reads as follows (see [START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations I: computation of stationary solutions[END_REF][START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations II: Dynamics and stochastic simulations[END_REF] for more details):

E GP E pu 1 , u 2 q " ż D 2 ÿ "1 " 1 2 |∇u | 2 `V |u | 2 ´Ωu ˚ L z u  `żD u ˚β 2
Fu, where V " ´ρ 2ε 2 is the confinement function associated to component number , β " 1 2ε 2 is the intra-component interaction, u ˚is defined as pĎ u 1 , Ď u 2 q, the matrix F is given by

F " ˆβ1,1 |u 1 | 2 `β1,2 |u 2 | 2 0 0 β 2,1 |u 1 | 2 `β2,2 |u 2 | 2 ˙, where ˆβ1,1 β 1,2 β 2,1 β 2,2
˙" ˆ1 δ δ 1 ˙, and the rotation operator is defined by L z " ´i`x B x ´yB x ˘.

The relation between the discrete analogue to E GP E and E ∆ ε,δ (see [START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates[END_REF]) is the following:

E ∆ GP E " E ∆ ε,δ ´δx δ y 4ε 2 N ÿ n"0 K ÿ k"0 `ρ2 n,k ˘ρą0 .
Observe that the difference between the two energies does not depend on u " pu 1 , u 2 q.

In the case of one component condensate (by taking u 2 " 0), we take β 1,1 " 1 and β 1,2 " β 2,1 " β 2,2 " 0, so that the last term of

E GP E reads ş D u ˚β 2 Fu " ş D β 2 |u 1 | 4 . 6.1. One component comparison.
In this Section, we consider a one component Bose-Einstein condensate. We compare both algorithms (GPELab's method and EPG) on the same computer. Each test is done separately. 6.1.1. Parameters for the EPG method. We take the following parameters for all the tests. The function ρ is the same as before (see [START_REF] Bao | Ground states of two-component bose-einstein condensates with an internal atomic josephson junction[END_REF]). We take ε " 5 ˆ10 ´2. For the initial datum, we choose ψ 1 " 1 5 expp´10x 2 ´10y 2 q. For the discretization parameters, we take the initial step h " 0.256, the tolerance h 0 " 10 ´12 , the length of the box 2L " 14, the radius of the disk R " 4, the species ratios N 1 " 1 and N 2 " 0, the number of points per direction N " 256 and the convergence tolerance K 0 " 10 ´2.

6.1.2. Parameters for GPELab's method. The equivalent of these parameters in GPELab (with the notations of [START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations I: computation of stationary solutions[END_REF][START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations II: Dynamics and stochastic simulations[END_REF]) are ' Ncomponents" 1, ' Type='BESP' (for Backward Euler pseudoSPectral scheme), ' Delta" 0.5 (the coefficient in front of the kinetic energy), ' Beta" 200, which corresponds to 1 2ε 2 for ε " 0.05, ' V px, yq " ´1 2ε 2 minr1, 10pR 2 ´x2 ´y2 qs (we had to add our own function), ' Initial datum: ψ 1 " 1 5 expp´10x 2 ´10y 2 q (we had to add it too), ' x min " y min " ´7 and x max " y max " 7. We had to modify the normalization step after each time step ∆ t . 6.1.3. Numerical results. The GPELab algorithm converges whenever the infinity norm of the difference between two successive iterations is less than a certain value G 0 " ∆ t ˆStop crit. It stops without convergence if the number of iterations done exceeds 10 6 . In this Section, the EPG method converges if either one of the following statements is true: at an iteration step n, }ψ 1 n ´ψ1

n´1 } 8 ă G 0 or K ∆ ă K 0 . For a first comparison, we set the rotation speed to Ω " 1. We also set ∆ t " 2 ˆ10 ´3 and Stop crit" 10 ´2 so that G 0 " 2 ˆ10 ´5. We display the results in Figure 23. For a second comparison, we set the rotation speed to Ω " 3. We also set ∆ t " 2 ˆ10 ´3 and Stop crit" 10 ´2 so that G 0 " 2 ˆ10 ´5. We display the results in Figure 24. For a third and final comparison, we set the rotation speed to Ω " 6. We also set ∆ t " 10 ´3 and Stop crit" 10 ´2 so that G 0 " 10 ´5. The results are displayed in Figure 25.

The numerical results displayed in Figures 23,24 and 25 show that the methods may converge to different minimizers, with different energies (see Figure 23), as well as they may converge to similar minimizers with similar energies (see Figure 25), taking into account suitable rescaling of the discrete energies (as indicated in Remark 7). Moreover, the qualitative behaviour of both methods is correct in the sense that higher rotations end up in creating more vortices, and the beginning of a central hole when Ω " 6 (see Figure 25). Note that the EPG method uses more interations in all cases. However, since the EPG method is explicit and the GPELab method is linearly implicit, the EPG method is faster in all cases. The order of magnitude of the speed up is 3. Finally, the stopping criterion K ∆ is much smaller in minimizers obtained by EPG than in that computed using GPELab. 6.2. Two components comparison. In this section, we consider two components Bose-Einstein condensates. We also compare both algorithms on the same computer. Each test is done separately. 6.2.1. Parameters for the EPG method. We take the following parameters for all the tests. We take ε " 5 ˆ10 ´2. For the initial datum, we choose ψ 1 " ψ 2 " 1 5 expp´10x 2 ´10y 2 q. For the discretization parameters we take the values N 1 " 0.55 and N 2 " 0.45, N " 256 and K 0 " 10 ´2. For the rest of the parameters, we refer to 6.1.1. (A) A minimizer of the energy E ∆ ε using the EPG method.

(B) A minimizer of the energy E GP E using the GPELab method. ' Initial datum: ψ 1 " ψ 2 " 1 5 expp´10x 2 ´10y 2 q, ' ∆ t " 5 ˆ10 ´4 and Stop crit" 2 ˆ10 ´2 (G 0 " 10 ´5), ' x min " y min " ´7 and x max " y max " 7. We also had to modify the normalization step after each time step ∆ t . 6.2.3. Numerical results. For a first comparison, we set the rotation speed to Ω " 3 and the interaction strength to δ " 0.7 (coexistence). We display the results in Figure 26. For a second comparison, we set the rotation speed to Ω " 3 and the interaction strength to δ " 1.5 (segregation). For the EPG method, to break off the symmetry, we chose the following initial datum ψ 1 " 1 5 expp´10px ´0.5q 2 ´10py `0.2q 2 q, ψ 2 " 1 5 expp´10px `0.5q 2 ´10py ´0.3q 2 q.

We display the results in Figure 27.

The results displayed in Figures 26 and27 both indicate that the methods may converge to different minimizers, even if their energies are similar (see Remark 7 for details on the rescaling). The qualitative behaviour of both methods is correct in the sense that they provide coexisting minimizers when δ " 0.7 and segregated minimizers when δ " 1.5. The number of iterations in both methods has the same order of magnitude (the maximal ratio is of order 2), in contrast to the one component simulations carried out in Section 6.1. The EPG method is much faster in both cases: The order of magnitude of the speed up is 3 when δ " 0.7 and 25 when δ " 1.5. Finally, the stopping criterion K ∆ is much smaller in minimizers obtained by EPG than in that computed using GPELab.

Conclusion

In this paper, we considered a Gross-Pitaevskii energy as a model for rotating one component and two components Bose-Einstein condensates in two dimensions in a strong confinement and strong rotation regime. First, we have introduced a discretization of this energy, using the FFT scheme and Plancherel's equality, which allowed a reasonable computation of the energy gradients. Second, in contrast to the literature, we proposed a minimization method for this discrete energy using an explicit L 2 gradient method with projection (EPG). This method allowed for the derivation of a stopping criterion. Third, we introduced two post processing algorithms for the numerical minimizers. One is aimed for the single vortices and the other for vortex sheets. Both allow to detect these structures and compute their indices. Fourth, we have ran our methods and algorithms for different physical regimes from one component condensates with high rotation to two components condensates in coexistence and segregation regimes. This validates recent theoretical results as well as it supports conjectures (as for example the existence of vortex sheets in the segregation regime). Last, we have compared the efficiency of the EPG method to that of the GPELab method [START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations II: Dynamics and stochastic simulations[END_REF]. On all the numerical tests presented, EPG appears to be, roughly speaking, at least 3 times faster than GPELab.

Figure 1 .

 1 Figure 1. S λ pn, kq with angles θ 0 to θ 8λ .

  (A) An example of the contours detected after the first part of the first step of the contour detection algorithm. (B) An example of the connex curves obtained at the end of the first step of the contour detection algorithm.

Figure 2 .

 2 Figure 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). We display just a part of the minimizer obtained in Fig 17A.

  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 Fig 2B. Our next step is to separate the connex components in K, to categorize the different sheets (since each minimizer can have more than one vortex sheet).

  (A) After Step 3: We can use a different colour of each contour (the colour changes with p). (B) After Step 4: We can use color shading to show the anticlockwise orientation of each contour.

Figure 3 .

 3 Figure 3. An example of the contours detected and orientated by the contour detection algorithm after the third and fourth steps.

( 3 )

 3 If the length of K sort p is less than 70% of the original size of K p , then we delete the last grid point added from K sort p and from K p . Next, we return to the (new) last grid point we added to K sort p before that and we repeat the second point from Step 4. If the length of K sort p exceeds 70% of the original size of K p , then we stop the algorithm. An example of what we obtain after Step 4 is displayed in Fig 3B.

Figure 4 .

 4 Figure 4. Ω as a function of ε and the associated zones.

5. 1 . 3 .

 13 Small rotational speed (Ω ε ă Ω 1 ε

  (A) The graph of the non-negative part of ρ.(B) The squared modulus of a minimizer.

Figure 5 .

 5 Figure 5. Comparison between the non-negative part of ρ (figure (A)) and the squared modulus of a minimizer for the energy E ∆ ε (figure (B)) with no rotation (Ω " 0) and strong confinement (ε " 10 ´2) computed with N " 512 and K ∆ " 10 ´2.

Figure 6 .

 6 Figure 6. Comparison of ρ and the squared modulus of a minimizer when ε " 0.1 and ε " 0.01 along the y-axis.

Figure 6

 6 Figure 6 numerically illustrates the convergence of the squared modulus of a minimizer for the energy E ∆ ε to the non-negative part or ρ when ε tends to 0 at fixed rotation speed Ω " 0. 5.1.4. Moderate rotational speed (Ω 1 ε ă Ω ε ă Ω 2 ε). We consider the case of a moderate rotation speed (Ω " 20) and strong confinement (ε " 10 ´2), according to Figure4. First, we compute a minimizer with N " 512. Then, we interpolate its real and imaginary parts and compute

( A )

 A The squared modulus of a minimizer. (B) The vortices' indices of a minimizer. (C) The index of the giant hole of a minimizer.

Figure 9 .

 9 Figure 9. The squared modulus (A) of a minimizer of the energy for Ω " 40 and ε " 10 ´2. The vortices' indices of a minimizer (B) detected with N min " 1, N max " 3, tol 1 " 0.05 and tol 2 " 0.02. The giant hole's index (C) with r " 1.523 (see remark 5).

( A )

 A The squared modulus of a minimizer. (B) The vortices' indices of a minimizer. (C) The index of the giant hole.

Figure 10 .

 10 Figure 10. The squared modulus (A) of a minimizer of the energy for Ω " 70, ε " 10 ´2. The vortices' indices of a minimizer of the energy (B) detected with N min " 1, N max " 3, tol 1 " 0.05 and tol 2 " 0.02. The index of the giant hole (C) detected with r " 2.85 (see Remark 5).

( A )

 A The squared modulus of the first component. (B) The squared modulus of the second component. (C) The sum of the squared moduli of both components.

Figure 11 .

 11 Figure 11. The squared modulus of the first component (A) and the second component (B) for a minimizer in the case of two components condensate with no rotation and δ ε " 4000. Fig (C) displays their sum.

( A )

 A The squared modulus of the first component. (B) The squared modulus of the second component. (C) The sum of the squared moduli of both components.

Figure 12 .

 12 Figure 12. The squared modulus of the first component (A) and of the second component (B) for a minimizer in the case of two components condensate with no rotation and δ " 1.5. Their sum is displayed in Fig. (C).

( A )

 A The squared modulus of the first component. (B) The squared modulus of the second component. (C) The sum of the squared moduli of both components.

Figure 13 .

 13 Figure 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 no rotation and δ ε " 1.02. Their sum is displayed in (C).

( A )

 A The squared modulus of the first component.(B) The squared modulus of the second component.

Figure 14 .

 14 Figure 14. 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 Ω " 1 in the segregation regime.

  (A) The squared modulus of the first component. (B) The squared modulus of the second component. (C) The vortices' indices of the second component of a minimizer for the energy detected with N min " 1, N max " 5, tol 1 " 0.05 and tol 2 " 0.02.

Figure 15 .

 15 Figure 15. The squared modulus of the first component (A) and the second component (B) of a minimizer in the case of two components condensate with moderate rotation Ω " 3. The vortices' indices of the second component are presented in Fig. (C).

( A )

 A The squared modulus of the first component. (B) The squared modulus of the second component. (C) The vortices' indices of the first component. (D) The vortices' indices of the second component.

Figure 16 .

 16 Figure 16. The squared modulus of the first component (A) and the second component (B) of a minimizer in the case of two components condensate with moderate rotation Ω " 6. The vortices' indices of the first and second component are presented in Fig. (C) and (D) (respectively) with N min " 1, N max " 5, tol 1 " 0.05 and tol 2 " 0.02.

( A )

 A The squared modulus of the first component. (B) The squared modulus of the second component. (C) The vortex sheet's indices of the first component. (D) The vortex sheets' indices of the second component.

Figure 17 .

 17 Figure 17. The squared modulus of the first component (A) and the second component (B) of a minimizer in the case of two components condensate with high rotation Ω " 15. The vortex sheets' indices of the first and second component are presented in Fig. (C) and (D) (respectively) with m " 0.4, M " 0.6 and tol 3 " 0.3.

( A )

 A The squared modulus of the first component.(B) The squared modulus of the second component.

Figure 18 .

 18 Figure 18. The squared modulus of the first component (A) and of the second component (B) of a minimizer of the energy E ∆ε,δ , in the case of no rotation (Ω " 0) and weak interaction strength (δ " 0.9).

( A )

 A The vortices indexes of the first component.(B) The vortices indexes of the second component.

Figure 19 .

 19 Figure 19. The vortices' indices of the first component (A) and second component (B) in the case of a two components condensate with high rotation (Ω " 7) and strong interaction strength (δ " 0.2) for the energy E ∆ ε,δ followed by red lines highlighting the triangular lattices.

  (A) The vortices' indices of the first component. (B) The vortices' indices of the second component.

Figure 20 .

 20 Figure 20. The vortices' indices of the first component (A) and second component (B) in the case of a two components condensate with high rotation (Ω " 8) and weak interaction strength (δ " 0.8) for the energy E ∆ ε,δ followed by red lines highlighting the square lattice.

( A )

 A The squared modulus and the vortices' indices of the first component for N min " 3, N max " 4. (B) The squared modulus and the vortices' indices of the second component for N min " 1, N max " 4.

Figure 21 .

 21 Figure 21. The squared modulus of the first component (A) and the second component (B) of a minimizer (in the case of double-core pattern) for the energy E ∆ ε,δ with Ω " 4 and δ " 0.995. The vortices' indices of each component are also presented in (A) and (B) (respectively) for tol 1 " 0.05 and tol 2 " 0.01.

  (A) The squared modulus and the vortices' indices of the first component for N min " 1, N max " 3. (B) The squared modulus and the vortices' indices of the second component for N min " 1, N max " 4.

Figure 22 .

 22 Figure 22. The squared modulus of the first component (A) and the second component (B) of a minimizer (in the case of stripe pattern) for the energy E ∆ ε,δ with Ω " 5 and δ " 0.98. The vortices' indices of each component are also presented in (A) and (B) (respectively) for tol 1 " 0.05 and tol 2 " 0.01.

(Figure 23 .

 23 Figure 23. Comparison between GPELab and EPG results for Ω " 1.

Figure 24 . 6 . 2 . 2 .Figure 25 .

 2462225 Figure 24. Comparison between GPELab and EPG results for Ω " 3.

( A )Figure 26 .

 A26 Figure 26. Comparison between GPELab and EPG results for Ω " 3 and δ " 0.7.

(A)Figure 27 .

 27 Figure 27. Comparison between GPELab and EPG results for Ω " 3 and δ " 1.5.

  and Λ n,k " λ k ,

	with
	ξ k "	´πpN `1q 2L	`kδ ξ , and λ n "	´πpN `1q 2L	`nδ λ ,

  1, we consider the values of |ψ| 2 on the squares S λ pn, kq "

	" pn ˘λ, k `jq ˇˇj P t´λ, . . . , λu	*	ď " pn `j, k ˘λq ˇˇj P t´λ, . . . , λu	* ,
	of length 2λδ

x , for λ P tN min , . . . , N max u. If for one of these λ " λpn, kq, the values of |ψ| 2 at all points of the square S λ pn, kq are such that |ψ| 2 ´|ψ n,k | 2 ą tol 2 , then we add the center pn, kq to the second list P, and we set λpn, kq as the caracteristic length of the potential vortex.

Acknowledgments

This work was partially supported by the Labex CEMPI (ANR-11-LABX-0007-01).