A variational mean-ﬁeld study of clusterization in a zero-temperature system of soft-core bosons

. We work out the ground-state diagram of weakly-repulsive penetrable bosons, using mean-ﬁeld theory with a Gaussian ansatz on the single-particle wave function. Upon compression, the ﬂuid transforms into a cluster supersolid, whose structure is characterized for various choices of the embedding space. In Euclidean space, the stable crystals are those with the most compact structure, i.e., triangular and fcc in two and three dimensions, respectively. For particles conﬁned in a spherical surface, as the sphere radius increases we observe a sequence of transitions between di ﬀ erent cluster phases, all having a regular or semiregular polyhedron as supporting frame for the clusters. The present results are relevant for the behavior of ultracold bosons weakly coupled to a Rydberg state.


Introduction
The ground state of a collection of N non-interacting identical bosons is a product of single-particle states, all equal to the ground state of one particle only. We say that the system state is a perfect condensate. The condensate is a phase-coherent state, since the wave function of an individual particle is in phase with the wave function of other particles. Upon heating, phase coherence is reduced and eventually destroyed by thermal fluctuations, as first predicted by Einstein for an ideal gas of bosons almost a century ago.
The concept of condensate can be extended to interacting many-body systems, if the highest eigenvalue of the one-body density matrix is a macroscopic number; then, there is a single-particle wave function that is occupied in a macroscopic way (a property equivalent to off-diagonal long-range order) [1]. For example, in simple models of bosons on a lattice, when the tunneling/hopping term dominates the interaction term in the Hamiltonian, then the system has a ground state consisting of every particle in the same state spread across the entire lattice.
We hereby focus on the behavior at zero temperature (T = 0), where many condensed quantum systems may still undergo (quantum) phase transitions driven by a non-thermal control parameter, such as pressure, magnetic field, or chemical composition. A simple example is the freezing of an ultracold fluid of soft-core bosons, as could be realized in bosonic atoms dressed with Rydberg states [2,3]. In purely classical terms, an interaction that is everywhere finite can stabilize cluster crystals at low temperature and high density (see, e.g., [4]). Indeed, it may be energetically more convenient to form isolated blobs of * Corresponding author: sprestipino@unime.it particles than having them distributed homogeneously in space [5,6]. Clusterization also occurs in quantum systems, as first shown theoretically by Pomeau and Rica in 1994 [7].
In the weak-interaction limit, an effective approach to the physics of ultracold bosons is mean-field theory (MF), which assumes that the many-body state is a perfect condensate. The "best" single-particle state ψ (that is, the one of minimum energy) obeys the Gross-Pitaevskii (GP) equation [1], where u is the pair potential and λ is the Lagrange multiplier enforcing the normalization condition Kunimi and Kato [8] have solved the GP equation in two dimensions, showing that the high-density phase is a triangular cluster crystal. Macrì and coworkers [9] have tested MF results against Monte Carlo simulation of a weaklyinteracting PSM system, finding that the exact freezing point lies indeed extremely close to the theoretical estimate. However, solving the GP equation for a threedimensional (3D) crystal is not as easy as in two dimensions, since the computational effort is much greater. This is a severe problem when many different crystalline structures compete for the ground state. A cheaper strategy, which we pursue here, is to take a reasonable form of con-EPJ Web of Conferences 230, 00008 (2020) https://doi.org/10.1051/epjconf/202023000008 FisMat 2019 densate wave function, dependent on parameters, and optimize its shape by minimizing the MF energy functional. We have checked in a few two-dimensional (2D) cases that the best variational state reproduces the GP wave function and energy very accurately [10]. Then, by examining a wide spectrum of 3D lattices, we have identified stable and metastable solid phases and characterized their melting behavior.

Three-and two-dimensional systems
The MF energy functional reads: representing the sum of zero-point kinetic energy and potential energy per unit particle. When the density is low, energy is minimized by spreading the wave function ψ over the system volume, and the stable phase is then fluid. On the contrary, for very high densities the energy is lower for a strongly inhomogeneous ψ peaked at the sites of a crystalline lattice (at least provided that the bounded potential u is "fatter" than Gaussian). Within the variational method, we write the condensate wave function as a sum of Gaussians: where the R's are direct-lattice vectors, while α and a are adjustable parameters. Clearly, α = 0 corresponds to the fluid phase. The normalization constant in (4) is given by [10] where v 0 is the volume of the primitive cell. To allow for the possibility of cluster crystals, the lattice constant a and the number density ρ are taken as independent variables. Denoting N c the number of lattice cells, the number of particles per cell is on average If in equilibrium ρv 0 > 1, then the solid phase is a cluster crystal.
On the other hand, ψ can also be written as a Fourier series, where the G's are reciprocal-lattice vectors. The Fourier coefficients of ψ are given by [10] ψ Unless α is exceptionally large, the Fourier series is better suited for the computations than the original ψ expression.
The main advantage of the Gaussian ansatz is an analytic simplification of the energy functional, allowing a considerable speed up in the numerical calculations [10]. Indeed, the kinetic energy admits a closed-form expression: where e 0 = 2 /(mσ 2 ) is a natural energy unit. At variance with E kin , the formula of the potential energy is explicitly lattice-dependent. For the triangular lattice, we find: where u(k) is the real-valued Fourier transform of u. In the first sum, only even-even combinations of reciprocallattice basis vectors b 1 and b 2 do occur; the other G's are exclusively involved in the second sum. Finally, .
(10) An expression similar to (9) holds for any Bravais or non-Bravais lattice, both in two and in three dimensions. Choosing the pressure P as the control parameter, the most stable phase at T = 0 is the one with the lowest enthalpy.
We confirm that, in two dimensions, the fluid freezes into a triangular cluster crystal (the coexistence densities are ρ f = 12.315 and ρ tr = 13.131, with ρ tr v 0 = 25.8 [10]). In three dimensions, the only relevant zero-temperature PSM phases besides the fluid are the compact cubic crystals (fcc, bcc, and hcp), see Fig. 1. Loosely packed crystals, such as the simple-cubic crystal and the diamond crystal, melt continuously at a common pressure P c that The transition pressure for the honeycomb and striped crystals is the same as for the square crystal [10]. Right: 3D case. The arrows mark the transition into the fcc crystal (left), the sh crystal (middle), and the sc crystal (right). The transition pressure for the diamond crystal is the same as for the sc crystal [10]. is much larger than, say, the fcc melting pressure. The fluid freezes into a fcc crystal (ρ f σ 3 /e 0 = 14.294 and ρ fcc σ 3 /e 0 = 16.599 are the density values at coexistence, with ρ fcc v 0 = 36.7 [10]), but the difference in enthalpy between fcc and hcp is very small. As the pressure increases, the hcp crystal eventually prevails, and a fcc-to-hcp transition then takes place.
The MF spectrum of fluid excitations can be obtained by solving, in the amplitude-phase representation, the time-dependent GP equation for a condensate wave function that deviates only slightly from homogeneity. It turns out (see e.g. [10]) that the perturbation has the form of plane waves with a dispersion relation of Bogoliubov type, In particular, if u(k) is negative in a range of k values then the fluid phase exhibits superfluid behavior (by Landau's argument). For the 3D PSM fluid, ω(k) shows a roton minimum for high enough densities (see Fig. 2). The roton minimum becomes unstable exactly at the transition density for continuous freezing, which will thus coincide with the upper stability threshold of the fluid. A further outcome of our theory is the supersolid character of the crystalline phases. A supersolid phase accommodates both diagonal and off-diagonal long-range order, meaning that particles are arranged in a crystalline structure but, at the same time, can flow without dissipation. Like a superfluid, also a supersolid can be characterized by the nature of its response to uniform axial rotations. Under a slow rotation, a fraction of the supersolid stands still, as witnessed by the value of its moment of inertia I being smaller than the value I 0 expected classically [11]. We have shown [10] that in MF theory the superfluid fraction 1−I/I 0 has a characteristic lower bound, ψ 2 min /ψ 2 max , where ψ min and ψ max are the minimum and maximum value of ψ in a cell. In our variational theory this lower bound is strictly positive, implying that the crystal is supersolid at every pressure. In fact, above a certain pressure the value of ψ in the interstitial region is so small that phase coherence will be destroyed by quantum fluctuations, and the system then enters a "normal" solid phase (see, e.g., [12]).

One-dimensional system
In the one-dimensional (1D) case, which we treated in Ref. [13], the analytic simplifications introduced by the Gaussian ansatz are even stronger than in higher dimensions. The only problem with MF theory is that a rigorous result by Pitaevskii and Stringari rules out, at zero temperature, any form of long-range order (including offdiagonal long-range order) in 1D bosonic systems with continuous group symmetries [14]. However, if we add a trapping potential in the axial direction, both the amplitude and phase fluctuations are suppressed at low T and one has a true condensate in one dimension, at least provided that N is large enough [15]. Another issue concerns the range of validity of the MF approximation. A hand-waving argument goes as follows: MF theory is expected to hold when the healing length / mρ u(0) (with u(0) ≈ σ d ), fixing the length scale above which collective physics dominates over singleparticle physics, is much larger than the average interparticle separation ρ −1/d . In equivalent terms, one might say that MF theory is valid whenever the kinetic energy per unit particle due to localization, 2 ρ 2/d /(2m), is much larger than the interaction energy u(0)ρ. In three dimensions, this leads to ρσ 3 /e 0 (e 0 / ) 2 , which corresponds to a density range that is wider the weaker the interaction strength. However, in one dimension the MF regime is rather ρσ /e 0 ( /e 0 ) 2 , and the approximation improves with increasing density.

Soft-core bosons on a sphere
Finally, we have extended our theory to the case of particles embedded in a spherical surface. The systems which more closely conform to this setting are ultracold bosons trapped in a thin spherical shell, which will be the subject of upcoming experiments under microgravity conditions [16].
Spherical boundary conditions are often employed in simulations as a means to discourage crystalline ordering at high density (as is well known, triangular order is frustrated on a sphere [17][18][19][20]). In practice, the curvature imposes a distinct excess of fivefold disclinations over sevenfold ones, which considerably complicates the search for optimal packings, even for small radii. Franzini and coworkers [21] have studied by density-functional theory a classical system of soft-core particles on a sphere, finding a rich catalog of cluster (pseudo)phases as a function of the sphere radius R (the suitability of the thermodynamic formalism is justified a posteriori by the large size of the system near the transition point). We expect the same to occur in a system of softly-repulsive bosons, which will also clusterize at high density, becoming solid-like inhomogeneous. As for the classical system, numerous structures would compete for stability, and the natural candidates for the dense phases are highly symmetric packings of clusters. Considering that for the 2D PSM the edge of the triangular lattice is about 1.51σ at the melting point, it is likely that the supporting frame of the clusters for a given R will be found among the regular or quasi-regular polyhedra with edge ≈ 1.51σ. For instance, since the radius of the sphere circumscribing the regular icosahedron equals an icosahedral cluster phase is most likely to occur for R 1.44σ. Expanding the condensate wave function in spherical harmonics, the specific energy E is again written as the sum of two terms [22], and where P l (x) are Legendre polynomials and d lm are the Fourier coefficients of |ψ| 2 . In practice, the l sum must be truncated, i.e., l ≤ l max , where l max is chosen in accordance with the spatial resolution adopted for the description [23]. We have first provided an analytic proof of clusterization for R = 1.45σ, i.e., where we expect that the symmetry of the stable phase is icosahedral. We employ a variational wave function that is a weighted sum of the two simplest icosahedral harmonics [24], with β a parameter to optimize. For this ψ the calculation of (14) and (15) is an analytic tour de force [22], and a graphical demonstration of the existence of the transition for the PSM fluid is provided in Fig. 3. Clearly, choosing the icosahedron as scaffolding for the clusters is just one possibility; as R increasingly departs from 1.45σ, other polyhedra will be better suited than the icosahedron for matching the condition ≈ 1.51σ. We have considered ten circumscribable polyhedra (see Table 1), and for each of them computed the grand potential as a function of R (for fixed radius and chemical potential, the stable phase is the one with the lowest grand potential). A one-parameter ψ adequate to represent the high-density patterns is a sum of Gaussians centered at the vertices R k (k = 1, . . . , n) of the given polyhedron: The fluid ground state (ψ = 1/ √ 4π) is recovered as a limiting case, that is for α = 0. We have followed two different but equivalent routes to compute the grand potential [22]: either i) a direct Monte Carlo estimate of the integrals originally defining E kin and E pot ; or ii) the calculation of a large number of Fourier coefficients for ψ and ψ 2 .
The resulting ground-state diagram of the PSM fluid, plotted in Fig. 4, includes as many as seven supersolid cluster phases in the R range from 0.6σ to 2.5σ. All transition lines are first-order, except for a portion of the melting line for the tetrahedral phase (marked as dashed in Fig. 4) which is second-order. Each cluster phase spans a certain  interval of R; when more phases compete for stability, the phase with the lowest value of grand potential turns out to be the one providing the most efficient occupation of the surface, or, equivalently, the highest coordination number.

Conclusions
We have used the variational method to evaluate, within MF theory and a Gaussian ansatz for the condensate wave function, the average quantum energy of a collection of fully penetrable bosons in spaces of various dimensionality. Using this information, we have identified the zerotemperature phases of the PSM system and clarified the nature of the pressure-induced freezing transition. In the dense phase particles are invariably gathered in clusters centered at the sites of a crystalline lattice. For the special case of particles embedded in a spherical surface, clusters instead fall at the vertices of a regular or semi-regular polyhedron inscribed in the sphere. Within our MF description, the dense phases are all supersolid, i.e., they are characterized by non-classical rotational inertia.