Ground state phase diagram of ultrasoft bosons

In 2D bosonic systems ultra-soft interactions develop an interesting phenomenology that ultimately leads to the appearance of supersolid phases in free space conditions. While suggested in early theoretical works and despite many further analytical efforts, the appearance of these exotic phases as well as the detailed shape of the ground-state phase diagrams have not been established yet. Here we develop a variational mean-field calculation for a generic quantum system with cluster-forming interactions. We show that by including the restriction of a fixed integer number of particles per cluster the ground-state phase diagram can be obtained in great detail. The latter includes the determination of coexistence regimes of crystals of different occupancy as well as crystals with super-fluid phases. To illustrate the application of the method we consider an ultrasoft potential for which the computational phase diagram is known. Our results show very good quantitative agreement with the simulations and suggest that the solid-superfluid coexistence could be a reliable marker to locate supersolidity.


I. INTRODUCTION
The theoretical study of the quantum phase transition occurring in the crystallization of soft-core bosons have received a particular attention over the last decade [1][2][3][4][5]. This interest is twofold. On the one hand ultrasoft bosons allow to understand exotic many-body effects in which quantum fluctuations weakens crystalline order to give rise to, e.g., supersolid states [6][7][8][9]. On the other hand, the advances in experimental Bose-Einstein condensation techniques involving cold atomic gases have boosted and fed back the need for theoretical insights on how to control these phases of matter [10][11][12][13].
In the supersolid phase, the systems still display a modulated density that breaks the rotational and translational symmetries, while at the same time develop superfluid transport properties. The characterization of the properties of these ground-state phase diagrams using specific interaction potentials have been a fruitful direction for numerical research [14][15][16]. Several analytical studies have focused as well on the supersolid and superfluid phases of systems with standard non-clusterforming soft interaction potentials. [17][18][19][20][21][22] Properties related to, e.g., the supersolid-superfluid transition [17,18] and the excitation spectra [21,22] are of great current interest for the scientific community. Many calculations use mean-field variants [20,21] and more sophisticated approaches [17,18] to tackle the problem in this scenario. In this context, the ground-state phase diagrams of ultra-soft bosons, in which the crystalline structures are formed by clusters instead of single particles, is a relatively new problem that is generating theoretical and experimental advances. [2,6,9] Bosons interacting via ultrasoft potentials crystallize in a lattice of clusters [23,24] which, at very small temperatures, melts due to the quantum fluctuations produced by the zero point motion. This quantum transition to a superfluid phase is always present in the ground state for the regime of small pressure and weak interaction strength. Near the melting point, however, a supersolid behavior have been predicted via numerical simulations and analytical approximations [2,25].
Recently, some important efforts have also been done to analytically determine these ground-state properties in the ultra-soft cluster-forming scenario, through formulations that can be applied to different potentials. In particular, by using generalizations of the Gross-Pitaevskii mean-field model [4] and variational formulations of the condensate wave function [3,9], a number of predictions have been put forward. These predictions include the order of the transitions and the dispersion relations of the supersolid lattice. Other important properties like the location of the phase boundaries are determined only approximately [3,4] and not in full agreement with the outcomes of Monte Carlo simulations.
In the present work we go a step further in the analytical understanding of the ground-state phase diagrams of ultrasoft bosons. By using a variational mean-field calculation we improve the predictions on the location of the phase boundaries of generic ultrasoft potentials. The main assumption used in the present formalism is the restriction of the number of particles within the clusters to be a constant integer value. We show that, for the model potential considered in Ref. [2], the ground-state phases and the coexistence regimes can be calculated to a degree of reliability that is consistent with Monte Carlo numerical simulations in moderate and high density. This interaction form is of relevance for ultracold atoms [18,26] and has been used to model a large variety of other softmatter systems [27][28][29].
We determine the solid-superfluid coexistence region (SSC) by looking at the chemical and mechanical equi-librium condition that exist between the pure solid and superfluid thermodynamical phases. Clearly, this regime can not be interpreted as a supersolid by any means. However, for the model potential used here as example, there is a coincidence regarding the shape and the location of the SSC region analytically calculated and the supersolid phase that have been observed via exact numerical simulations [2].

II. METHOD
Consider a system of interacting bosons in 2D with a pair potential V ( r) = U v( r). The units of energy and length are such that the kinetic energy operator of a particle is given byT = −( ∇) 2 /2. Under these conditions the Hamiltonian of the system is given bŷ Here the potential v(r) is bounded and its Fourier transform has a negative minimum at some finite wave vector. These conditions determine the ultrasoft character of the interactions and ensure the cluster-forming character of system in the in the classical regime. The characterization of the ground-state phase diagram is done by means of a variational approach. It is thus expected that our results will better describe the system in the regime where the interparticle distance is is small than the length scale of the potential. A variational ground state function is proposed consistent with the mean-field approximation [3,9] where the single particle wave function is Here A stands for the area of the system and the set of vectors k j are selected in such a way that the field φ 0 ( x) reproduce the symmetries of a triangular lattice of particles. While other structures could be stable for specific potentials, we focus here on triangular lattices as this is the most common structure observed in simulations of ultrasoft potentials [2,30,31]. Consequently k j = (p u 1 + q u 2 )k 0 , with p and q integers, u 1 = (0, 1) and u 2 = ( √ 3/2, −1/2). We set c 0 = 1 without loss of generality.
These Ansatz for the single-particle ground state function are usually used in the literature to model states in which particles form a triangular crystal or even triangular crystal of clusters [3,32]. A central constraint, that in the case of cluster-crystals has not been taken into account to the best of our knowledge, is the relation between k 0 , the lattice spacing a and the occupancy number n of the clusters in the crystal, which is a positive integer. To consider n as a positive integer does not rule out states with fractional average occupation, since these states can be seen as the result of the phase coexistence of solids with consecutive integer occupation number.
Considering our definition of φ 0 ( x), it can be shown that for this function to be periodic with a spatial period a, the main modulation wave vector should be k 0 = 4π/( √ 3a). At the same time the average density ρ is related to the lattice spacing a and the occupancy number n by the relation This relation implies that cluster-crystals with different occupancy number, at the same density, have different lattice spacing. As a consequence This relation is at the core of numerous properties of the cluster-forming particle system. The disregard of this constraint, in the assumption that n can be treated as a variational real parameter, explains absence of transitions between different cluster-crystal phases in meanfield approaches when describing the zero and finite temperature properties of classical and quantum particles systems [3,33].
In numerical simulations cluster-forming particle systems have been observed to remain with a fixed occupancy number while density is increased [34,35]. While this fact have been already observed, it is generally considered that the presence of hopping between neighboring clusters justifies the consideration that n and consequently k 0 are variational parameters [30,36,37].
In this work we assume that n is an integer value. This means that the cluster-crystal phases consist in lattices of equally occupied clusters of bosons. We use this assumption to make a complete analysis of the ground-state phase diagram indicating the regions corresponding to the different clusters phases. At some special regions, of course, coexistence between states with occupancy number n and n + 1 exists, but only in these regions the occupancy number can be interpreted as varying continuously between n and n + 1.
Proceeding with the construction of the energy functional, the energy per particle of the cluster crystal with n particles per cluster, at an average density ρ, is given by where once we take into account the rotation symmetry of the pair interaction potential. By minimizing the energy per particle Eq. (8) in terms of the coefficients c j the ground-state (U versus ρ) phase diagram of the system can be obtained in this mean-field approximation.

III. RESULTS
For the quantum regime it has to be considered not only the different cluster configurations, characterized by its occupancy number, but also the superfluid homogeneous phase. This homogeneous configuration can be recovered from the modulated solution φ 0 ( x) by setting c j = 0 for all j = 0. So it is straightforward that the energy per particle of the homogeneous phase is given by To minimize analytically the function of Eq. (8) considering the whole set of coefficients c j can be a very difficult task. Consequently, finding closed expressions for the boundaries between the different phases is fairly complicated. However, the more interesting transition happening in the quantum regime, i.e. the quantum melting of the solid phases into the superfluid phase, is in fact tractable within certain approximations. Around this transition, quantum fluctuations characterized by the amplitude of the zero point motion are expected to be strong enough to destabilize the crystalline order. In this scenario the particles are less localized and the profile is smoother, making the Fourier amplitudes of the profile smaller for all non-zero wave vectors. In order to analytically advance we consider the single mode approximation, c j = 0 for j > 1. Such as an approximation can not be fully justified with qualitative arguments, and its ultimate validity can be verified by comparing with general mean-field numerical results and with previous Monte Carlo simulations.
A. Analytical results on the quantum melting Following the single mode approximation c j>1 = 0 lets proceed with the analytical study of the quantum melting of the solid phase. In the single mode approximation the energy per particle of the system (Eq.8) yields where k 0 (n, ρ) is given by Eq. (6). Consequently, the difference between the energy per particle of the modulated and the homogeneous phase (c 1 = 0) is The boundary between each type of modulated phase with the homogeneous phase defines an energy-crossing line through the condition where c * 1 represent the value of c 1 minimizing the energy of the modulated state n (c 1 ) and consequently ∆ n (c 1 ), once we take into account that the energy of the homogeneous phase is independent of the value of c 1 . Our definition for c * 1 implies then d∆ n dc 1 (c * 1 ) = 0, Solving Eq. (12) for c * 1 we concluded that: whereŨ = U n|v(k 0 (n, ρ))| and U and ρ represent. This result can now be substituted in Eq. (13) concluding that, at the transition,Ũ Eq. (15) implies that, in the single mode approximation, the energy crossing line between the solid and the superfluid phase is given by the conditionŨ = 4 √ 3π 2 /7. This is an universal condition, independent of the occupancy number of the solid and valid for any cluster forming potentialv(k). At the energy crossing line, the amplitude c 1 takes the universal value c 1 = 1/3, which can be interpret as a quantum analog of the Lindemann criterion.
At this point it is possible to calculate the value of the interaction strength U at the transition point from the n-particle cluster solid to the superfluid phase by using the first condition in Eq. (15). This leads to where k 0 (n, ρ) is given by Eq. (6).
For different values of n, Eq. (16) represents the melting curves of the different solid phases, and naturally they are valid in the region of densities in whichv(k 0 (n, ρ)) is negative and contains the absolute minimum ofv(k). The union of the areas delimited by the curves U n (ρ) represents the whole solid region. Interestingly, the boundary of the solid region is not a simple monotonous curve as obtained in mean-field calculations. This is in agreement with reported results of ground-state numerical simulations. [2] In order to test the single mode approximation we additionally calculated the solid-to-fluid phase boundary resulting from the minimization of the full energy function Eq. (8) for the potential v(r) = (1 + r 6 ) −1 . This is a repulsive potential expressed through the dimensionless quantities r/r c → r, v(r)/U 0 → v(r), where r c and U 0 are the units of length and energy respectively. The potential approaches a constant value for small interparticle distances r < 1, and drops to zero for r → ∞ with a repulsive van der Waals tail [26].
In Fig. 1 the result from this numerical minimization (blue dots) is plotted together with the curves U n (ρ) for occupancy numbers 1 ≤ n ≤ 5. As can be observed, the numerical results and the single mode melting curves compare very well, which indicates that in the melting region the single mode approximation is essentially correct. This is a quite useful validation because it largely simplifies all analytical calculations in the melting region. The envelope of the family of curves U n (ρ) can be also calculated and yields This envelope correspond to the result that would be obtained by considering n as a real variational parameter, which is the assumption in previous mean-field descriptions. Not to take into account the variation of the stability of the different cluster-crystal phases with fixed occupancy number can explain why many theoretical calculations find a monotonous behavior of the solid-to-superfluid phase boundary. The obtained scaling U (ρ) ∝ ρ −1 of the envelope curve is also a result that have been observed with numerical simulations [2].

B. Extent of the SSC region
For a given occupancy n, with varying density each cluster phase is limited by the inverted dome-shaped region of the curve U n (ρ). This transition from the solid to the homogeneous or superfluid state is expected to be of first order type and consequently occurs through a crossover along the coexistence region [3]. In the coexistence regions the pressure P (ρ) and the chemical potentials µ(ρ) of each phase are equal and remain the same when density is increased from the beginning to the end of the coexistence regions. The mathematical condition determining the densities ρ 1n and ρ 2n at the beginning and end of the coexistence regions is given by P n (ρ 1n ) = P n+1 (ρ 2n ) µ n (ρ 1n ) = µ n+1 (ρ 2n ).
Within the canonical ensemble formalism, pressure and chemical potential can be calculated using the relations P n (ρ) = ρ 2 ∂ n(ρ) ∂ρ and µ n (ρ) = n (ρ) + ρ ∂ n (ρ) ∂ρ [38]. To pursue the calculation of the coexistence region analytically we need to calculate the energy per particle of the cluster state as a function of the density. Even within the single-mode approximation, no simple analytical closed expression exist for the ground-state energy of the system. In order to obtain such analytical expression from Eq. (10) further approximations have to be done. First, the general relation betweenŨ and the value of c 1m minimizing de energy function Eq. (10) is found from the condition d∆ dc1 (c 1 ) = 0. This leads to the relatioñ Then the value ofŨ is used to write a closed expression for the minimum energy per particle in terms of c 1m , eliminating the interaction strength parameter U . This procedure yields Now a series expansion is made forŨ around the value of c 1m at which the energies of the cluster and the homogeneous configuration equal, fulfilling Eq. (15). Inverting this expansion we can rewrite c 1m as a power series of (Ũ − 4 7 √ 3π 2 ). This expansion can be then substituted in Eq. (20) to finally obtain a series expansion of the energy ∆ n in powers of (Ũ − 4 7 √ 3π 2 ). This procedure leads to the expression which is valid forŨ ≥ 4 7 √ 3π 2 . ForŨ < 4 7 √ 3π 2 the minimum of ∆ n is reached at c 1 = 0 and consequently ∆ n = 0.
Examining the system of equations in Eq. (18) for the densities at the beginning and the end of the coexistence region, we realize it needs as an input the quantities n (ρ) and its derivatives with respect to ρ. A general analytical solution of this equation system for the n-cluster SSC phase, for an arbitrary potentialv(k), is impossible even in the single mode approximation. We estimate it using an expansion of n (ρ) up to first order around the density at which the energy of the cluster and fluid phases equals ρ 0 (n, U ). We proceed considering where h (ρ 0 ) represents the energy of the homogeneous state. Additionally let us note that n (ρ 0 ) can be calculated using the obtained result for ∆ n (Ũ ) in Eq. (21). A straightforward calculation allows us to conclude that Taking into account the above simplifications, we directly obtain the densities at the beginning and the end of the coexistence regions To asses the accuracy of the analytical estimations on the coexistence region, in Fig. 2 a comparison is shown between the approximate solution of Eq. (24) and the exact values of the density at the beginning and end of the SSC. The exact values are determined from the exact mean-field value of n (ρ) following the procedure described in this section. In Fig. 2 the exact coexistence regions corresponding to different occupancy numbers are represented by the blue dots while the analytical estimation of Eq. (24) is represented by the continuous red curves.
As can be seen, the analytical approximation underestimates the extent of the SSC. It is possible to show that this is a general feature resulting from neglecting higher order corrections in ρ − ρ 0 (n, U ) for the expansion of n (ρ). Nevertheless, the result serves as a simple analytical lower bound for the extent of the SSC, which only exists in a narrow region.

C. Numerical phase diagram
The numerical minimization of the energy function allows us to construct the full phase diagram for the potential under consideration v(r) = (1 + r 6 ) −1 . This groundstate phase diagram is shown in Fig. 3, including not only the pure cluster phases (nCC), but also the coexistence regions between cluster solids with different occupation (SC) and the SSC. The method followed to determine the coexistence regions in this case is the same sketched in the previous section. As expected, we observe that increasing density the pure cluster phase with occupancy number n is limited from below by the coexistence region of the cluster transition from n − 1 to n and from above by the coexistence region of the cluster transition from n to n + 1.
As can be seen in Fig. 3, the SSC only exist in a narrow region between the solid and the superfluid phases. Thus, the SSC region follows the non-monotonic behavior of the energy crossing line of the solid-to-fluid transition (see Fig. 1). This behavior is well documented with Monte Carlo simulations for the supersolid phases of this potential [2]. Moreover, it was also reported that the density intervals corresponding to the coexistence of two solid phases are more likely to develop supersolidity, which is consistent with the phenomenology observed for the SSC in Fig. 3.
For comparison, we also show in Fig. 3 the superfluidsupersolid boundary as numerically obtained in Ref. [2] with dark blue dots. The correspondence between the exact localization of these points and our SSC boundary without any extra fit is remarkable, although for small density this correspondence is lost. The superfluid phase boundary is thus quantitatively and qualitatively better described by looking at the coexistence regions, which are a direct consequence from considering that clusters phases have a fixed integer number of particles per cluster.
An interesting feature regarding the phase diagram is related to the boundary of the coexistence regions. From the thermodynamical point of view a system with only one component can not present an extended region in the phase diagram in which three phases coexist. This constraint is a result of the Gibbs phases rule. At the same time, when the three coexistence regions solid(n)solid(n+1), solid(n)-superfluid and solid(n+1)-superfluid meet there are a number of thermodynamical constraints that need to be fulfilled. This makes that the topology of the phase diagram, in the region where coexistence between different phases coincides, needs to reproduce the kind of behavior shown in Fig. 4, independently of the level of approximation of the calculation.
To understand why this is a a general result we have to realize that, within each coexistence region, the chemical potential and the pressure have constant values along the density axis. At the same time, the pressure and the chemical potential need to be continuous functions of the density ρ and of the strength of the interaction U . Consider now that at some value U 0 the coexistence region of the transition from a phase I to a phase II, increasing density, meets the coexistence region of the transition from the phase II to a phase III. This means that at U 0 the pressure and the chemical potentials of the pure I and III phases, at the boundary of the whole I-II-III coexistence region, are equals. And this is exactly the conditions for the establishment of the I-III coexistence region. The above arguments imply that at U 0 the energy curves resulting from the condition ρ 2 ∂ (ρ) ∂ρ = constant, associated with the coexistence I-II, II-III and I-III, coincide. Consequently the overlap of the I-II and II-III coexistence regions with the I-III coexistence region have to be vertical, as shown in panels A and B of Fig. 4.
It is worth noting, from Eq. 16 and Eq. 24, that for a given cluster state the extent in densities of the SSC region decreases as we approach to the bottom of the dome shape curve separating the solid from the superfluid phase (Fig. 2). This mathematical result is a consequence of the fact that, at a given U , the angle between the curves n (ρ) and h (ρ) in its crossing point decreases to zero as we approach the limit of stability of the solid region. Lowering this angle implies that the curves are progressively more parallel and, consequently, the region in which we can further minimize the energy by creating a coexistence state, i.e. following an energy-density relation of the form c (ρ) = (ρ 1 ) − P c (ρ −1 − ρ −1 1 ), is consistently diminished. The functionality c (ρ) previously stated, is just a consequence of demanding that the pressure of the system within the coexistence region remains constant as the density is varied. This explanation clarifies why the extent of the SSC regions increases as we move away from the bottom of each solid dome.
The above discussion clarifies the connection between the solid coexistence states (SC) with the SSC. In short, the coexistence regions between different pure cluster phases always evolve into a SSC phase as the interaction strength is decreased. It can be ensured that this transition occurs at a fixed value of U in the whole density range of coexistence between different cluster phases and, at the moment of such transition, the extent of the SSC is the biggest possible. In this sense it can be understood that clusters coexistence favors the development of the SSC, as previous works have pointed out [2].

IV. CONCLUDING REMARKS
We have shown that the analytical understanding of the ground-state phase diagrams of ultrasoft bosons can be achieved by a variational approach when considering an integer number of particles per cluster. With this assumption, the detailed shape of the phase diagrams can be improved with respect to mean-field calculations where a real number is considered for the cluster occupation. The analytical approach introduced here is generally applicable to any ultrasoft potential.
Among the determination of the regions of coexistence of solid phases, the solid-superfluid coexistence region (SSC) is particularly relevant. The comparison with previous Monte Carlo outcomes of bosonic systems revealed that the SSC region is visually coincident with the supersolid phase. The method can be used to find nontrivial information like detailed phase diagrams and possible markers for the emergence of supersolid regimes.