Ground state of weakly repulsive soft-core bosons on a sphere

TL;DR

This paper investigates the ground state phases of weakly interacting soft-core bosons on a spherical surface, revealing a transition from fluid to various cluster-based supersolid phases as density increases.

Contribution

It introduces a variational approach to analyze finite-size phase crossovers and cluster formations in a spherical geometry for weakly repulsive bosons.

Findings

Quantum fluctuations stabilize the fluid phase at low density.

Compression leads to inhomogeneous cluster phases.

All identified cluster phases exhibit supersolidity.

Abstract

We study a system of penetrable bosons embedded in a spherical surface. Under the assumption of weak interaction between the particles, the ground state of the system is, to a good approximation, a pure condensate. We employ thermodynamic arguments to investigate, within a variational ansatz for the single-particle state, the crossover between distinct finite-size "phases" in the parameter space spanned by the sphere radius and the chemical potential. In particular, for radii up to a few interaction ranges we examine the stability of the fluid phase with respect to a number of crystal-like arrangements having the symmetry of a regular or semi-regular polyhedron. We find that, while quantum fluctuations keep the system fluid at low density, upon compression it eventually becomes inhomogeneous, i.e., particles gather together in clusters. As the radius increases, the nature of the high-density aggregate varies and we observe a sequence of transitions between different cluster phases ("solids"), whose underlying rationale is to maximize the coordination number of clusters, while ensuring at the same time the proper distance between each neighboring pair. We argue that, at least within our mean-field description, every cluster phase is supersolid.

Figures

Tables

Table 1. Table 1: Ratio between the edge length ℓ ℓ \ell and the circumscribed radius R 𝑅 R for the polyhedra depicted in Fig. 3. For the snub-cube case, t = ( 1 + 19 − 3 ​ 33 3 + 19 + 3 ​ 33 3 ) / 3 𝑡 1 3 19 3 33 3 19 3 33 3 t=\left(1+\sqrt[3]{19-3\sqrt{33}}+\sqrt[3]{19+3\sqrt{33}}\right)/3 and β = 26 + 6 ​ 33 3 𝛽 3 26 6 33 \beta=\sqrt[3]{26+6\sqrt{33}} . In the last two lines, the quoted ℓ / R ℓ 𝑅 \ell/R refers to the biscribed form of the polyhedron, and ℓ ℓ \ell is the short edge (the most numerous one).

polyhedron	$\ell /R$	$D=2R^2/{\ell }^2$
tetrahedron	$2\sqrt{2/3}=1.632\dots$	$3/4$
cube	$2\sqrt{3}/3=1.154\dots$	$3/2$
octahedron	$\sqrt{2}=1.414\dots$	$1$
dodecahedron	$4/\left(\sqrt{3}+\sqrt{15}\right)=0.713\dots$	$\left(9+3\sqrt{5}\right)/4=3.927\dots$
icosahedron	$4/\sqrt{10+2\sqrt{5}}=1.051\dots$	$\left(5+\sqrt{5}\right)/4=1.809\dots$
cuboctahedron	$1$	$2$
rombicuboct.	$2/\sqrt{5+2\sqrt{2}}=0.714\dots$	$\left(5+2\sqrt{2}\right)/2=3.914\dots$
snub cube	$\sqrt{2\left(2-8/\beta +\beta \right)/3}/\sqrt{t^2+t^{-2}+1}=0.744\dots$
tetrakis hex.	$\sqrt{6\left(3-\sqrt{3}\right)}/3=0.919\dots$	$\left(3+\sqrt{3}\right)/2=2.366\dots$
pentakis dod.	$\sqrt{30\left(15-\sqrt{15\left(5+2\sqrt{5}\right)}\right)}/15=0.640\dots$

Equations

$$u(x)=v(R\arccos x)\mathrm{or}v(s)=u\left(\cos \frac{s}{R}\right).$$

$$\Psi ({\Omega }_1,\dots ,{\Omega }_N)=\prod _{i=1}^N\psi ({\Omega }_i),$$

$$-\frac{{\hslash }^2}{2m}{\nabla }^2\psi +(N-1)\int {\mathrm{d}}^2{\Omega }^{\prime }|\psi ({\Omega }^{\prime }){|}^{2}u(\hat{\mathbf{r}}\cdot {\hat{\mathbf{r}}}^{\prime })\psi (\Omega )=\lambda \psi (\Omega ),$$

$${\nabla }^2=\frac{1}{R^2}\left\{\frac{1}{\sin \theta }\frac{\partial }{\partial \theta }\left(\sin \theta \frac{\partial }{\partial \theta }\right)+\frac{1}{{\sin }^2\theta }\frac{{\partial }^2}{\partial {\varphi }^2}\right\}$$

$$\lambda =-\frac{{\hslash }^2}{2m}\int {\mathrm{d}}^2\Omega {\psi }^{\ast }{\nabla }^2\psi +(N-1)\int {\mathrm{d}}^2\Omega {\mathrm{d}}^2{\Omega }^{\prime }|\psi (\Omega ){|}^{2}u(\hat{\mathbf{r}}\cdot {\hat{\mathbf{r}}}^{\prime })|\psi ({\Omega }^{\prime }){|}^2.$$

$$\mathcal{E}[\psi ]=-\frac{{\hslash }^2}{2m}\int {\mathrm{d}}^2\Omega {\psi }^{\ast }{\nabla }^2\psi +\frac{N-1}{2}\int {\mathrm{d}}^2\Omega {\mathrm{d}}^2{\Omega }^{\prime }|\psi (\Omega ){|}^{2}u(\hat{\mathbf{r}}\cdot {\hat{\mathbf{r}}}^{\prime })|\psi ({\Omega }^{\prime }){|}^2.$$

$$\psi (\Omega )=\sum _{l=0}^{\infty }\sum _{m=-l}^l{c}_{lm}{Y}_l^m(\Omega )\mathrm{with}\sum _{lm}|c_{lm}{|}^2=1$$

$${\mathcal{E}}_{\mathrm{kin}}=\frac{{\hslash }^2}{2m{R}^2}\sum _{lm}l(l+1)|c_{lm}{|}^2.$$

$${\mathcal{E}}_{\mathrm{pot}}$$

$$\mathcal{E}=\frac{N-1}{4}{\int }_{-1}^1\mathrm{d}xu(x)=\frac{N-1}{4}{\int }_{-1}^1\mathrm{d}xv(R\arccos x)=\frac{N-1}{8\pi }\int {\mathrm{d}}^2\Omega v(R\theta )⟶\frac{\rho }{2}\tilde{v}(0),$$

$$R=\frac{\ell }{4}\sqrt{10+2\sqrt{5}},$$

$$\psi (\hat{\mathbf{r}})=C_{\alpha }\sum _{k=1}^n\exp \left\{-\alpha {\left(\frac{R\hat{\mathbf{r}}-{\mathbf{R}}_k}{\ell }\right)}^2\right\}=C_{\alpha }\sum _{k=1}^n\exp \left\{-\frac{2R^2}{{\ell }^2}\alpha (1-\hat{\mathbf{r}}\cdot {\hat{\mathbf{R}}}_k)\right\},$$

$$c_{lm}=\int {\mathrm{d}}^2\Omega Y_l^{m\ast }(\Omega )\psi (\Omega )$$

$$T_l^n(\Omega )=\sum _{m=-l}^l{b}_{nlm}{Y}_l^m(\Omega ).$$

$$T_0^0$$

$$T_{10}^0$$

$$\psi =C_{\beta }(T_0^0+\beta T_6^0)\left(\mathrm{with}{C}_{\beta }=\frac{1}{\sqrt{1+{\beta }^2}}\right)$$

$${\mathcal{E}}_{\mathrm{kin}}=21{\beta }^2{C}_{\beta }^2\frac{{\hslash }^2}{m{R}^2},$$

$$\frac{N-1}{2}{C}_{\beta }^4\int {\mathrm{d}}^2\Omega {\mathrm{d}}^2{\Omega }^{\prime }{\left(T_0^0+\beta T_6^0(\hat{\mathbf{r}})\right)}^{2}u(\hat{\mathbf{r}}\cdot {\hat{\mathbf{r}}}^{\prime }){\left(T_0^0+\beta T_6^0({\hat{\mathbf{r}}}^{\prime })\right)}^2.$$

$$\int {\mathrm{d}}^2\Omega {\mathrm{d}}^2{\Omega }^{\prime }(T_0^0{)}^{2}u(\hat{\mathbf{r}}\cdot {\hat{\mathbf{r}}}^{\prime })(T_0^0{)}^2=\frac{1}{2}{\int }_{-1}^1\mathrm{d}xu(x);$$

$$\int {\mathrm{d}}^2\Omega {\mathrm{d}}^2{\Omega }^{\prime }2\beta T_0^0{T}_6^0(\hat{\mathbf{r}})u(\hat{\mathbf{r}}\cdot {\hat{\mathbf{r}}}^{\prime })(T_0^0{)}^2=0;$$

$$\int {\mathrm{d}}^2\Omega {\mathrm{d}}^2{\Omega }^{\prime }{\beta }^2(T_6^0(\hat{\mathbf{r}}){)}^{2}u(\hat{\mathbf{r}}\cdot {\hat{\mathbf{r}}}^{\prime })(T_0^0{)}^2=\frac{{\beta }^2}{2}{\int }_{-1}^1\mathrm{d}xu(x);$$

$$\int {\mathrm{d}}^2\Omega {\mathrm{d}}^2{\Omega }^{\prime }2\beta T_0^0{T}_6^0(\hat{\mathbf{r}})u(\hat{\mathbf{r}}\cdot {\hat{\mathbf{r}}}^{\prime })2\beta T_0^0{T}_6^0({\hat{\mathbf{r}}}^{\prime })=2{\beta }^2{\int }_{-1}^1\mathrm{d}xu(x)P_6(x).$$

$$I_1=\int {\mathrm{d}}^2\Omega {\mathrm{d}}^2{\Omega }^{\prime }{\beta }^2(T_6^0(\hat{\mathbf{r}}){)}^{2}u(\hat{\mathbf{r}}\cdot {\hat{\mathbf{r}}}^{\prime })2\beta T_0^0{T}_6^0({\hat{\mathbf{r}}}^{\prime })$$

$$I_2=\int {\mathrm{d}}^2\Omega {\mathrm{d}}^2{\Omega }^{\prime }{\beta }^2(T_6^0(\hat{\mathbf{r}}){)}^{2}u(\hat{\mathbf{r}}\cdot {\hat{\mathbf{r}}}^{\prime }){\beta }^2(T_6^0({\hat{\mathbf{r}}}^{\prime }){)}^2.$$

$$\left(\begin{array}[]{ccc}6&6&6\\ m&m_{1}&-m_{2}\end{array}\right)$$

$$I_1=\frac{20\sqrt{11\cdot 13}}{17\cdot 19}{\beta }^3{\int }_{-1}^1\mathrm{d}xu(x)P_6(x).$$

$$\left(\begin{array}[]{ccc}l&6&6\\ 0&0&0\end{array}\right)\,,\,\,\,\left(\begin{array}[]{ccc}l&6&6\\ m&m_{1}&-m_{2}\end{array}\right)\,,\,\,\,\left(\begin{array}[]{ccc}l&6&6\\ -m&m_{3}&-m_{4}\end{array}\right)\,.$$

$$m=m_1=m_2=0;$$

$$m=0,m_1=m_2=5;$$

$$m=0,m_1=m_2=-5;$$

$$m=5,m_1=-5,m_2=0;$$

$$m=5,m_1=0,m_2=5;$$

$$m=-5,m_1=5,m_2=0;$$

$$m=-5,m_1=0,m_2=-5;$$

Full text

I Introduction

In the last few decades, thanks to the continued advance in the preparation and manipulation of ultracold atomic gases, the investigation of quantum correlation effects has reached a level of sophistication which would have been simply unimaginable before Dalfovo ; Leggett1 ; Bloch . By confining atoms in optical and magnetic traps also the properties of low-dimensional quantum systems can be analyzed in detail, making it possible to test theoretical predictions and approximations (see, e.g., Kinoshita ). Even the range and strength of interatomic forces can be tuned to a certain extent Buechler (e.g., by the technique of Rydberg dressing), which has ultimately allowed to realize atomic systems characterized by an effective isotropic repulsion saturating to a finite value for zero separation Henkel ; Pupillo .

Finite-strength interactions are frequently encountered in the classical realm as models for polymer coils or dendrimers dispersed in a good solvent (see, e.g., Refs. Louis ; Mladek ). The phase behavior of such fluids can be very rich, featuring in equilibrium any sort of mesoscopic structures (like clusters, micelles, and filaments — see Rendiconti and references cited therein). For purely repulsive particles, a distinctly universal behavior emerges at low temperature Likos , where, depending on the shape of the potential, the high-density phase is either fluid or cluster solid. The latter phase can be described as a crystalline system with multiply-occupied cells, each hosting on average the same number of particles (see examples in Refs. Zhang ; Prestipino1 ; Prestipino2 ). Cluster crystals are characterized by a marked single-particle diffusion Moreno , which keeps the interstitial density at a non-zero value in equilibrium. It is clear that clustering, as a self-assembly phenomenon, can only occur when the formation of bunches of fully overlapping particles is energetically preferred over diffuse partial overlap Mladek .

In recent years, various quantum models of softly-repulsive bosons have been considered, whose phase diagram was worked out at zero temperature ($T=0$) both in mean field (MF) Henkel2 ; Ancilotto ; Kunimi ; Prestipino3 ; Prestipino4 and by Monte Carlo (MC) simulation Saccani ; Cinti1 ; Cinti2 ; Macri . In these systems the fluid-to-solid transition is the necessary outcome of the softening of roton-like modes in the fluid. The mechanism promoting quantum crystallization at $T=0$ is different from the freezing of hard-core fluids at high $T$, which is typically an entropy-driven (rather than an energy-driven) phase transition. Moreover, quantum cluster crystals may be supersolid, a feature which lacks a classical analog. Supersolidity has to do with an anomalous decrease of rotational inertia Leggett2 ; Sepulveda ; Kuklov ; Boninsegni , as if a fraction of the system remains stationary when the crystal is set into slow rotation around an axis.

In the present study we explore by MF theory the low-temperature physics of penetrable bosons in a setting which apparently has not been considered so far, i.e., that provided by confinement to a spherical surface. The systems which more closely resemble our model system are ultracold dilute gases trapped in a thin spherical shell, which have been the subject of a few experimental studies Zobay ; Garraway . Other examples of real systems bearing some similarity to our model are multi-electron bubbles in liquid helium Tempere , arrangements of protein subunits on spherical viruses Zandi , and colloidal particles in colloidosomes Fantoni . In such systems the interparticle forces depend on the Euclidean distance rather than on the arc-length distance, which is the intrinsic metric for particles embedded in the surface of a sphere. However, this difference is immaterial as long as in our theory the pair potential is expressed in terms of the angular separation between the particles (see Section II). Spherical boundary conditions have often been used in numerical experiments Prestipino5 ; Prestipino6 ; Prestipino7 ; Guerra ; Post ; Vest ; Bozic as a means to discourage crystalline ordering at high density (as well known, triangular order is frustrated on a sphere). In practice, the sphere curvature imposes a distinct excess of fivefold disclinations over sevenfold ones, which considerably complicates the search for optimal packings, even for small radii. Very recently, Franzini et al. Franzini have studied by density-functional theory a system of classical particles interacting through a generalized-exponential repulsion (GEM-4), finding a rich catalog of cluster phases as a function of the sphere radius $R$.

It is reasonable that, as the spherical surface gets more and more filled with particles, it will be found more convenient also for a quantum system of penetrable disks to clusterize, thus becoming solid-like inhomogeneous; moreover, as for a classical system, it is likely that numerous aggregates will compete for stability as a function of $R$. The most symmetric ones, i.e., those sharing the symmetries of a regular or semi-regular circumscribable polyhedron, are natural candidates for the high-density phases. There is a limiting case where the theoretical analysis of the quantum system at $T=0$ is greatly simplified, that is weak interparticle repulsion. Then, MF theory becomes an effective method, as practically demonstrated for a specific instance of soft-core bosons by the “exact” phase diagram reported in Cinti2 . As already made in Prestipino3 ; Prestipino4 , we further simplify our treatment using an educated guess of the condensate wave function, to be optimized by the variational method. By taking advantage of a well-established theoretical framework, we aim to gain insight into the self-organization principles underlying structure selection in a quantum many-body system characterized by a wealth of possible ground states.

The outline of the paper is as follows. In Sec. II we introduce the model and outline the variational MF theory employed to study its ground-state behavior. To give a flavour of the results obtained, in Sec. III we work out analytically a simpler exercise, which is nonetheless capable to predict the onset of clusters at high density in a specific range of $R$ values. Afterwards, in Sec. IV we present the full phase diagram of the system as a function of $R$ and chemical potential. In Sec. V we devote special attention to the issue of supersolidity of the spherical cluster phases. We show that, within our theory, all such phases are indeed supersolid. Concluding remarks are postponed to Sec. VI.

II Model and theory

We investigate a system of $N$ identical spinless bosons, living on a sphere of radius $R$ and interacting with each other via a bounded potential $v(s)$, function of the arc-length distance $s$. A paradigmatic case of finite repulsion is the penetrable-sphere model (PSM) potential, $v(s)=\epsilon\vartheta(\sigma-s)$, $\vartheta$ being the Heaviside step function (PSM bosons will be our favourite case study later). It is convenient to introduce another parametrization of the potential, written in terms of the scalar product between the unit vectors $\hat{\bf r}$ and $\hat{\bf r}^{\prime}$ identifying the positions on the sphere of the interacting pair. Using $\hat{\bf r}\cdot\hat{\bf r}^{\prime}=\cos(s/R)$, we define:

[TABLE]

If the interaction potential were given in terms of the 3D Euclidean distance $r$, the definition of $u$ in (2.1) would be modified into $u(x)=v(2R\sqrt{(1-x)/2})$, but no change will occur in the subsequent analysis.

It is not a priori obvious how to quantize a system of particles living in a curved space. Canonical quantization rules are inconsistent and a way out is to quantize angular momentum directly — see this point thoroughly discussed in Kleinert . In case of a free particle on a sphere, this entails taking the Hamiltonian (kinetic energy) to be $L^{2}/(2mR^{2})=-\hbar^{2}/(2m)\nabla^{2}$, where $m$ is the particle mass and $\nabla^{2}$ is the Laplace-Beltrami operator on the sphere (see Eq. (2.4) below). This same approach was followed by many authors Ezra ; Seidl ; Loos ; Yang .

In the MF (Hartree) approximation, particles are treated as they were independent of each other and the $N$-boson ground state is therefore a pure condensate:

[TABLE]

where $\Omega_{i}=(\theta_{i},\phi_{i})$ are the angular coordinates of the $i$-th particle, i.e., the spherical variables specifying its 3D position ${\bf r}_{i}=R\hat{\bf r}_{i}$ (in the following, $\Omega$ and $\hat{\bf r}$ are used interchangeably as argument of $\psi$). The best choice of single-particle wave function is that minimizing the expectation value of the Hamiltonian in $\Psi$, which corresponds to a normalized function obeying the (time-independent) Gross-Pitaevskii (GP) equation Gross1 ; Pitaevskii1 ; Gross2 (see Appendix A):

[TABLE]

where

[TABLE]

is the spherical Laplacian and

[TABLE]

The value of $\lambda$ is consistent with $\psi$ being a solution to (2.3): indeed, multiplication of both sides of (2.3) by $\psi^{*}$ and subsequent integration over angles immediately leads to (2.5) under the assumption of unit norm for $\psi$. In Appendix A, we discuss problems related with the numerical solution of Eq. (2.4). We argue that a faster and physically more grounded method, which has proved effective in our exploration of the ground state of the planar system, is to minimize the following MF energy functional (kinetic energy per particle plus potential energy per particle) via the optimization of a parametric wave function:

[TABLE]

For a short-range potential, the ground-state energy in the planar limit $R\gg\sigma$ is only controlled by the dimensionless quantity $\rho\sigma^{2}\epsilon/e_{0}$ (where $e_{0}=\hbar^{2}/(m\sigma^{2})$), which we hereafter refer to as the “density”. When $R$ is finite the sphere radius is an additional control parameter, i.e., the properties of the system depend separately on $R$ and $\rho$. In the following, we take $\sigma$ and $e_{0}$ as units of length and energy, respectively.

We now provide a more explicit expression of ${\cal E}[\psi]$ that applies for any normalized wave function written as an expansion in spherical harmonics:

[TABLE]

(if $\psi(\Omega)$ is real then $c_{l,-m}=(-1)^{m}c_{lm}^{*}$). Computing the kinetic energy is straightforward; since $R^{2}\nabla^{2}Y_{l}^{m}=-l(l+1)Y_{l}^{m}$, we readily find:

[TABLE]

As to the potential energy ${\cal E}_{\rm pot}$ (second term in Eq. (2.6)), in Appendix B we derive the following result:

[TABLE]

where $P_{l}(x)$ are Legendre polynomials and the matrices are Wigner 3-j symbols.

The total energy per particle is the sum of (2.8) and (LABEL:eq2-9). In practice, the $l$ sum must be truncated, i.e., $l\leq l_{\rm max}$, where $l_{\rm max}$ is chosen in accordance with the spatial resolution adopted for the description (see this point discussed, e.g., in Ref. Prestipino8 ). To check consistency, let us consider the homogeneous fluid, corresponding to $c_{lm}=\delta_{l0}\delta_{m0}$ or $\psi=Y_{0}^{0}=1/\sqrt{4\pi}$. It then follows from Eqs. (2.8) and (LABEL:eq2-9), as well as directly from Eq. (2.6), that

[TABLE]

where the last step follows in the planar limit, that is for $N\rightarrow\infty,\,R\rightarrow+\infty$, and $N/(4\pi R^{2})\rightarrow\rho$. As expected, the limiting value of ${\cal E}$ is nothing but the specific energy of the planar fluid Prestipino3 .

III Cluster formation at high density: Proof of concept

As argued in the Introduction, in a spherical quantum system the most stable $T=0$ configuration would not necessarily be fluid. Depending on the radius $R$, other phases may be expected to arise as ground states in a system of softly-repulsive particles at high density. In particular, we guess a primary role for cluster-crystal-like arrangements having the symmetry of a regular (i.e., Platonic) or semi-regular (i.e., Archimedean or circumscribable Catalan) polyhedron. It is easy to conjecture that the stable phase at high density would crucially depend on the value of $R$, since the latter determines the edge length $\ell$ of the inscribed polyhedron and consequently also the geodesic distance between two neighboring clusters. Considering that for the PSM interaction the edge of the triangular-crystal lattice is about $1.51\sigma$ at melting Prestipino3 , we expect that the structure of the high-density phase will be found among those polyhedra having $\ell\lesssim 1.51\sigma$. For instance, since the sphere circumscribing the regular icosahedron has a radius of

[TABLE]

an icosahedral cluster phase is most likely to occur for $R\approx 1.4\sigma$.

A real, one-parameter form of $\psi$ adequate to represent the pattern expected at large $\rho$ on the sphere is a sum of Gaussians centered at the vertices ${\bf R}_{k}$ $(k=1,\ldots,n)$ of the inscribed polyhedron:

[TABLE]

where $\alpha$ is a variational parameter and $C_{\alpha}$ is a (real) normalization constant (notice that $R/\ell$ is a pure number, specific of the given polyhedron). The fluid ground state ($\psi=1/\sqrt{4\pi}$) is recovered from (3.2) as a limiting case, i.e., for $\alpha=0$. In a solid-like system, $\alpha>0$ represents the inverse square width of the local-density peaks. We point out that a real $\psi$ is not a limitation whatsoever; indeed, we show in Appendix C that the true single-particle wave function of minimum energy is, as already known from three dimensions Pomeau , a real function.

We now describe the method followed to draw the “melting line” of a given cluster phase as a function of $R$ for $T=0$. To this purpose we employ a thermodynamic framework, with the idea that when a crossover (a rounded phase transition) occurs at fixed $R$ from one ground state to the other the number of particles is very large (we shall see a posteriori that this is always a safe assumption). In this respect, the chemical potential $\mu$ is a more meaningful control parameter than the pressure $P$ since the surface area is fixed. In brief, we first determine the energy per unit particle $e$ as a function of $\rho$, taking $N=4\pi R^{2}\rho$ in Eq. (2.6). Called $\overline{\alpha}(\rho)$ the point of absolute minimum of ${\cal E}([\psi(\alpha)];\rho)$ for the fixed $\rho$, we have $e(\rho)={\cal E}([\psi(\overline{\alpha}(\rho)];\rho)$. Once the energy has been computed, the transition point $\mu_{c}(R)$ is where the fluid and the solid have the same grand potential per unit area (i.e., where the minimum of $\rho(e(\rho)-\mu)$ is the same for both phases). In the fluid phase, where Eq. (2.19) holds, the relation between $\mu$ and $\rho$ is thus found to be $\mu=(2+(N-1)^{-1}){\cal E}(\rho)\simeq 2{\cal E}(\rho)$.

In order to compute ${\cal E}[\psi]$, two different roads can be followed: either we evaluate Eq. (2.6) numerically, or we attempt an estimate of the Fourier coefficient

[TABLE]

for all $l\leq l_{\rm max}$, and then use Eqs. (2.8) and (LABEL:eq2-9). Indeed, in the following we will pursue both routes; however, before that we show the feasibility of our approach by providing a fully analytic demonstration of clusterization in a system of spherical bosons at $T=0$.

For our proof we make use of a variational wave function simpler than (3.2), but still endowed with the symmetries of the high-density phase we aim to describe (we will focus on the icosahedral cluster phase). In this regard, it is useful to recall an important paper by Zheng and Doerschuk Zheng where they explain how to construct a basis in the subspace of square-integrable $\Omega$ functions that are invariant under every rotation of the icosahedral group (see related comments on this subject at Baez ; Egan ). These basis functions, denoted $T_{l}^{n}(\Omega)$ and dubbed icosahedral harmonics, are real and orthonormal, and given by

[TABLE]

For fixed $l$ there are $N_{l}$ icosahedral harmonics $T_{l}^{n}$ ($n=0,\ldots,N_{l}-1$) that are linear combinations of the $Y_{l}^{m}$ for $m=-l,\ldots,l$ (hence $N_{l}\leq 2l+1$). Zheng and Doerschuk have derived recursive formulae for the coefficients $b_{nlm}$, for arbitrary $n,l,m$, including the cases where $N_{l}>1$ (which only occurs for $l\geq 30$). In particular, the first three icosahedral harmonics turn out to be:

[TABLE]

Of particular interest to us is the function $T_{6}^{0}$, which has twelve maxima of same height centered at the vertices of a regular icosahedron (see Fig. 1). Hence, our problem becomes one of determining whether, at sufficiently high density, the energy of

[TABLE]

reaches its minimum for some $\beta>0$. In this circumstance, the fluid phase (represented by $T_{0}^{0}$) is doomed to transform upon compression into an icosahedral cluster phase. The only caveat is that $\beta<\beta_{\rm max}\simeq 0.7$ in (3.6), if we want to exclude the appearance of spurious maxima in $\psi^{2}$.

While it is immediate to compute the kinetic energy per particle (cf. Eq. (2.8)), equal to

[TABLE]

it is much harder to determine ${\cal E}_{\rm pot}$, given by the double integral

[TABLE]

After expanding each square, Eq. (3.8) becomes the sum of nine terms, not all distinct, most of which can be simply evaluated by using the orthonormality property of spherical harmonics:

[TABLE]

More cumbersome is the calculation of

[TABLE]

and

[TABLE]

In the former case, we are required to compute 3-j symbols of the kind

[TABLE]

with $m,m_{1},m_{2}=0,\pm 5$. The only non-zero symbols are those for which $m+m_{1}-m_{2}=0$, which can occur in one of seven ways. The end result is

[TABLE]

As for $I_{2}$, we directly start from Eq. (LABEL:eq2-9) with $\psi=\beta T_{6}^{0}$. We now need to compute the following 3-j symbols:

[TABLE]

The first symbol is non-zero exclusively for $l$ even and not larger than 12. For each allowed $l$, the second symbol in Eq. (3.14) is non-zero in at most nine cases (three cases for $l=0,2,4$, seven for $l=6,8$, and nine for $l=10,12$):

[TABLE]

and similar considerations apply for the third 3-j symbol.

After a lengthy series of steps we are eventually led to:

[TABLE]

with

[TABLE]

We underline that the individual integrals contributing to the total energy have been successfully checked, for a few $\beta$ values, against MC integration (see Sec. IV).

For the PSM potential,

[TABLE]

it holds:

[TABLE]

Using the property (valid for any $l>0$):

[TABLE]

we easily obtain:

[TABLE]

In other words, the PSM energy for $\psi=C_{\beta}(T_{0}^{0}+\beta T_{6}^{0})$ admits an explicit expression in closed form.

We show in Fig. 2 some data obtained from Eqs. (3.7), (3.16), (LABEL:eq3-17), and (3.21). The plotted quantity is the excess energy $\Delta{\cal E}={\cal E}[\psi(\beta)]-{\cal E}_{f}$ with ${\cal E}_{f}=(N-1)E_{0}/2$, namely the energy of the solid relative to the fluid, which has been computed for $R=1.45\sigma$ and three distinct values of $\rho$. As $\rho$ increases, the solid energy falls eventually below the fluid energy, implying a transition from the fluid to the icosahedral phase upon varying $\mu$ at fixed $R$ (in this case, each cluster hosts roughly 30 particles at melting). This transition has the nature of a first-order phase change, accompanied by metastability of both phases beyond the transition point. It is worth noting the resemblance of this phenomenon to the onset of icosahedral ordering of disordered aggregates of disclinations in a system of hard calottes on a sphere Prestipino5 ; Prestipino6 ; Guerra . Even in that case a geometric pattern emerges when tuning a control parameter (the density), although its mechanism is purely entropic rather than energy-promoted as in the present case.

IV Results

The main value of the simplified variational calculation carried out in Sec. III is to make it evident that in a system of weakly-repulsive spherical bosons a sharp crossover occurs at $T=0$, as a function of $\mu$ and for $R\approx 1.4\sigma$, from the fluid to a cluster phase of icosahedral symmetry. We emphasize that the choice of the icosahedron as supporting frame for the clusters is just one possibility; in fact, as $R$ increasingly departs from $1.4\sigma$, other polyhedra will be better suited than the icosahedron to match the condition $\ell\lesssim 1.51\sigma$. In Table I we report the value of $X=\ell/R$ for ten different solids, depicted in Fig. 3, which we have selected among tens of regular or semi-regular circumscribable polyhedra VisualPolyhedra as the reference structures that likely underlie the cluster phases for $R$ up to $\approx 2.5\sigma$; each such geometry would become relevant in a $R$ interval centered about $1.51\sigma/X$. Besides the five regular (Platonic) solids, we focus our attention on three Archimedean solids and two biscribed Catalan solids (we recall that Archimedean solids have regular faces — not all of the same type — meeting in identical vertices, while Catalan solids are dual to Archimedean solids and not all vertices are equivalent). In selecting the solids in Table I our criterion was to rule out all semi-regular polyhedra with too many vertices (more than 32) or too large faces (which would correspond to big surface “holes” devoid of particles). Clearly, we have no argument to exclude that other structures will also come into play (in fact, we have good reasons to think that some low-symmetry structures are actually relevant, see more below), but there is anyway no hope to identify all local minima in what is likely to be a rugged free-energy landscape.

To draw the $T=0$ phase diagram of the spherical PSM system by the variational method, we need to evaluate the energy per particle (2.6) with $\psi$ given as in Eq. (3.2). A viable method is to resort to MC integration. For any of the structures in Table I, the specific energy $\cal E$ at fixed $R$ is immediately obtained for any $\rho$ once the kinetic energy per particle ${\cal E}_{\rm kin}$ and the potential energy per pair, $2{\cal E}_{\rm pot}/(N-1)$, are given. Using a standard algorithm to generate points $\hat{\bf r}_{i}$ distributed at random over the unit sphere Krauth , we can estimate the two above-cited energies by the following formulae:

[TABLE]

with $u_{ij}=u(\hat{\bf r}_{i}\cdot\hat{\bf r}_{j})\,,\psi_{i}=\sum_{k=1}^{n}\exp\left\{-D\alpha(1-\hat{\bf r}_{i}\cdot\hat{\bf R}_{k})\right\}$, and

[TABLE]

In the above expressions, $n$ denotes the number of vertices of the given polyhedron while $D$ is a shorthand for $2R^{2}/\ell^{2}$. To make the error on the kinetic-energy estimate comparable to that on the potential energy, the total number ${\cal N}_{\rm tot}$ of points $i$ in the first of Eqs. (4.1) has been taken equal to the number of random pairs in the second (typically ${\cal N}_{\rm tot}\approx 10^{10}$).

Figure 4 gives an idea of the type of results obtained. These data are relative to a pure condensate with snub-cube symmetry and refer to $R=1.9\sigma$. On the left panel of Fig. 4, we have plotted the excess energy $\Delta{\cal E}={\cal E}-{\cal E}_{f}$ as a function of $\alpha$ for a few densities; we see that, for $\rho\gtrsim 14$, a minimum develops for a non-zero $\alpha$, signaling the onset of an inhomogeneous phase at high density. On the right panel, the chemical potential has been adjusted so that the grand potential of the fluid (hence, its thermodynamic pressure) equals that of the cluster phase. This condition defines the snub-cube transition point $\mu_{c}$ for the given $R$, whereas the abscissae of the two equal minima are the coexistence densities, $\rho_{f}$ and $\rho_{s}$. For this case, the average number of particles per cluster at melting is $4\pi R^{2}\rho_{s}/24=30.38$.

Before presenting the full phase diagram we introduce an alternative method to draw the function $\Delta{\cal E}(\alpha)$ for fixed values of $R$ and $\rho$, which may also serve to check consistency with MC data. This method is fully analytic and consists in reconstructing the variational energy through the exact calculation of a sufficiently large number of $\psi$ and $\psi^{2}$ modes (cf. Eq. (A.16)). Once the latter quantities have been computed, the energy will be determined as a function of $\alpha$ through the formula (by far more compact than the sum of (2.8) and (LABEL:eq2-9)):

[TABLE]

Let us rewrite the variational wave function (3.2) as:

[TABLE]

Using the Funk-Hecke formula Estrada ,

[TABLE]

which is a remarkable integral identity holding for any sufficiently regular function $A$ of $x\in[-1,1]$, we first obtain $c_{lm}$ (up to the still unknown $C_{\alpha}$ constant) as:

[TABLE]

As for $d_{lm}$, one observes that

[TABLE]

with ${\bf v}_{k,k^{\prime}}=\hat{\bf R}_{k}+\hat{\bf R}_{k^{\prime}}\,,v_{k,k^{\prime}}=|{\bf v}_{k,k^{\prime}}|$, and $h_{k,k^{\prime}}(\hat{\bf r}\cdot\hat{\bf v}_{k,k^{\prime}};\alpha)=\exp\{-D\alpha\left[2-v_{k,k^{\prime}}(\hat{\bf r}\cdot\hat{\bf v}_{k,k^{\prime}})\right]\}$. In the double sum above, the contribution from those pairs of vertices that are diametrically opposite on the sphere has been taken into account separately. Upon inserting the above expression into the second of Eqs. (A.16), we get:

[TABLE]

with $h_{k,k^{\prime},l}(\alpha)=2\pi\int_{-1}^{1}{\rm d}x\,h_{k,k^{\prime}}(x;\alpha)P_{l}(x)$. Finally, we obtain $C_{\alpha}$ by imposing $\psi$ normalization:

[TABLE]

In the above equation, the integral of $h_{k,k^{\prime}}$ over the full solid angle was evaluated in a coordinate system where $\hat{\bf v}_{k,k^{\prime}}$ is aligned with the $z$ axis (moreover, notice that $C_{0}^{-2}=4\pi n^{2}$).

Equations (4.6), (4.7), and (4.9) allow to determine $c_{lm}$ and $d_{lm}$ exactly for all $l$ and $m$, and then the energy ${\cal E}$ from Eq. (4.3). Apparently, this route to ${\cal E}$ has always to be preferred to MC integration. In practice, a limitation comes from the rate of convergence of the two series (A.16), which is slower the larger $R$ and $\rho$, and this entails computing a lot of Fourier coefficients (the calculation of Legendre polynomials and spherical harmonics for large $l$ would not be a problem, as it can be carried out to any desired precision using the recurrence relations obeyed by these functions).

For the sake of clarity, let us consider the icosahedral case. For $R=1.4$ and $\rho=14$, we have computed the Fourier coefficients of $\psi$ and $\psi^{2}$ up to $l=16$. Giving the icosahedron an orientation such that two of its vertices lie on the $z$ axis, the only non-zero coefficients are those for $l=0,6,10,12,16$ and $m=0,\pm 5,\pm 10,\pm 15$, and are all real (a different orientation would imply different coefficients, but the weights $\sum_{m}|c_{lm}|^{2}$ and $\sum_{m}|d_{lm}|^{2}$ of each $l$ sector will be invariant note ). By truncating both series in Eq. (4.3) at $l_{\rm max}=6,10,12$ (which, we stress, is not equivalent to truncating the $\psi$ series (2.7) at $l=l_{\rm max}$), one obtains the energy plots in Fig. 5. We see that the MC data are already well reproduced with $l_{\rm max}=12$ (up to $\alpha\approx 20$), while a smaller $l_{\rm max}$ is insufficient to obtain good results unless $\alpha$ is low (indeed, with the exception of $l=m=0$, $|c_{lm}|$ and $|d_{lm}|$ are all increasing functions of $\alpha$).

Let us finally present our variational MF results for the ground-state diagram of spherical PSM bosons, which has been constructed by only considering the possibility of inhomogeneous phases with the symmetries of the polyhedra listed in Table I. Looking first at the fluid-solid transition lines in the $R$-$\mu$ plane (Fig. 6), we see that each particular cluster phase can only exist in a finite range of $R$ values. In all cases, the transition occurs for a $\mu$ value larger than on the infinite plane (meaning that solid-like order is discouraged by the curvature of the sphere, as expected). We find a difference in behavior between those cases (tetrahedron, cube, and dodecahedron) where each cluster has only three other clusters around, and the other phases with a “cluster coordination number” $Z$ larger than 3. While for the latter phases the transition is invariably first-order, it is of mixed type for $Z=3$, i.e., continuous for small radii and first-order otherwise, with both characters coexisting in a narrow interval of radii. For a $R$ in this range, on increasing $\mu$ the fluid first freezes continuously; then, a secondary minimum develops in $\Delta{\cal E}(\alpha)$, until the system eventually undergoes a second, now discontinuous, isostructural transition (see an example in Fig. 7 and another one in Fig. 8). This scenario is only in part reminiscent of the behavior in flat space, where the transition is always continuous for open lattices, while being first-order for the compact ones Prestipino3 .

In the same Fig. 6 we have highlighted in pink the region of $R$-$\mu$ plane where the fluid is mechanically unstable. In flat space the loss of fluid stability above a certain density is heralded by the softening of roton-like excitations, which signals an instability towards the formation of a solid-like density wave. The same happens on a sphere, and we discuss at length in Appendix D how this phenomenon precisely occurs as a function of $R$. A remarkable finding is that, similarly as on a plane Prestipino3 , continuous freezing falls exactly at the upper stability threshold of the fluid. For example, solid-like fluctuations with $l=3$ become costless right at the continuous transition to a tetrahedral phase, whose lowest non-zero modes beyond $l=0$ are indeed $l=3$ and $m=0,\pm 3$ (when one vertex of the tetrahedron lies at the north pole of the sphere). A further message from Fig. 6 is the existence of $R$ intervals where, on increasing the density, the fluid becomes unstable before freezing. This is clearly impossible, and the reason why this occurs is that we have actually missed to identify all the relevant phases of the system — since, probably, the underlying polyhedra have non-equivalent vertices, i.e., they have low symmetry or no regularity at all (suggestions on where to search may come from numerical studies of the Thomson problem Wales ; Neubauer ).

The freezing and melting lines of spherical PSM bosons at $T=0$ are shown in Fig. 9. Here we can appreciate the difference in “transition strength” between the various phases, which not for nothing is higher for the polyhedra having large faces (i.e., for the solids whose vertices are less efficiently spread over the surface). For completeness, we report in Fig. 10 (left panel) the position, denoted $\alpha_{\rm min}$, of the negative minimum in $\Delta{\cal E}(\alpha)$. We see that $\alpha_{\rm min}$ typically increases with $R$ and is larger the less stable the cluster phase (but there are anyway exceptions). Instead, in the right panel of Fig. 10 we show the value at melting of the mean cluster size $N_{\rm cl}$ as a function of $R$. We see that, in the $R$ range where each cluster phase is maximally stable, $N_{\rm cl}$ lies between 20 and 30 for most of the phases, i.e., near the value ($25.83$) characteristic of the triangular cluster crystal Prestipino3 .

Finally, the full phase diagram of the system at $T=0$ is presented in Fig. 11. It includes as many as seven cluster phases in the $R$ interval from roughly $0.5\sigma$ to $2.5\sigma$. However, as commented before, this list of thermodynamically stable phases is far from exhaustive. Even in the quoted range of $R$, we can safely say that not all ground states of the system have been identified; on the other hand, the analysis made already allows to draw some conclusions: a) Each cluster phase spans a certain interval of $R$; if in a given range of radii there are many phases competing for stability, the winner is the one providing the most efficient occupation of the surface, or, equivalently, the highest cluster coordination number $Z$. b) As the sphere radius increases, the particles find it convenient to re-organize, adjusting the number of clusters on the surface so as to keep the distance between neighboring droplets near the magic value of $1.51\sigma$ (for the PSM); in turn, this implies an increase of $Z$ with $R$ towards the asymptotic value of 6. c) While the number of clusters is determined by $R/\sigma$, it is nevertheless nearly independent of the density (indeed, the solid-solid loci in Fig. 11 are almost vertical).

V Supersolidity of the spherical cluster phases

We conclude our analysis by showing that the cluster phases identified in the previous Section are all supersolid. In the supersolid phase of matter, still elusive in 4He but found in numerous lattice models Batrouni ; Wessel ; Pollet and, eventually, also observed in a quantum system with continuous ground-state degeneracy Leonard , the periodic density modulation typical of a solid coexists with the dissipationless flow of a superfluid.

As proposed by Leggett Leggett2 , a supersolid can be characterized by its response to uniform axial rotations: under a slow rotation a fraction of the quantum solid may stand still, with the result that its moment of inertia is smaller than expected from classical mechanics. Leggett has called superfluid fraction of a quantum solid the quantity

[TABLE]

where $I$ is the moment of inertia around the axis of rotation and $I_{0}$ its classical value. Supersolidity occurs when $f_{s}>0$.

When a system of $N$ particles is subject to rotation, say, around the $z$ axis, the infinitesimal change in energy due to rotation is $\omega{\rm d}L_{z}$, where $\omega$ is the angular velocity and $L_{z}={\cal O}(N)$ the $z$-component of the total angular momentum. For bosons in the product state (2.2), the average $L_{z}$ per unit particle is given by the symmetrized expression

[TABLE]

(where we are assuming that the sphere is immersed in 3D space). At $T=0$ we require that the state $\psi$ minimizes the free-energy functional ${\cal E}[\psi]-\omega\left\langle\Psi|L_{z}|\Psi\right\rangle/N$ Leggett3 , which for small $\omega$ is the specific energy in the absence of rotation minus $(I/N)\omega^{2}/2$ Aftalion . In other words:

[TABLE]

When $\omega$ is non-zero the quantum state $\psi$ acquires a phase, $\Theta(\Omega)=\omega S(\Omega)+{\cal O}(\omega^{2})$, whereas the square amplitude changes from $\eta_{0}$ (viz., the square amplitude for $\omega=0$) to $\eta=\eta_{0}+\omega\eta_{1}+{\cal O}(\omega^{2})$. Inserting Eq. (C.44) into (5.2), we eventually obtain:

[TABLE]

leading in turn, up to terms of order $\omega^{3}$, to

[TABLE]

${\cal E}_{0}[\eta]$ being the energy functional for $\Theta={\rm const}$. Therefore:

[TABLE]

Upon considering that

[TABLE]

we finally obtain:

[TABLE]

While the calculation of $f_{s}$ is difficult, finding a lower value is easier:

[TABLE]

where $\eta_{0,{\rm min}}$ and $\eta_{0,{\rm max}}$ are the minimum and maximum values of $\eta_{0}$ on the sphere. Hence, a finite density contrast $\eta_{0,{\rm max}}/\eta_{0,{\rm min}}$ is the fingerprint of supersolidity. To obtain the estimate in Eq. (5.9), we have considered that the minimum of $\int{\rm d}^{2}\Omega\left[\nabla S-(m/\hbar)\hat{\bf z}\wedge{\bf r}\right]^{2}$ is reached for $\nabla S=0$ (indeed, the Euler-Lagrange equation for the functional in (5.9) is $\nabla^{2}S=0$, because $\nabla\cdot(\hat{\bf z}\wedge{\bf r})=0$, and the only regular $S(\Omega)$ with $\nabla^{2}S=0$ is proportional to $Y_{0}^{0}$, hence it is a constant). Since $\eta_{0}$ is positive definite for our variational solution, we conclude that $f_{s}$ is strictly positive — in other words, in our theory all cluster phases are supersolid for every $\mu$. It cannot be excluded that, for very large $\mu$, a $\psi$ function vanishing in the interstitial region between the polyhedron vertices can have a lower free energy than the Gaussian ansatz. Even in this case, supersolidity will occur at least in the vicinity of the freezing point.

VI Conclusions

An efficient method to study weakly-interacting bosonic particles at zero temperature is mean-field theory, which assumes a perfect condensate for the system state. As a further simplification, the single-particle wave function can be accurately modeled through some physically-motivated ansatz Prestipino3 ; Prestipino4 , which is then optimized by use of the variational method. In the present study, we have employed variational mean-field theory to investigate a (finite) system of penetrable bosons confined to a spherical surface, essentially with the aim to follow the evolution of ordering tendencies with the radius $R$ in a genuinely quantum system.

The ground-state diagram of the system is very rich, featuring many distinct high-density “phases” as a function of $R$, all characterized by the presence of clusters of overlapping particles. The mechanism behind the onset of cluster phases on a sphere is the same as in flat space, and is purely energetic in character. We have found that many, but definitely not every, of these phases have clusters distributed at the vertices of a regular or semi-regular polyhedron inscribed in the sphere, and the stable phase at a given $R$ is the one ensuring the maximum possible (cluster) coordination number that is consistent with a distance between neighboring clusters of about $1.51\sigma$ (for PSM bosons Prestipino3 ), i.e., roughly the same as in the triangular cluster crystal. The existence of cluster phases is intimately rooted in the characteristics of the interaction between particles, and can be anticipated from the nature of the elementary excitations of the (super)fluid phase, which are roton-like for sufficiently large densities. When the roton mode eventually softens, the fluid becomes unstable towards a solid-like density modulation, and that marks the upper boundary of the homogeneous phase. In practice, unless the cluster coordination number is very small, the phase transition occurs before reaching the instability threshold, and in this case freezing is first-order.

Finally, we have given an analytic argument showing that, at least according to our variational analysis, all cluster phases are supersolid, i.e., they exhibit a reduced moment of inertia compared to its classical value. We ascribe this property to the finite strength of interparticle forces, which, by allowing particles to diffuse freely within the surface, can sustain a superfluid component in the cluster phases.

Our results can find application for the behavior of ultracold gases of bosonic atoms confined in spherically-symmetric bubble traps Zobay ; Garraway , as will be made available in future experiments carried out in a microgravity environment Elliott .

Acknowledgements

One of us (S. P.) wishes to thank A. Sergi and E. Bruno for stimulating discussions at an early stage of this research.

Appendix A Gross-Pitaevskii equation on a sphere

In this Appendix we give a variational derivation of the GP equation, different from the one provided in Rogel-Salazar and adapted to the sphere case.

The starting point is the MF energy functional, written for a general (i.e., normalizable but not necessarily of unit norm) single-particle wave function $\psi$:

[TABLE]

with

[TABLE]

and

[TABLE]

In Eqs. (A.11) and (A.12), the solid-angle element is ${\rm d}^{2}\Omega=\sin\theta\,{\rm d}\theta{\rm d}\phi$. According to the variational principle, the “best” approximate ground state is such that the average energy $\langle H\rangle$ in (A.10) is minimum. For the latter to occur, a necessary condition is $\delta\langle H\rangle=0$. Upon observing that

[TABLE]

we readily arrive at:

[TABLE]

with

[TABLE]

It is hardly necessary to observe that, when the radius $R$ is of order $\sigma$, $N-1$ may not be quite the same as $N$. Now observe that, if $\psi$ obeys (A.14) also $c\psi$ is a solution, for all $c\neq 0$. In particular, we can always choose $c$ such that $\psi$ is normalized to 1. In this case, Eq. (A.14) becomes the GP equation in its standard form (see Eqs. (2.3) and (2.5)), with $\lambda=(\langle K\rangle+2\langle U\rangle)/N$.

The most natural way to solve the GP equation is to expand the solution and its square modulus in a series of spherical harmonics,

[TABLE]

where the $d_{lm}$ coefficients can clearly be expressed in terms of the $c_{lm}$ themselves (see Eq. (A.26) below). Observe that the $d_{lm}$ do not fulfil any particular sum rule (because, at variance with $\psi$, the function $|\psi|^{2}$ is not subject to any specific normalization). Moreover, since $|\psi(\Omega)|^{2}$ is real, it is generally $d_{l,-m}=(-1)^{m}d_{lm}^{*}$. The rest of the derivation follows the one provided in Appendix B for the MF energy functional, and we finally arrive at:

[TABLE]

with

[TABLE]

Equation (A.21) is a non-linear set of equations, similar to that reported in Eq. (2.6) of Ref. Prestipino3 . Therefore, it is tempting to solve it self-consistently, assuming arbitrary initial values for the $c_{lm}$ and filtering out the minimum energy eigenvector of unit norm at each iteration step. However, while this method works well in the planar case it dramatically fails to converge in the present case (unless the density is very small), and the problem is not alleviated by the use of a mixing scheme. To guarantee that iteration of (A.21) becomes a contractive fixed point iteration, we may think to replace spherical harmonics with a different basis of functions, possibly a different one for any specific radius and polyhedral symmetry. However, besides the difficulty of devising specific basis functions for each case, the price to pay is losing the good properties of spherical harmonics that allow simplifying the final form of the GP equation. In view of this, a more viable procedure is to resort to a variational ansatz for $\psi$, which also offers the advantage of better elucidating the physics behind the minimum-energy state.

Appendix B Proof of Eq. (LABEL:eq2-9)

Here we obtain an expression for the potential-energy term in Eq. (2.6). We first observe that any bounded $u$ with finite support can be written as a Fourier integral:

[TABLE]

(notice that $k$ is a dimensionless variable). Upon considering that

[TABLE]

and using the expansion of a plane wave in spherical harmonics, we get:

[TABLE]

$j_{l}(k)$ being a spherical Bessel function. Alternatively, and also more conveniently, $u(x)$ can be written as a series of Legendre polynomials:

[TABLE]

By noting that

[TABLE]

we promptly obtain:

[TABLE]

In the following we prefer using (B.32) in Eq. (2.6), rather than (B.29), while deferring a direct proof of

[TABLE]

to below in this Appendix.

Next, we note that

[TABLE]

Putting all things together:

[TABLE]

The integral of a product of three spherical harmonics (also called Gaunt coefficient) has a known value:

[TABLE]

where the $2\times 3$ tables are Wigner 3-j symbols Stone . As a result, Eq. (B.35) gets simplified, eventually transforming into Eq. (LABEL:eq2-9).

Finally, we provide a proof of Eq. (B.33), which establishes the equivalence between (B.29) and (B.32). We seek for a different expression of

[TABLE]

Taking advantage of the formula Ludu

[TABLE]

which is valid for all $l$ and $k$, we first obtain (for any $x$ satisfying $-1<x<1$):

[TABLE]

where in the last step we used the property $P_{l}(-x)=(-1)^{l}P_{l}(x)$. Similarly,

[TABLE]

The case of $x=\pm 1$ needs a different treatment, since the delta argument vanishes at one of the extrema of the $t$ integration interval. Indeed, one can prove that

[TABLE]

(rather than $\pi(\pm i)^{l}$, as it would follow from Eq. (B.39) and (B.40)).

Coming to the calculation of (B.37), we have:

[TABLE]

where, in consideration of $j_{l}(-k)=(-1)^{l}j_{l}(k)$, the last step holds for both even and odd $l$. The inner integral in (LABEL:b-16) can be evaluated for any $-1<x<1$ using Eqs. (B.39) and (B.40):

[TABLE]

Despite the case $x=\pm 1$ would need a separate treatment, the value of the integrand at the $x$ boundary in (LABEL:b-16) is irrelevant for the value of the same integral. By plugging Eq. (B.43) in (LABEL:b-16), we readily arrive at the desired Eq. (B.33).

Appendix C The condensate of minimum energy is real

In this Appendix, we show that the true ground state of a system of spherical bosons with specific energy given as in Eq. (2.6) is necessarily represented by a real wave function.

Using the Madelung representation,

[TABLE]

with $\eta\geq 0$ and a real $\Theta$ defined up to an arbitrary additive constant, the potential-energy term in Eq. (2.6) immediately reads (pulling out the homogeneous-system energy):

[TABLE]

Aside from a constant factor, the kinetic-energy term is given by

[TABLE]

where the last step follows after partial integration. Upon inserting (C.44) into (C.46), and taking account of the expression

[TABLE]

of the gradient of a scalar function $f(\theta,\phi)$, we obtain:

[TABLE]

whence finally:

[TABLE]

Equation (C.49) represents, mutatis mutandis (i.e., with $4\pi$ in place of the volume), the same result holding in flat space. Looking at Eqs. (C.45) and (C.49), it is clear that the minimum (kinetic) energy is attained for a constant $\Theta$. As a global phase in the wave function cannot affect the results, we can always choose $\Theta=0$. This implies that the MF ground state is real and non-negative.

Appendix D Elementary excitations of the fluid

We hereafter investigate the collective excitations of the system in the fluid phase, using a method similar to that employed in Ref. Prestipino3 . The starting point is the time-dependent GP equation,

[TABLE]

describing the MF quantum dynamics of identical bosons at $T=0$. Multiplying Eq. (D.50) by $\psi^{*}$ and subtracting the complex conjugate of the resulting equation, we arrive at the following continuity equation:

[TABLE]

where we have considered that the divergence of a vector field ${\bf A}=A_{\theta}\hat{\bm{\theta}}+A_{\phi}\hat{\bm{\phi}}$ is

[TABLE]

In (D.51), the velocity field is clearly

[TABLE]

where the last step follows after inserting the Madelung form (C.44).

Another equation connecting $\eta$ and $\Theta$ is obtained by plugging Eq. (C.44) in the time-dependent GP equation, with the result that:

[TABLE]

While the imaginary part of (D.54) gives back the continuity equation, the real part reads:

[TABLE]

Taking the gradient of (D.55), we arrive at a Navier-Stokes-like equation without viscosity term:

[TABLE]

where it is intended that

[TABLE]

We now derive an approximate equation valid for a $\psi$ function departing only slightly from the homogeneous-system solution, $\eta=1$ and $\nabla\Theta=0$. Such perturbed solutions are the elementary excitations of the fluid phase. Inserting $\eta=1+\delta\eta$ and $\nabla\Theta=\delta{\bf u}$ into the continuity equation (D.51), and simply ignoring every term that is not linear in $\delta\eta$ or $\delta{\bf u}$, we first get:

[TABLE]

Moreover, we have:

[TABLE]

which eventually allow us to simplify Eq. (D.56) as:

[TABLE]

Inserting Eq. (D.60) into the second of Eqs. (D.58), we finally arrive at:

[TABLE]

This equation admits solutions in the form $\delta\eta=\varepsilon{\rm Re}\left\{Y_{l}^{m}(\Omega)e^{i\omega_{l}t}\right\}$, where $\varepsilon$ is a small dimensionless amplitude. The dispersion relation of these waves can be obtained by observing that:

[TABLE]

In particular, we derived the last equation above from the Funk-Hecke formula (4.5) for $A=u$. Substituting Eqs. (D.62) into (D.61), we finally obtain:

[TABLE]

which can be viewed as the spherical version of the Bogoliubov spectrum.

As long as the r.h.s. of (D.63) is positive, the fluid phase is stable and (by the same argument exposed in Sec. V) superfluid; conversely, if $\omega_{l}$ is purely imaginary, quantum dynamics will drive the system arbitrarily far from $\eta=1$. It turns out that, for each fixed $l$, $\omega_{l}^{2}$ turns from positive to negative at a density of

[TABLE]

Hence, when the density exceeds a certain value the fluid becomes destabilized. The upper threshold for fluid stability in terms of chemical potential is

[TABLE]

It turns out that $\mu_{\rm inst}(R)$ shows cusps where the $l$ value providing the minimum in (D.65) jumps by one (see Fig. 12). Quite remarkably, the oscillatory behavior of $\mu_{\rm inst}(R)$ is similar to that computed within density-functional theory for the $\lambda$-line of a classical fluid of spherical soft-core particles Franzini . In Fig. 12 the $\mu_{\rm inst}(R)$ locus is reported for three distinct models of interaction, namely the PSM, GEM-8, and GEM-4 potentials. We see that the region of fluid stability extends more and more the smoother is the generalized-exponential interaction (same as found in flat space Prestipino3 ). As $R$ grows to infinity, the instability line flattens out until it finally equates the planar threshold (e.g., $46.2979\ldots$ for the PSM potential Prestipino3 ).