Bose-Einstein Condensation on Curved Manifolds

Here we describe a weakly interacting Bose gas on a curved manifold, which is embedded in the three-dimensional Euclidean space.~To this end we start by considering a harmonic trap in the normal direction of the manifold, which confines the three-dimensional Bose gas in the vicinity of its surface.~Following the notion of dimensional reduction as outlined in [L.~Salasnich et al., Phys.~Rev.~A {\bf 65}, 043614 (2002)], we assume a large enough trap frequency so that the normal degree of freedom of the condensate wave function can be approximately integrated out. In this way we obtain an effective condensate wave function on the quasi-two-dimensional surface of the curved manifold, where the thickness of the cloud is determined self-consistently. For the particular case when the manifold is a sphere, our equilibrium results show how the chemical potential and the thickness of the cloud increase with the interaction strength.~Furthermore, we determine within a linear stability analysis the low-lying collective excitations together with their eigenfrequencies, which turn out to reveal an instability for attractive interactions.


Introduction
Bose-Einstein condensates (BECs) are today among the most explored many-body quantum systems. Ranging from more simple thermodynamics [1,2] to the more complex turbulent regime [3,4,5], BECs represent an interesting workbench for a large variety of both experimental and theoretical phenomena. In particular, there are many tuning knobs to change the respective system properties in a highly controllable way, among them most importantly its dimensionality. Three-dimensional (3D) systems allow to investigate many fundamental macroscopic quantum phenomena as, for instance, superfluidity [1,2,6]. The two-dimensional (2D) case enables both the formation and the dynamics of vortices or the realization of transitions such as the one of the Berezinskii-Kosterlitz-Thouless (BKT) type [7,8,9], which is characterized by the binding or unbinding of vortex-antivortex pairs. But so far all experimental confinements of 2D Bose gases took place only in flat space. However it was predicted that interesting topological effects on 2D curved manifolds may occur due to the presence of vortices and their dynamics [10,11,12].
Some time ago it was proposed to use the coupling between the internal atomic structure and external electrodynamic fields in order to produce a spherical shell [13,14]. With this it is possible to imagine the realization of a BEC on a curved manifold. Since then there have already been several attempts to confine a Bose gas on the surface of a sphere, which leads to the so-called bubble trap [15,16,17]. However, due to the gravitational sag, the gas tends to concentrate on the bottom of the sphere. How to realize a Bose gas to occupy the whole sphere is currently an experimental challenge. One way to avoid the effects of gravity is to perform experiments within microgravity settings [18,19,20]. The recent installation of the NASA Cold Atom Laboratory (CAL) at the International Space Station (ISS) nourishes the prospect to realize soon such a bubble trap in microgravity, so that a BEC on a 2D curved manifold will become an experimental reality [21,22,23]. The planned experimental upgrade BECCAL at the ISS will even allow for a binary Bose mixture of rubidium and potassium to be confined in such a bubble trap [24].
Inspired by this notion of a bubble trap, several theoretical predictions have already been obtained for both static and dynamic properties of a Bose gas confined on a spherical geometry. For instance, the collective modes of such a system were already investigated for different regimes [25]. Furthermore, a crossover between 3D and 2D was considered, where the respective limits correspond to a completely filled and a hollow sphere [26,27]. Additional studies on a sphere deal with the critical temperature for the onset of Bose-Einstein condensation [28] and superfluidity [29]. In both cases the limit of an infinitely large radius of the sphere reproduces the corresponding 2D Euclidean results, i.e., in the former case the Mermin-Wagner-Hohenberg theorem [30,31] and in the latter case the BKT phase transition [7,8]. Moreover, the impact of quantum fluctuations and the formation of clusters were studied [32]. Quite recently, even the ground state and collective excitations were analysed for a dipolar Bose-Einstein condensate in a bubble trap [33]. Furthermore, the free expansion of a hollow condensate was also studied [34].
Obviously a sphere does not represent the only possible convenient geometry for confining a Bose gas, therefore a more general consideration of a curved manifold is necessary. For instance, one can use as a starting point a generic smooth manifold and proceed then with the modelling by taking into account the influence of asymmetries and deformations on the static and thermodynamic properties. Following this notion, we restrict ourselves in the present paper to the case of such a smooth manifold. To this end we start in section 2 by defining the relevant basic mathematical objects such as the Gaussian normal coordinate system to describe a manifold and the Laplace-Beltrami operator, which extends the Laplacian to a curved manifold. Afterwards, we follow the notion of reference [35] to implement a dimensional reduction. Starting from a 3D mean-field description, we introduce in section 3 a confining potential, which restricts the condensate in the direction perpendicular to the manifold. Correspondingly, the 3D wave function factorizes into a 2D part, representing the condensate on the manifold, and a one-dimensional Gaussian function, which describes the confinement. With this we derive in section 4 from a 3D Gross-Pitaevskii equation a self-consistent set of equations for both the 2D condensate wave function on the manifold and its width. Then, section 5 discusses the equilibrium case for a sphere. In section 6 we formulate the dynamical equations on a generic smooth manifold, while in section 7 we specialize this dynamics to a sphere and perform subsequently a linear stability analysis. With this we calculate the low-lying frequencies, which turn out to be stable for repulsive interactions. Finally, section 8 analyses the modes corresponding to these frequencies. It turns out that oscillations with a higher frequency predominantly occur in the direction of the confinement, whereas oscillations with lower frequencies are mainly restricted on the sphere.

Differential Geometrical Preliminaries
Due to the fact that in this work we pursue a more general approach to the problem of Bose-Einstein condensation in special geometries, it is necessary to introduce relevant aspects of differential geometry that are used later on. To this end we consider M to be the manifold [36] to which the Bose gas is confined. The way to describe a manifold is not uniquely defined, i.e., there are a lot of coordinate systems that could be used. Furthermore, it is often not possible to describe the whole manifold with only one coordinate system. In those cases we would have a local coordinate system for each portion of the manifold and we should put all these pieces together to have a global description. The way of choosing these portions and their number depend on the respective manifold and on the chosen coordinate systems. Anyhow, we can assure that this number can be chosen to be finite for compact manifolds, such as spheres and ellipsoids. For not compact manifolds it could be that this number of portions is infinite, but it can always be chosen to be countable. For simplicity, we suppose in the following that the manifold is described by only one coordinate system. Even though it is not true for many examples, it will not affect our main result. Furthermore, all the arguments presented here could be straightforwardly reproduced for the case of needing to describe the manifold with more than one piece.
Here we choose to use a Gaussian normal coordinate system [37,38]. It is always possible to do that: any manifold can be locally described by a Gaussian normal coordinate system. According to the illustration in figure 1, these coordinates are defined as follows. Consider a small portion of M, where all the points can be characterized by two real variables (x 1 , x 2 ) belonging to an open subset of the real plane IR 2 . So we can denote the points on this manifold by the vector p(x 1 , x 2 ). Let v 0 (x 1 , x 2 ) be a unit normal vector to the manifold at the point p(x 1 , x 2 ). Now, we describe the neighbouring points, which do not belong necessarily to M, by Note that the points belonging to M are described by x 0 = 0 and that ones, which do not belong to M, are described by x 0 = 0. Furthermore, fixing any constant value for x 0 , a locally parallel manifold to M is defined, which we denote by M(x 0 ). In particular, M(0) simply coincides with the manifold M. One detail to be pointed out is that equation (1) only represents a local description, so, in principle, x 0 cannot be arbitrarily large. We consider that the Gaussian normal coordinate system is well-defined in the interval |x 0 | < R/2, where R is the minimum, over all points p belonging to M, of the smallest curvature radius of each point p. In order to be more precise, for any fixed point p we denote by R 1 (p) and R 2 (p) the two principal curvature radii of the manifold at that point. Choose as R(p) the minimum between these two radii, i.e., Then, define R as the minimum value of R(p) over the whole manifold M, i.e., For the special case of a sphere, its radius coincides with R. Due to global properties, it could be for some manifold that the Gaussian normal coordinate system describes twice the same point. We exclude such manifolds in the following, restricting ourselves only to manifolds where this situation does not occur. Now we present a heuristic recipe for deriving equation (1). To this end, we suppose to know the manifold equation, which determines the points p belonging to a manifold portion as a function of x 1 and x 2 , that is p = p(x 1 , x 2 ). In order to find the respective tangent vectors we evaluate Then the unit normal vector to the manifold at p(x 1 , x 2 ) is defined by the cross product between these tangent vectors and a subsequent normalization, i.e., With this, equation (1) is well-defined and |x 0 | = |(q − p) · v 0 | defines the distance of the point q from the manifold.  (1) by setting x 0 = 0. The vectors v 1 and v 2 are two tangent vectors and v 0 (x 1 , x 2 ) is a unit normal vector to the manifold at the point p(x 1 , x 2 ). The neighbouring point q(x 0 , x 1 , x 2 ) does not necessarily belong to the manifold M, and the Gaussian normal coordinate system defines this point via equation (1) by identifying x 0 with the distance between q and the manifold. In the illustration, p represents the point belonging to M, which is the nearest to q. Now, we can define the tangent vectors to the surface M(x 0 ), which is locally parallel to M, using an analogous procedure. From formula (1) we calculate the corresponding tangent vectors as v i (x 0 , x 1 , x 2 ) = ∂q(x 0 , x 1 , x 2 )/∂x i , for i = 1, 2. The normal vector to the manifold M(x 0 ) at the point q( and turns out to be v 0 (x 0 , x 1 , x 2 ) = v 0 (x 1 , x 2 ). The latter statement can be seen from the property v 0 (x 1 , x 2 ) · v i (x 0 , x 1 , x 2 ) = 0, which follows due to for i = 1, 2. With that, a basis for the 3D space at each point q(x 0 , x 1 , x 2 ) is given by . This allows us to define a covariant metric of the 3D space in the neighbourhood of the manifold M, since each component of the metric is defined by the scalar product between two respective basis vectors, i.e., where µ and ν range from 0 to 2. From the definition of the normal vector v 0 in equation (5) and from the above discussion, we conclude both G 00 (x 0 , x 1 , x 2 ) = 0, for i = 1, 2, for all x 0 , x 1 and x 2 where the coordinate system is well-defined. With that and using the Gaussian normal coordinate system, we obtain that the covariant metric for the 3D space in the surrounding of the manifold M can be represented by the matrix for µ and ν ranging from 0 to 2, while i and j range only from 1 to 2. For x 0 fixed, g ij (x 0 , x 1 , x 2 ) denotes the covariant metric of the manifold M(x 0 ). Each entry of this metric represents the scalar product of the respective tangent vectors v 1 (x 0 , x 1 , x 2 ) and v 2 (x 0 , x 1 , x 2 ). For the special case x 0 = 0, g ij (0, x 1 , x 2 ) represents the metric of the manifold M and we use the abbreviated notation g ij (x 1 , x 2 ). Note that the metric g ij (x 0 , x 1 , x 2 ) can be Taylor expanded around the metric g ij (x 1 , x 2 ) and written in terms of other local properties of the manifold M, as is summarized in appendix A. There, it is also shown that calculating the square root of the determinant of this metric yields for |x 0 | < R/2, where we have introduced the notation det g = det g ij (x 1 , x 2 ). This result turns out to be useful in the next sections. We know that in many-body quantum theory the kinetic energy of a 3D BEC is given by the Laplacian of the condensate wave function. A generalization of the Laplacian for the 3D space in a generalized coordinate system is called Laplace-Beltrami operator [39]. It is expressed by with η and κ ranging from 0 to 2 and det G denoting the determinant of the covariant metric. Note that we use the Einstein notation in (11), so that a summation over the same co-and contravariant indices is implicitly assumed. For a flat 3D space described by a Cartesian coordinate system, the metric is given by an Euclidean metric with the Kronecker symbol as its respective components and the above formula recovers the standard representation of the Laplacian. On the other side, still for the flat 3D space represented via the Gaussian normal coordinate system, using the metric (9), the Laplace-Beltrami operator (11) reduces to with the abbreviation Note that (13) denotes also a Laplace-Beltrami operator, but this time in the context of the manifold M(x 0 ) for each fixed value of x 0 , and the indices i, j range from 1 to 2.
In particular, when x 0 = 0, the operator (13) represents the Laplace-Beltrami operator of the manifold M, denoted in the following simply by ∆ M .

Normalization of Condensate Wave Function
Now that we have our mathematical objects well-defined, we introduce two important physical quantities, which we need to consider when dealing with a BEC on a manifold: the potential which confines the Bose gas to the manifold and a particular ansatz for the condensate wave function. We suppose that the Bose gas is confined in the immediate vicinity of the manifold. Such a confinement could be realized, for instance, by a harmonic oscillator potential in the normal direction to the manifold, which has its minimum on the manifold. In the Gaussian normal coordinate system introduced in the previous section, this potential has the form Here the particle mass M and the frequency ω define a length scale in terms of the oscillator length σ osc = h/M ω, which represents the order of magnitude of the thickness of the Bose gas cloud surrounding the manifold. In the following we assume that the frequency ω is so large that the oscillator length σ osc is much smaller than the minimum R of the respective local curvature radii on the manifold, i.e., σ osc R. This corresponds to the physical situation that the Bose gas forms a thin shell around the manifold.
In addition to this harmonic potential, we also allow the Bose gas to be affected by a potential U depending on the manifold coordinates x 1 and x 2 . Thus, the total potential is given by Note that V harm (x 0 ) and U (x 1 , x 2 ) are 3D potentials, even though each one does not depend on all three variables. Let us now compute as a physical quantity the number of particles. Denoting the condensate wave function Ψ(x 0 , x 1 , x 2 ), the number of particles is given by an integral over the whole 3D space of its squared norm, i.e., N = dV |Ψ(x 0 , x 1 , x 2 )| 2 . But since the gas is confined in the vicinity of the manifold M, the integral becomes naturally restricted to the manifold neighbourhood We consider that this neighbourhood is defined by the points q(x 0 , x 1 , x 2 ) of the 3D space described by the Gaussian coordinate system according to equation (1) with |x 0 | < R/2, with R σ osc . Thus, the number of particles in the gas can be expressed as In order to describe the confinement of the Bose gas in a thin shell, one is tempted to follow the arguments of the reference [35] and choose a trial wave function of the form We plug this trial function into the above normalization integral (16) and expand the term det g(x 0 , x 1 , x 2 ) in a power series using equation (10). After this expansion, we are able to approximately perform the integrals in the limit of large R, which leads to an exponentially small error O(e −R 2 /σ 2 ) [40, (8.25)]. The integral over x 0 of each term of the power series times a Gaussian is given by for even non-negative integer values of n, while it vanishes for odd non-negative values of n. Therefore, only the even order terms survive and provide results at least of the order of σ 2 /R 2 . Thus, the number of particles can be calculated through the normalization of the 2D function ψ(x 1 , x 2 ) apart from a polynomial error, i.e., But now we argue that a better choice of the trial wave function is provided by where we have defined For x 0 = 0 we have that ψ(0, x 1 , x 2 ) represents the 2D wave function of the gas. In order to calculate the number of particles, we can follow the same procedure as the one above, but now we do not need to perform a Taylor series expansion, since the denominator of the term (21) matches with the one from the volume element. To perform the integral it is only necessary to approximately take the limit of large R, leading to an exponentially small error. Thus, the number of particles is given by the integral of the squared norm of the 2D wave function apart from only an exponentially small error: We conclude that the ansatz (20), (21) is a much better approximation than (17), as the normalization of the wave function in (22) is more accurate than (19). Therefore, we investigate in the following the consequences of (20), (21) in view of reducing the original 3D problem to an effective 2D one.

Reducing Dimensionality
In the previous sections we introduced the necessary differential geometrical notation and the physical ideas on how to confine a weakly interacting Bose-Einstein condensate in the neighbourhood of a manifold. Now we proceed with the description of this system by considering its grand-canonical energy where g int = 4πh 2 a s /M denotes the interaction strength, determined by the s-wave scattering length a s , and µ is the chemical potential of the system. Since the confinement frequency ω is supposed to be large enough, we can follow the same procedure as in the last section and perform this integral only on the manifold neighbourhood N (M). Then, the energy (23) is well approximated by Inserting the trial function (20), (21) into the energy functional (24), we expand it into a Taylor series with respect to x 0 and perform the resulting integral with respect to x 0 approximately in the limit of large R as in section 3. Note that the Laplacian (12) with the trial function (20), (21) reads explicitly Approximately integrating the energy (24) with respect to x 0 , we obtain Based on the ideas of reference [35], we extremize the resulting energy (26) with respect to both ψ * and σ. In the first case we obtain the time-independent 2D Gross-Pitaevskii equation where ψ = ψ(x 0 = 0, x 1 , x 2 ) denotes the 2D wave function, while ∆ M is given in equation (13) in the case when x 0 = 0, and represents an effective potential due to the non-trivial metric. Moreover, g 2D = g int /(σ √ 2π) turns out to be the 2D interaction parameter, which gets larger for smaller values of the width σ. It shows that a strong confinement leads to stronger effective two-particle interactions on the manifold.
Extremizing instead the energy (26) with respect to the cloud width σ yields Thus, in equilibrium, one has to solve both equations (27) and (29) for ψ and σ by taking into account the particle number in equation (22).

Equilibrium on a Sphere
The simplest case to be studied is the ground state of a sphere with radius R. For simplicity, we suppose that we do not have any external potential, i.e., U (x 1 , x 2 ) = 0. Furthermore, from (28) we conclude that the effective potential V eff vanishes for such a sphere. Due to the rotational symmetry of the sphere, the gas in the ground state is described by a uniform distribution, thus the 2D wave function ψ 0 satisfies ∆ M ψ 0 = 0. Moreover, from the normalization (22) we obtain With this the time-independent 2D Gross-Pitaevskii equation (27) reduces to an algebraic equation for the chemical potential µ, i.e., the equation of state where represents the dimensionless interaction strength. Note that P can be tuned by changing the particle number N , the s-wave scattering length a s , as well as by changing the oscillator length σ osc or the radius R. Correspondingly, for a sphere also the ground state thickness σ 0 is uniform, so (29) reduces to In figure 2 we plot the results for the dimensionless Gaussian width σ 0 /σ osc and the dimensionless chemical potential µ/hω as functions of the dimensionless interaction strength P . The width of the Gaussian is positive for any value of P . It coincides with the harmonic oscillator length σ osc for vanishing interactions and increases for repulsive interaction strengths (P > 0). For strong repulsive interactions, its asymptotic behaviour is of the form σ 0 /σ osc = 3 √ P . For attractive interaction strengths (P < 0) it decreases as |P | increases, tending to zero with an asymptotic behaviour of the form σ 0 /σ osc = −1/P .
The dimensionless chemical potential coincides with 1/2 for vanishing interactions and increases for positive values of P . Its asymptotic behaviour for large positive values of P is µ/hω = (5/4)P 2/3 . It reaches zero at about P ≈ −0.44 and turns to be negative for smaller values of P . Its asymptotic behaviour for large negative values of P is given by µ/hω = −3P 2 /4. In order to have some intuitive notion of these dimensionless values, suppose that one is arranging to perform a bubble trap experiment in microgravity with about N = 10 5 rubidium atoms on a sphere with a radius about R = 10 µm and a harmonic trap such that the harmonic oscillator length is of the order of σ osc = 1 µm [41]. As the swave scattering length a s is about 100 times the Bohr radius, we obtain a dimensionless interaction of about P = 2.1. For later figures we always use those experimentally realistic parameters.

Dynamics on a Curved Manifold
In the last section, we calculated equilibrium results for a gas confined on a sphere, finding that the cloud has a positive width for all interactions strengths. In order to determine the stability of the system, we now embark upon a linear stability analysis. To this end, we extend the equilibrium consideration of section 4 and treat the Bose gas dynamically.
We begin a dynamical analysis deriving the temporal evolution of the 2D wave function and of the cloud width. To do that, instead of the energy, we use the corresponding action and a more general ansatz than that of equation (20), which includes an imaginary width term in the Gaussian exponent, according to [42,43] Ψ( Here B represents the variational parameter conjugated to the cloud width, which is necessary to be included in order to properly describe the dynamics of the system. Using the same procedure as in the last section, we insert this ansatz into equation (34) and integrate approximately the variable x 0 . With this we get the following expression for the action: Following the standard approach, we extremize the action (37) with respect to ψ * , σ and B. In this way, we obtain at first the evolution equation for ψ as well as the corresponding equation for σ A subsequent extremization of the action with respect to B yields We remark that in the above three equations (38)- (40), the value of the first variable is fixed x 0 = 0, such that ψ stands for ψ(0, x 1 , x 2 , t). These equations describe the dynamics of a Bose gas on a curved manifold by determining the evolution of the 2D wave function, as well as the real and imaginary cloud width self-consistently.

Collective Modes of BEC on a Sphere
Collective modes of a condensate are of great value for experimental studies, since they allow a quantitative characterization of the underlying system, even when an optical absorption projection is made in the data collection. In addition, collective modes are associated with the equilibrium state around which they occur. Being able to analyse the collective modes of a confined condensate creates the possibility of understanding the influence of various system parameters on the hydrodynamics of the system. In order to determine the low-lying collective modes, we now study small perturbations of the ground state for a Bose gas confined on the surface of a sphere of radius R. To this end, we perform a linear stability analysis of the evolution equations derived in the previous section. Note that the effective potential V eff (x 1 , x 2 ) in equation (38) vanishes for the case of a sphere. We suppose a small perturbation of the ground state in the form where ψ 0 is given in equation (30), the chemical potential µ follows from equation (31), and σ 0 is defined via equation (33).
Inserting the perturbed quantities (41)-(43) into equations (38)- (40) and considering only the first order terms of δψ, δσ and δB, we obtain Here we have used that for a sphere the Laplace-Beltrami operator is proportional to the square of the angular momentum operator L 2 via Note that the eigenvalues of L 2 are given byh 2 l(l+1), for l = 0, 1, 2, ... being the angular momentum quantum numbers and its eigenfunctions are proportional to the spherical harmonics Y lm , with m = 0, ±1, ..., ±l.
The technical details on how to solve equations (44)-(46) are relegated to appendix B. There it is shown that these equations can be straight-forwardly solved by decomposing the functions δψ, δσ, δB for all l = 0, 1, 2... in terms ofȲ lm = Y lm + Y * lm for m = 0, ..., l and proportional toȲ lm = −i(Y lm − Y * lm ) for m = −l, ..., −1. Here we restrict ourselves to summarize and discuss the respective results.
For l = 0, the colective oscillation mode frequency is given by which is of the order of the transversal confinement frequency ω. It coincides with the formula obtained in reference [25], where the mode associated to this frequency was called accordion mode.
For l ≥ 1, irrespective of the sign of the dimensionless interaction strength P we find two branches of collective oscillation mode frequencies, a larger one and a lower one, but degenerate with respect to the magnetic quantum number m. The branch with larger oscillation frequencies is approximately given by while the branch with smaller oscillation frequencies reads approximately Here represent smallness parameters, since the width σ osc is supposed to be much smaller than the radius R of the sphere. As discussed at the end of section 5, we consider σ 2 osc /R 2 = 0.01 to be realistic for a bubble trap in microgravity. A plot of the frequencies Ω l and Λ l as functions of the dimensionless interaction strength P for the lowest values of l can be found in figure 3. These graphics are made for the realistic range of P values, according to the parameters given at the end of section 5.
From equations (33) and (48) we read off that Ω 0 /ω equals to 2 for vanishing interaction and approaches asymptotically to √ 3 for a large dimensionless interaction strength P , as can be seen in panels 3(a) and (b). For l = 1, 2, 3, the frequencies Ω l /ω have a similar behaviour as Ω 0 /ω, but they turn out to be larger than Ω 0 /ω. From the plots we can also see that frequencies Ω l /ω increase with the angular momentum quantum number l.
The dimensionless frequencies Λ l /ω are positive for vanishing interaction, even though they are much smaller than 1. The frequencies Λ l /ω reach zero for some negative value of P and monotonically increase with P , see figure 3(c) and (d). From equation (50) we see that these frequencies decrease for a smaller value of the parameter δ l . On the other hand, for some negative value of the dimensionless interaction strength P the lower frequencies could become imaginary, such that the corresponding solution exhibit an exponential behaviour. We stress that only for quite small negative values of P we still have a stable solution, as can be seen in figure 3(d), and they turn out to be unstable as soon as P is decreased even to relatively moderate values. Note that for l = 0, the lower frequency is not defined due to the conservation of the particle number, as is discussed in more details in appendix B.

Analysis of Modes
Within a linear stability analysis of small perturbations on a sphere, we have derived analytic expressions for two types of collective oscillation mode frequencies Ω l and Λ l , as well as understood their dependences on the angular momentum quantum number l and on the dimensionless interaction strength P . Now, we analyse the respective density profiles of these oscillations on the sphere, whose calculations are relegated to appendix B. We first discuss the accordion mode, which occurs for the angular momentum quantum number l = 0, and afterwards we analyse the modes for larger angular momentum quantum numbers l, and also illustrate some examples. Finally, we discuss in detail the direction of the oscillations, which can be in the confinement direction, along the surface of the sphere, or even have a mixed behaviour. For l = 0 only the mode with frequency (48) appears. In this case, the temporal evolution of each component of the wave function associated with the frequency Ω 0 turns out to be Here r, θ and ϕ denote the spatial variables in terms of the spherical coordinates.
where C 00 is a proportionality constant defined by the intensity of the perturbation, which has to be small. Note that we have Reδψ 00 = 0 in order to satisfy the conservation of the particle number, see appendix B for further details.
To illustrate the accordion mode, the evolution of its density profile given by the squared norm |Ψ 00 (r, θ, ϕ, t)| 2 of the 3D wave function (35) is pictured in figure 4(a)-(c) at different times in x-z plane, for the chosen parameters at the end of section 5. The proportionality constant is chosen to be C 00 = 0.1. The Gaussian width of the state in (b) coincides with the one of the equilibrium state, while the width of the states in (a) and in (c) are, respectively, larger and smaller than the one of the equilibrium state. From the initial state in (a) which has the largest width, the system evolves to the state shown in (b) and finally reaches the state with the smallest width in (c). Then it returns to the state in (b), to the state in (a), and so on. This oscillation happens with frequency Ω 0 . The radial density profile of these three stages for fixed ϕ = 0 and θ = 0 are plotted in figure 4(d). The thinnest stage in green shows a higher peak and the thickest in orange shows a smaller one. This happens since the number of particles N is conserved.
For the cases where l ≥ 1, there are two types of oscillation frequencies Ω l and Λ l , given in equations (49) and (50), respectively. For Ω l , the temporal evolution of the associated components of the wave function are given by  where C 1 lm is a proportionality constant. For Λ l the solutions are δψ lm (t) = C 2 lm ψ 0 with C 2 lm also being a proportionality constant. We illustrate the density profiles of these modes in figure 5 for l = 1 and l = 2, with m = 0. Figures 5(a) and (c) show the density profile |Ψ(r, θ, ϕ, t)| 2 in the x-z plane. Figures 5(b) and (d) show the condensate density on the surface of the sphere, i.e., |Ψ(R, θ, ϕ, t)| 2 . The proportionality constants are chosen to be C 1 10 = C 1 20 = 0.1. In these cases, both the width σ lm (t) = σ 0 + δσ lm (t) and the density of the 2D wave function |ψ lm (t)| 2 = |ψ 0 | 2 + 2|ψ 0 |Reδψ lm (t) turn out to oscillate in time.
From panels 5(a) and (c) we read off that the regions on the sphere where the density maxima are located have the minimal width. The oscillations of δσ change the shape of the Gaussian, similarly to what happens in figure 4(d). The difference is that for the accordion mode this happens in a spherically symmetric way, while for larger values of l there is an angular dependence, as is illustrated in figures 5(a) and (c). Oscillations of Reδψ do not change the width of the Gaussian nor its shape, and only the amplitude of the wave function is changed. For l = 0 such oscillations do not occur, because they would change the number of particles in the system. For l ≥ 1 these oscillations change the distribution of particles along the sphere according to the For the modes with frequencies Ω l given in equation (53), the amplitude of δσ lm is much larger than the amplitude of the real part of δψ lm , since the latter one is proportional to δ l , which is small. In this sense, we say that the density oscillations predominantly occur in the direction perpendicular to the sphere. For the modes corresponding to Λ l given in equation (54), we see that the oscillations of Reδψ lm have a fixed amplitude, while the amplitude oscillations of δσ lm depend on the interaction strength P . This dependence occurs explicitly and implicitly, since both σ 0 and Ω 0 depend on P according to (33) and (48). If P is zero, the amplitude oscillations of δσ lm vanish, so we can say that the mode is in a parallel direction to the sphere. If P increases, the amplitude oscillations of δσ lm increase and it can even be comparable to the amplitude oscillations of Reδψ lm . In this sense, we say that the direction of oscillations turns out to have a mixed behaviour with both the parallel and perpendicular components, which we call a diagonal oscillation.
In order to illustrate the above statements, in figure 6 we plot the normalized vectors in the direction of the vectors (|Reδψ lm |/ψ 0 , δσ lm /σ osc ) for various values of l and P in order to illustrate the difference in the contributions of Reδψ lm and δσ lm in the modes associated to the frequencies Ω l and Λ l , respectively. Note that the amplitude of the collective modes have a degeneracy on the values of m for l fixed. Moreover, the directions of these vectors coincide with the notion given above, of oscillations being in a direction perpendicular, parallel or diagonal to the sphere.
In panel 6(a) the angular momentum quantum number l = 1 is fixed and we vary the dimensionless interaction strength P from 0 to 5. We see that the amplitude oscillations for the modes with larger frequency almost do not change, while the amplitude oscillations for the modes corresponding to the smaller frequency are quite sensitive to the values of P . The amplitude oscillation that is parallel to the sphere for P = 0 turns to be diagonal as P increases. In panel 6(b) the dimensionless interaction strength P = 2.1 is fixed and we vary the angular momentum quantum number l. For the modes corresponding to Ω l , l is varied from 0 to 10, while for the modes corresponding to Λ l , l is varied from 1 to 10, since there is no such mode for l = 0. We see that for both branches the vectors increase their angle with the vertical for increasing l, but this is more pronounced for the modes associated with the frequencies Λ l .
From figure 6 we see that the modes corresponding to the Ω l branch predominantly oscillate in the perpendicular direction, with a negligible dependence on P and a small dependence on l. On the other hand, the modes from the branch corresponding to Λ l depend on both P and l, but have a dominant parallel component. With this, we conclude that the oscillations with higher frequencies are in the direction of the confinement trap, while the oscillations with smaller frequencies occur along the sphere, i.e., in the not confined direction.

Conclusions
Motivated by recent experimental advances in the field, we have studied both the static and the dynamic properties of a Bose-Einstein condensate confined on the surface of a curved manifold. To this end, we provided a general formulation of the problem and derived a self-consistent set of equations for the 2D condensate wave function on the manifold and its width. In particular, we found an effective potential in the resulting 2D Gross-Pitaevskii equation, which depends on the metric of the manifold in a non-trivial way but vanishes for a sphere. For the latter special case we determined in equilibrium how both the width of the condensate and its chemical potential increase with the repulsive interaction strength. Moreover, we found via a linear stability analysis two branches of collective excitations, with distinctly different frequencies. The larger branch frequencies turned out to be of the order of the harmonic confinement frequency and, thus, correspond to oscillations predominantly in the direction of the confinement, i.e, in the perpendicular direction to the sphere. The lower branch frequencies are much smaller and represent oscillations predominantly on the sphere. Our results represent concrete predictions, which presumably could be confirmed in upcoming bubble trap experiments in the NASA Cold Atom Laboratory at the International Space Station [24].
Note that the collective modes analysed in this paper differ from the ones identified in reference [33]. This is due to the fact that our ansatz (35), (36) for the wave function neglects the possibility that the condensate levitates below or above the minimum of the confinement potential.

From the symmetry of the second derivatives and (A.3), (A.4), we then conclude that
2) can be rewritten as s ij = 1 2 (v i ·∂ j v 0 +v j ·∂ i v 0 ), proving the above statement that s ij is symmetric. Now, consider the operator s · g −1 , which is known in the literature as the Gauss map. The directions of its eigenvectors e 1 and e 2 are called principal directions, and these directions correspond to the minimal and maximal curvatures. The matrix representing s · g −1 is given in this basis by where κ 1 and κ 2 denote the curvatures of the manifold in the direction of e 1 and e 2 , respectively. The mean and the Gaussian curvatures of this manifold are defined by Furthermore, there is a third fundamental form [45] defined as It can be shown that this form turns out to be a combination of the previous two, via the relation [45] h ij = −Kg ij + 2Hs ij , (A.8) and the matrix which represents h · g −1 in the principal directions is given by (A.9) Now, let us consider a manifold M(x 0 ) parallel to the manifold M, as explained in section 2. The metric, i.e., the first fundamental form of this manifold M(x 0 ), for |x 0 | < R/2, is defined by and from expression (1) we conclude Combining this equation with the definitions (A.1), (A.2), and (A.7), and taking into account that the fundamental forms are symmetric, we read off that this metric of the manifold M(x 0 ) can be expressed in terms of the three fundamental forms of the manifold M as follows: From the corresponding matrix forms (A.5) and (A.9), it is then possible to deduce Thus, we obtain where R is defined as the minimum mean radius over all points p belonging to M according to equation (3). The square root of this formula can be expanded in a Taylor series, yielding equation (10) from the main text.

Appendix B. Determining the Collective Modes
In this section we present the detailed calculations for the respective results of sections 7 and 8. For the sake of simplicity, we will use a dimensionless form for our equations.
To this end, we define the following dimensionless variables and the dimensionless derivative operators With that, the dimensionless form of the linearized equations (44)-(46) reads where the coefficients k lm , h lm , s lm and b lm are real functions of t for all l = 0, 1, 2, ... and m = 0, ±1, ... ± l. With this, it is possible to decouple the equations of different indices l and m, sinceȲ * lm =Ȳ lm . Thus, the following equations turn out to be much simpler than they would have been with the standard decomposition. Now we insert the ansatz (B.7) into equations (B.3)-(B.5), by taking into account that L 2 Y lm = l(l + 1)Y lm . With this, we obtain that the real part of equation (B.3) is given by Correspondingly, equation (B.4) reduces to Rearranging these equations they turn out to be of the form of a system of linear differential equations: where we introduce the vector X lm (t ) T = (k lm (t ), h lm (t ), s lm (t ), b lm (t )), (B.13) and the matrix Q l = Q 0 + δQ l with and The general solution of this system is given by e Q l t X lm (0), where X lm (0) denotes the initial condition. The problem of computing the exponential matrix e Q l t applied to some vector is standard when the matrix Q l has four different eigenvectors. Let X i , for i = 1, 2, 3, 4, be the eigenvectors associated to the eigenvalues λ i , and X lm (0) to be written in the eigenbasis, i.e., X lm (0) = 4 i=1 v i X i . Then the solution of (B.12) with initial condition X lm (0) is given by X lm (t ) = 4 i=1 v i e λ i t X i . This procedure is going to be used for solving equation (B.12) when l ≥ 1 in appendix B.2 and when l = 0 and P = 0 in the next subsection.
But when l = 0 and P = 0, it turns out that the matrix Q 0 has only three eigenvectors, meaning that the eigenbasis is incomplete. In this case one has to consider generalized eigenvectors, as is shown in detail in the following subsection. To find the solution of (B.12) when l = 0, we have to compute e Q 0 t X 00 (0). The eigenvalues of Q 0 are λ 1± 0 = ±iΩ 0 with multiplicity one and λ 2 0 = 0 with multiplicity two. Here, Ω 0 = Ω 0 /ω and Ω 0 is given in equation (48) of the main text. We discuss now separately the cases P = 0 and P = 0.
When P = 0 the respective eigenvectors are where X 1± 0 are associated to the eigenvalues ±iΩ 0 and X 2± 0 correspond to the eigenvalue 0. Then the initial condition is decomposed according to where v i , for i = 1, 2, 3, 4 are constants, so the exponential matrix X 00 (t ) = e Q 0 t X 00 (0) yields the solution Let us consider now P = 0. In this case it turns out that the matrix Q 0 has only three eigenvectors. In order to have a basis of the four-dimensional vectorial space, we choose the vectors where X 1± 0 are the eigenvectors associated to the eigenvalues ±iΩ 0 , X 2+ 0 is the eigenvector associated to the eigenvalue 0 andX 2 0 is a generalized eigenvector with the property Q 0X 2 0 = X 2+ 0 . Remember that σ 0 = 1 and Ω 0 = 2 when P = 0, meaning that the vectors X consists of the vectors (B.19) as column vectors. With that the exponential of Q 0 t is given by With this, for an initial condition given by the vector the solution of (B.12) reads . (B.23) Thus, there is a secular term in the last term of the right-hand side, which leads to linear growth in time. In the following we show that this solution is eliminated due to the conservation of the number of particles, which is described by are real vectors. With this we obtain which is a real solution. Analogously, we could compute the solution e Q 0 t i , but it turns out to lead to the same dynamics as above apart from a phase. So, we will consider only the solution (B.31) which can be written as with the initial condition X 0 (0) T = (0, 0, C 00 , 0), (B.33) where C 00 is a proportionality constant. This solution means that, when one performs at t = 0 a small perturbation of the ground state given by then the temporal evolution of this perturbation is given by δσ 00 (t ) = C 00 cos(Ω 0 t )Ȳ 00 , (B.35) These equations are the dimensionless form of equations (52) in the main text.
Appendix B.2. Solving Differential Equations for l ≥ 1 In this section we solve equation (B.12) for l ≥ 1. The solution of this system is given by where X lm (0) denotes the initial state. In order to evaluate the matrix exponential e Q l t we have to determine at first the eigenvalues of the matrix Q l . They turn out to be the roots of the polynomial The roots λ 1± l are always imaginary, so they correspond to oscillating solutions. Instead, the roots λ 2± l are imaginary for positive values of the radicand, leading also to oscillating solutions, but they are real when the radicand is negative, leading to exponentially increasing and decreasing solutions. The radicand is positive for positive and small negative values P , but it turns out to be negative for most negative values of P . Thus, we can conclude that the system is stable for positive interactions and small enough negative interaction strengths, while it becomes unstable for most negative interaction strengths.
The frequencies of oscillation are given by Ω l = Ω l /ω = |Re(iλ 1± l )| and Λ l = Λ l /ω = |Re(iλ 2± l )|, and read up to the first order of δ l Ω l = Ω 0 + 1 Ω  with the initial condition Analogously, we also obtain with the corresponding initial condition