The inverse-square interaction phase diagram: unitarity in the bosonic ground state

Ground-state properties of bosons interacting via inverse square potential (three dimensional Calogero-Sutherland model) are analyzed. A number of quantities scale with the density and can be naturally expressed in units of the Fermi energy and Fermi momentum multiplied by a dimensionless constant (Bertsch parameter). Two analytical approaches are developed: the Bogoliubov theory for weak and the harmonic approximation (HA) for strong interactions. Diffusion Monte Carlo method is used to obtain the ground-state properties in a non-perturbative manner. We report the dependence of the Bertsch parameter on the interaction strength and construct a Pad\'e approximant which fits the numerical data and reproduces correctly the asymptotic limits of weak and strong interactions. We find good agreement with beyond-mean field theory for the energy and the condensate fraction. The pair distribution function and the static structure factor are reported for a number of characteristic interactions. We demonstrate that the system experiences a gas-solid phase transition as a function of the dimensionless interaction strength. A peculiarity of the system is that by changing the density it is not possible to induce the phase transition. We show that the low-lying excitation spectrum contains plasmons in both phases, in agreement with the Bogoliubov and HA theories. Finally, we argue that this model can be interpreted as a realization of the unitary limit of a Bose system with the advantage that the system stays in the genuine ground state contrarily to the metastable state realized in experiments with short-range Bose gases.


Introduction
The amazing progress in the field of ultracold Bose and Fermi atoms has provided a versatile tool for a highly controlled investigation of the properties of quantum systems. The enormous advantage [1] of ultradilute quantum gases as compared to solid state systems and quantum liquids is the system high purity and controllability of the parameters, permitting one to tune the interaction strength by using a Feshbach resonance. In particular, this permits one to revisit a number of old problems in new setups. The Bardeen, Cooper, and Schrieffer theory [2][3][4] (Nobel prize in Physics in 1972), developed in the 60's to describe superconductivity in metals was more recently used to describe the BCS -BEC crossover [5][6][7][8][9][10][11][12][13]. In the original model the Cooper pairs are formed due to attraction between electrons and there is limited control over the interactions and the purity of the system. This is not the case in ultracold atoms where the Cooper pairs are formed between two spin components of alkali atoms. Another characteristic example is that of a polaron [14], a quasiparticle introduced by Lev Landau to describe properties of an electron moving in a crystal, what was recently experimentally and theoretically studied in Bose [15][16][17][18][19][20][21][22][23][24] and Fermi [25][26][27][28] gases. On the other hand, ultracold atoms also provided a way to experimentally address some problems which were otherwise unaccessible. In particular, the unitary gas has been realized in a number of experiments [29][30][31] and some of its properties are like those found in neutron stars [32].
The main idea behind universality is that two-body interactions in a dilute Fermi or Bose gas can be described by a single parameter [33], namely the s-wave scattering length a, as the mean interparticle distance is large compared to the range of the interaction potential and its details are not important. Under such conditions, the properties of the gas are governed by a single dimensionless quantity, the gas parameter na 3 . When this parameter is small, na 3 1, usual perturbative methods can be applied [1,34]. In the opposite non-perturbative regime, na 3 1, the s-wave scattering is large compared to the mean interparticle distance and it drops out of the considerations. As a result, the only physically relevant scale which is left at zero temperature is the density. Such conditions, known as a unitary regime, lead to a number of peculiar properties. For example, the energy per particle can be written as E/N = ξE F where E F is the Fermi energy, fully defined by its density, and ξ is a dimensionless constant, often referred to as a Bertsch parameter. The formal similarity between Bose and Fermi systems with density n is that the only physical length scale is defined by the interparticle distance ∝ n −1/3 , the momentum scale by the Fermi momentum k F ∝ n 1/3 and the energy unit is the Fermi energy E F =h 2 k 2 F /(2m) ∝h 2 n 2/3 /m. As a result the scaling of the properties of a Bose system with density is similar to the scaling of an ideal Fermi gas. Such a system will have thermodynamics similar to that of an ideal Fermi gas, even if there are strong interactions in the system and it lacks a perturbative expansion parameter. Thus, knowledge of a single parameter ξ provides a rather exhaustive description of the system, while its value cannot be found perturbatively. The use of Quantum Monte Carlo methods has turned out to be very fruitful for the calculation of ξ in two-component Fermi gases as such methods provide an ab initio non-perturbative approach (see Ref. [35] and references within for a historical overview).
The requirement for the unitary regime is that R l a, where R denotes the range of the potential, l ∝ n −1/3 is the mean interparticle distance and a the s-wave scattering length. For purely repulsive interactions the s-wave scattering length is smaller than the range of the potential for any finite repulsion, a < R, and a = R for the infinitely high interaction potential (hard spheres of diameter R). It means that the condition R a cannot be satisfied for a repulsive interaction. Instead, for an attractive interaction, the s-wave scattering can be large compared to the range of the potential. For a zero-range potential, the energy of a two-body bound state is related to the s-wave scattering length as E = −h 2 /(ma 2 ) where m is mass of a single atom. In the limit of a vanishing bound state, E → 0, the s-wave scattering length diverges. This means that the unitary regime essentially corresponds to the resonant scattering and it is experimentally created by tuning the magnetic field to the Feshbach resonance position. This method works well for fermions where the attraction between different components is compensated by the quantum pressure originating from the Pauli exclusion principle. The situation is drastically different for bosons where attraction between atoms leads to a many-body collapse (a bright soliton) which means that a homogeneous gas state is not thermodynamically stable. Although, recently, a unitary Bose gas was observed in a number of experiments [36][37][38], it is not the true ground state but can probably be interpreted as a metastable state with a finite lifetime. Moreover, the simple description in terms of a single parameter, the two-body s-wave scattering length, is not sufficient as Efimov three-body physics comes into play [39] and it provides an additional length scale. As a result, it has not been possible to create a ground-state unitary Bose gas. From this perspective, it is interesting to investigate the possibility of using the inverse square pair interaction in order to realize the unitary regime in the ground state of a Bose system. The key feature of the 1/r 2 interaction is that it is scale free in a system composed of quantum particles. The s-wave scattering length calculated for this interaction diverges, similar to what happens in the unitary regime for short-ranged gases. As a result, there is no any other length scale except for the mean interparticle distance. Consequently, the energy per particle can be expressed as E/N = ξE F where the Bertsch parameter ξ does not depend on the density, leading to the same thermodynamics as in an ideal Fermi gas with the Fermi energy E F .
The effects of strong correlations eventually leading to the quantum crystallization typically can be understood from analysing the value of the dimensionless quantum parameter q, defined as the ratio of the characteristic interparticle interaction energy V to the characteristic quantum kinetic energȳ h 2 /ml 2 , where l is the average distance between particles. An important feature of the inverse-square interaction Qh 2 /mr 2 is that the dimensionless quantum parameter q ∝ Q is independent of the interparticle distance or of the density n. We note that an analogous situation holds for electrons in graphene and in other Dirac materials. In these systems, due to the linear dependence of the kinetic energy from the momentum operator, the characteristic quantum kinetic energy ∝hV F /r (where V F is the concentration-independent Fermi velocity, V F ≈ c/300 with c the speed of light). Therefore, for graphene and other Dirac materials, the characteristic quantum kinetic energy scales linearly with the average distance between electrons, similar to the linear scaling of the electron-electron Coulomb interaction V(r) = e 2 /r. As a result, the governing quantum parameter is q = e 2 /(hV F ) ≈ 2 for graphene. This leads, in particular, to the impossibility of Wigner crystallization for graphene (in the absence of strong magnetic fields).
There is a crucial difference between classical and quantum mechanical systems interacting with an attractive inverse-square potential, known as a quantum anomaly [40,41]. While there are no qualitative changes in a classical system as Q is changed, in a quantum system there exists a critical value of Q = −1/4. For a weaker attraction, Q > −1/4, two interacting particles have a finite energy.
Instead, for more attractive interactions, Q ≤ −1/4, the energy becomes infinite and the ground state ceases to exist, a phenomenon known as "the fall of a particle to the center" [33]. Suppression of the quantum collapse, i.e. of the non existence of the ground state has been studied for a central 1/r 2 interaction potential in a number of systems by using mean-field methods [42][43][44] and a quantum Monte Carlo technique [45]. Very recently the quantum anomaly has been observed in two-dimensional Fermi gases [46,47].
The inverse-square potential is a famous example of exactly-solvable many body quantum systems in one dimension (see the book [48] and references within). Both the wave function and energy can be explicitly written as a function of the interaction parameter both in homogeneous and trapped geometries [49][50][51][52]. It was observed that a similar construction of the ground state wave function in two dimensions leads to equivalent statistical averages as for a system of classical charges, while the quantum Hamiltonian in addition possesses three-body interaction terms [53][54][55][56].
In the following we study the zero-temperature phase diagram of particles interacting via inverse square interactions and interpret the system properties in relation to the Bose system at unitarity.

Homogeneous Calogero-Sutherland model
We consider the following model Hamiltonian describing a system of N particles of mass m interacting via the repulsive inverse square interaction potential In numerical simulations we impose periodic boundary conditions in a box of linear size L so that the thermodynamic limit is obtained by increasing N at a fixed density n = N/L 3 . The dimensionless parameter Q defines the strength of the interaction potential. A peculiarity of the Calogero-Sutherland interaction potential is that its dependence on the interaction distance, r −2 , scales similarly to the kinetic energy term. As a result, there is no other length scale than the density. This means that by changing the density it is not possible to change the phase of the system which instead is controlled by the dimensionless parameter Q. For small values of Q it is possible to develop a perturbative theory of a weakly interacting Bose gas while in the opposite regime of a classical system with large Q we use the harmonic crystal theory.

Bogoliubov theory
The scattering problem for the potential Q = −1/8 + ν 2 /2 results in a constant phase shift δ = −(π/2)(ν + 1/2), independent of the incident momentum k [40]. In this sense, while the scattering phase is well defined, it is not possible to introduce the s-wave scattering length a, as the standard definition result in a divergent value, a = − lim k→0 tan(δ(k))/k → ∞. The independence of the phase shift δ of momentum k is another manifestation of the impossibility to introduce a unit of scale from the two-body problem. In turn, this means that the standard perturbative theory developed for Bose gases is not applicable to the inverse square interaction and appropriate Bogoliubov theory should be developed.
In the limit of Q → 0 the gas becomes weakly interacting. We develop Bogoliubov theory in this regime, based on the assumption that the condensate is macroscopically occupied. In the first quantization the Hamiltonian (1) is written aŝ In order to calculate the Fourier transform of the interaction potential we first remove the short-range divergence of the potential by screening where the screening length a defines the value of the potential at zero according to V scr (r = 0) = Qh 2 /(ma 2 ). The Fourier transform of (3) can be now evaluated and is related to the screened-Coulomb (Yukawa) function. In a homogeneous system the field operators can be expanded in the basis of plane waves, whereâ k is an operator which annihilates a particle in the single-particle state with momentum k. In a finite-sized box, only momenta satisfying the periodic boundary conditions are allowed. Substitution of (5) into (2) gives the expression of the Hamiltonian in the second-quantization form Within the Bogoliubov theory, it is assumed that as the k = 0 state is macroscopically occupied, the corresponding operator can be treated as a number,â 0 = √ N 0 and it can be used to construct a perturbative series, Here and in the following the primed summations stand for summations over all k except k = 0. In addition to the textbook expression [1,34] derived for a short-range potential, here the Fourier transform of the interaction potential at finite momentum, V k , appears. Terms of a similar nature have appeared in Bogoliubov theory for dipolar interactions [57] which are also long-ranged. A more important difference is that we do not have renormalization terms, which include correction beyond the lowest-order Born approximation to the coupling constant in terms of the s-wave scattering length to remove a divergence in the energy. The interaction potential in Eq. (1) would have a divergent s-wave scattering length and, indeed, is not a short-range potential. No second-order correction appears in our case and we immediately obtain convergent results for the Lee-Huang-Yang correction, as will be shown below.
Conceptually this is important, as there is a similarity between the 1/r 2 interaction potential and the δ-interacting potential in two dimensions as both scale as one over the distance squared. The difference is that the theory with δ-interaction potential suffers from the infrared divergency which is cured by the renormalization procedure which breaks the symmetry between classical and quantum systems. This anomalous symmetry breaking known as the conformal anomaly was predicted for two-dimensional gases in Ref. [58,59] and was very recently experimentally observed in two experiments Ref. [46,47]. Thus, it is important to note that the 1/r 2 interaction potential on the contrary does not have the conformal anomaly.
Quadratic inâ k form (7) can be diagonalized using the Bogoliubov transformation, By setting the coefficients in the off-diagonal terms to zero one obtains the Bogoliubov amplitudes where n 0 = N 0 /V is the density of the condensate. The Bogoliubov excitation spectrum is given by Bogoliubov transformation (9) is used to diagonalize Hamiltonian (6) which leads to the diagonal form,Ĥ where the ground-state energy has two contributions, the first coming from the mean-field interactions and the beyond-mean-field contribution stemming from quantum fluctuations (analogous to the Lee-Huang-Yang corrections for short-range interactions) For the inverse-square potential with no screening (a = 0), the Fourier transform V k for k = 0 diverges in the thermodynamic limit, N → ∞, resulting in an infinite energy per particle. At the same time, it is expected that the physical correlation functions must be well defined. The physically correct gapless spectrum is ensured by the following choice of the chemical potential which is consistent with the mean-field energy (15). The resulting Bogoliubov spectrum, features a square-root dependence on the dimensionless momentumk = k/n 1/3 for low-energy excitations. Such excitations can be interpreted as plasmons similar to the Coulomb case and appear due to the long-range nature of the interactions. This situation should be contrasted with the linear phononic behavior typical for short-range potentials. Quantum depletion of the condensate can be obtained by summation over the momentum distribution of particles excited out of the condensate. It follows from Eq. (9) that the zero-temperature momentum distribution is proportional to the Bogoliubov amplitude (11) as n k = â † kâ k = v 2 k . In the thermodynamic limit the sum can be approximated by an integral over the volume v, and the condensate fraction is suppressed linearly by the interaction strength Q in the weakly interacting limit. The quantum fluctuation correction (16) can be evaluated in the thermodynamic limit and is equal to where C BMF = 3 5 π 5/6 Γ( 1 3 )Γ( 7 6 ) = 3.871. Remarkably, the n 2/3 density dependence of the bosonic energy (20) is similar to that of an ideal Fermi gas.

Quantum Monte Carlo approach
Diffusion Monte Carlo method is used to study numerically the ground state properties. The statistical noise in the method can be greatly reduced by using a physically sound guiding wave function. Thus, before its construction we analyze in detail the long-and short-range behavior of the pair correlations in the system.

Long-wavelength part of wave function
Here we derive the long-wavelength limit of the many-body wave function following the recipe proposed by Reatto and Chester in Ref. [60] based on hydrodynamic grounds. Within the harmonic approximation, the long-range part of the many-body wave function can be written as where m k = m/(Nk 2 ) is the effective mass, ω k is the long-wavelength plasmonic excitation spectrum and ρ k = N ∑ j=1 e ik·r j is the Fourier transform of the density operator. It means that the long-wavelength asymptote has a Bijl-Jastrow form, where the correlation function is given by The low-momentum excitation spectrum for the screened potential (3) is obtained from the Bogoliubov excitation spectrumhω k = nV kh 2 k 2 /m with V k given by Eq. (4), resulting in the thermodynamic limit result In the limit of the unscreened potential, a → 0, Eq. (24) takes the simple form Interestingly, the asymptotic decay does not containh, as it does for the usual phononic systems, where χ(r) = (mc)/(π 2h nr 2 ).

Short-range part of wave function
When two particles approach each other, r → 0, the divergence in the interaction potential implies that the wave function cannot diverge faster than 1/ √ r in order to have finite potential energy. For corresponding distances, contribution from other particles can be neglected resulting effectively in a two-body problem. The Schrödinger equation for two particles for the relative motion is where r denotes the distance between the particles, m r = m/2 is the reduced mass and we are searching for a spherical solution. For convenience we have expressed the interaction parameter as Q = λ(λ + 1). Equation (26) can be explicitly solved for k = 0 (zero-energy scattering problem) with the solution written as a linear combination of r λ and r −1−λ terms. The latter term leads to a diverging wave function and infinite potential energy, and is to be discarded. Thus, the scattering at zero momentum results in with C some normalization constant. Consequently, the Bijl-Jastrow terms should behave like (27) for short distances. In a many-body problem, the reduced mass m b used in Eq. (26) to generate the variational wave function can be treated as a variational parameter [61] which parameterizes different variational states. It was shown in Ref. [61] that by doing so one obtains a better accuracy for the description of liquid helium which is a strongly correlated system. In our case the correlations are not so strong and we find that keeping the two-body value, m b = m/2, is sufficient for performing calculations.
The scattering solution for zero momentum, k = 0, is given by Eq. (27). Taking into account the first correction due to finite value of k > 0 we get the subleading correction, The described procedure of searching for the solution in terms of a sum at short distances by cancelling the leading powers in the potential and kinetic energies is equivalent to the cusp condition commonly used in strongly diverging potentials (for example, Lennard-Jones one). The cusp condition (28) can be recovered in a formal way from the exact solution of the Schrödinger equation (26) at a finite momentum, in terms of the spherical Bessel function of the first kind J n (x). Taylor expansion of Eq. (29) gives which, indeed, has the same functional form as the expansion (28) obtained from the cusp condition.

Guiding wave function
In order to improve the convergence we use an importance sampling technique and we chose the guiding wave function in the Nosanow-Jastrow pair form, where r c i denotes the radius vector of the each of the crystal lattice sites and α is a variational parameter governing particle localization near the crystal lattice site. When the localization is absent, α = 0, the wave function is translationally invariant and it describes the gas phase. Instead, for a finite localization strength, the wave function has broken translational symmetry and it is used to obtain properties in the solid phase.
We construct the Jastrow pair function f (r) in such a way that it combines the physics of the two-body scattering at short distances and approaches Eq. (30) as r → 0, while at large distances it approaches the hydrodynamic asymptote given by Eq. (25). We find it convenient to use the following choice of f (r), where coefficients c 1 , · · · , c 5 are defined by the conditions of the continuity of f (r) and its first derivative f (r) at the matching point R match and at L/2. The latter condition ensures that the guiding wave function (31) satisfies the periodic boundary conditions imposed on the box of size L × L × L. The matching position R match is a free parameter and its value is optimized by minimizing the variational energy.

Mean-field contribution
While the beyond-mean-field contribution (16) to the energy is finite, the leading contribution (15) is divergent. The reason is that the BMF term comes from well-behaved quantum fluctuations while the diverging term comes from the long-range contribution to the potential energy. This divergence is similar to that of a Coulomb gas, in which the divergence coming from n 0 V 0 terms is removed by imposing the condition of the charge neutrality. A reasoning similar to that of a jellium model can be applied to our case. The energy of a uniform charge will be subtracted from the total and potential energies.
It is important to note that the use of the jellium model does not change the physical properties of the original Hamiltonian (1). All correlations in the system remain exactly the same. At the same time it becomes possible to quantify the potential energy which otherwise diverges in the thermodynamic limit.
The standard Bogoliubov theory does not distinguish the total n and the condensate n 0 densities. That is, the subtraction of n 0 V 0 from the chemical potential, appearing in the Bogoliubov theory (17), is actually equivalent to the subtraction of nV 0 in the jellium model.

Direct summation
In the jellium model for Coulomb charges, energy of the opposite uniform charge is subtracted from the energy of a finite system, Eq. (17). The simplest way to implement this condition is to truncate the interaction potential spherically to V(r) = 0 for r > L/2. The finite size effects on the energy are significantly reduced by adding the missing tail energy By assuming that the pair correlation function has reached its asymptotic value, g 2 (r) = n 2 , at half the size of the box The implied condition of constant g 2 (r) is indeed satisfied in the gas phase as verified a posteriori in Fig. 4. Instead, in the crystal phase the self averaging over the peaks in g 2 (r) makes Eq. (34) applicable for calculation of the mean energy correction. The jellium model assumes a uniform charge distribution, so that its contribution to the energy exactly coincides with the tail energy for distances r > L/2. Instead, the contribution to the jellium energy coming from smaller distances definitely differs from the physical contribution at the same distances as it ignores correlations between particles. It can be explicitly evaluated resulting in Jellium contributions (34-35) will be subtracted from the total and potential energy eliminating the divergence in the thermodynamic limit. The jellium term (35) overestimates the potential energy which as it will be shown later is the dominant contribution to the total energy. As a result, the energy after the subtraction will be effectively negative while physically the total energy diverges to plus infinity in the thermodynamic limit.

Smooth version of the long-range potential
It is known that numerical simulations of systems with long-range interactions suffer from a number of technical issues. In particular, the following two challenges arise.
The first challenge to deal with is that the total energy per particle diverges in the thermodynamic limit. For example, the same problem arises in the calculation of the energy of an electron gas. Physically, the Coulomb energy of a same-charge system is diverging as the total charge becomes infinitely large in the thermodynamic limit. A standard procedure in this case consists in using the jellium model which adds the opposite charge uniformly distributed in the volume, so that the total charge remains equal to zero in the thermodynamic limit. In Monte Carlo simulations we subtract the jellium contribution (35) which permits us to calculate the energy in the thermodynamic limit. Within the Bogoliubov theory a gapless spectrum is recovered when a related term (17) is subtracted.
The second challenge is that once the diverging contribution is subtracted, the remaining integral converges slowly to the thermodynamic limit. A possible way out consists in using the Ewald summation technique [62,63] which, essentially, converts a summation over images (slow power-law convergence) to a sum over coordinate and momentum space (Gaussian convergence). Instead, in the present work we use a method of "smooth cutoff" 1 . To do so, we substitute the original interaction potential V(r) in Hamiltonian (1) by its "smooth" version V sm (r) = V(r)θ sm (r/R c ), (36) where the smoothing function is taken as The "tail" energy corresponding to a uniform average for distances r > R c is calculated analytically resulting in the following contribution which is added to the energy per particle. There are three parameters which need to be specified. Two of them correspond to the shape of the smoothing function (k 1 and k 2 ) and the third one to the cutoff distance R c . Parameters k 1 and k 2 should be large enough so that the cut-off method works effectively. In practice we find that the choices k 1 = 6 and k 2 = 2 provides good results. The cut-off parameter is taken as R c = 0.73336L where L is the linear size of the box side. The exact value of this parameter is calculated from the considering that the tail energy obtained with the smooth cutoff is equal to that of a sharp cutoff at R c = L/2, that is, lim A→∞ ( A 0 (Q/r 2 )(1 − θ sm (r/0.73336L))4πr 2 dr − A L/2 (Q/r 2 )4πr 2 dr) = 0. The obtained number (0.73336) depends on the particular choice of the other parameters (k 1 = 6 and k 2 = 2) and on the degree ν = 2 of the power-law potential V(r) = Q/r ν . We have seen that the sharp cutoff for the gas gives a fast convergence of energy to the macroscopic limit.
To summarize, we use a smooth interaction potential (37) which coincides with the original potential in the thermodynamic limit, but has a faster convergence to it. In the following we will focus on the properties of the original Hamiltonian (1) in the N → ∞ limit, obtained by the Monte Carlo method with the smooth interaction potential.

Classical limit
In the limit of large Q the potential energy dominates and the properties of the system can be analyzed by using a semiclassical approach. The system crystalizes, as happens in classical systems at zero temperature.

Equilibrium energy
The main contribution to the energy comes from the potential energy which to leading order can be estimated by considering an ideal crystal energy, where R denotes all equilibrium positions of the particles in the crystal (both within the same and other elementary cells). The energy per particle (39) calculated with the 1/r 2 interaction potential diverges leading to an infinitely large result in the thermodynamic limit. The jellium model can be conveniently used. 1 To be presented in details in a dedicated publication A possible way to do so is to limit summation in Eq. (39) to a sphere of diameter L and subtract the jellium contribution Eq. (35) from it. As a result, a finite classical energy is obtained.
The classical energy (39) is slightly different for various possible crystal packings. The optimal one corresponds to face-centered cubic (fcc) packing with the energy which is only marginally below the energy E bcc /N = −4.5369Q¯h 2 n 2 m of the body-centered cubic lattice and E hcp /N = −4.5066Q¯h 2 n 2 m of the hexagonal close-packed lattice. Such tiny differences are hard to resolve in Monte Carlo calculations. In our simulations we always assume that the crystal packing is that of a fcc lattice as it is energetically preferable in the classical limit.

Harmonic approximation
The next correction to the energy and the excitation spectrum can be calculated by using a harmonic crystal theory [64]. It is assumed that the particle positions r i are close to the corresponding lattice sites R i of the ideal crystal so that the u i = R i − r i describes small deviations from it. The potential energy can be expanded to a quadratic form in which the constant contribution is given by Eq. (39), the linear terms vanish as they should at in the minimum of the potential energy, and quadratic terms are described by the Hessian matrix Φ µν (R) = (∂ 2 U(R)/∂r µ ∂r ν ), giving, The can be conveniently introduced so that the energy can be interpreted as a quadratic form known as the harmonic approximation,

Excitation spectrum
The Newton equation of motion mü = −∂E harm /∂u can be solved according to the Bloch theorem as a wave with some dispersion law, u(r, t) = e ik·r−ω(k)t . The frequencies of the excitations can be found by diagonalizing the corresponding matrix. The correction to the leading term given by the equilibrium energy (39) is then obtained by integrating the phonon energyhω(k) over the first Brillouin zone. Here we limit ourselves to establishing the leading dependence on the interaction strength Q in the energy and on momentum k in the excitation spectrum.
Apart from some numerical coefficients, the Newton equation for the longitudinal mode can be recast in the following form 13 of 27 The low-momentum limit k → 0 corresponds to the summation over large distances. Even if the crystal has an irregular radial structure, for large distances it can be effectively neglected so that the sum can be approximated by an integral over the average density.
where R cuto f f , is some cutoff possibly needed to remove ultraviolet divergence (if present), that does not affect the low-momentum properties. For the inverse square interaction, V(r) = Qh 2 /(mr 2 ), this leads to mω 2 (k) = constQh 2 |k|/m and the low-lying excitation spectrum is A number of important properties can be noted: (i) the low-lying excitations are not linear in the momentum but rather follow a square-root dependence (ii) the strength of the excitation scales as √ Q with the interaction parameter (iii) the energy can be expressed in terms of the Fermi energy and k in Fermi momentum. Interestingly, properties (i-iii) found here in a classical system also hold true in a quantum weakly-interacting Bose gas, as given by the Bogoliubov spectrum (18). It is important to note that Planck's constant appearing ash 2 in the excitation spectrum of a classical system has no profound meaning. Indeed, according to Eq. (44) the frequency is proportional to the square root of the interaction potential which already containsh 2 . The conversion of the frequency ω to the energy E(k) requires mutliplication by anotherh. In this way the result contains the square of the Planck's constant although the underlying processes are entirely classical.
It is also interesting to check what is the dependence of the low-lying excitation spectrum (45) on momentum for other power-law potentials. It can be explicitly checked, that Coulomb 1/r charges in two dimensions also follow the law, E ∝ Q|k|. In this sense we refer to the low-lying excitations as plasmons. Instead, 1/r 3 interaction in 3D is similar to one-dimensional chain of Coulomb 1/r charges featuring a "weak" logarithmic prefactor in front of the linear dispersion term [65][66][67], E(k) ∝ | ln(k)||k|, being still a long-range potential. Instead, for the short-range potentials (1/r 4 , 1/r 5 , etc) the dispersion relation obtained with Eq. (45) is linear.
Also, we note that from Eq. (45) it is possible to conclude that the leading correction to the ground-state energy scales as √ Q and is small compared to the leading term (39) linear in Q.

Thermodynamic properties
One of the main results of the present work is the prediction for the dependence of the Bertsch parameter ξ on the interaction strength Q. In practice we extract ξ(Q) from the energy per particle obtained by diffusion Monte Carlo calculations and extrapolated to the thermodynamic limit. For large Q the Bertsch parameter grows linearly with Q and it is graphically more convenient to present ξ(Q)/Q ratio instead. This scaled quantity is reported in Fig. 1. Note that the diverging mean-field (jellium) contribution (35) is subtracted. In this way we obtain an extensive quantity which is additive with the system volume. Due to this substraction, the resulting energy is negative even if it physically corresponds to a purely repulsive system.
The gas properties in the weakly interacting regime, Q → 0, are correctly described by the Bogoliubov theory, as shown by a dashed line in Fig. 1. We find that for Q ≈ 0.1 the predictions of the Bogoliubov theory are very close to the DMC results and we conclude that the Bogoliubov theory can be safely used for obtaining the energy at smaller values of Q. The beyond mean-field correction (20) scales as ξ(Q) ∝ Q 5/3 and originates from the quantum fluctuations. It is interesting to note that the corresponding correction for a short-range Bose gas, given by the Lee-Huang-Yang term [68][69][70], is positive while we find a negative correction. The reason for this is that in the second order theory for short-range interaction one needs to take into account the renormalization of the coupling constant as opposed to the first Born approximation. Instead, the inverse-square potential corresponds to an infinite s-wave scattering length which does not need to be renormalized. The second-order perturbative theory lowers the energy as reflected by the negative sign of the correction. As Q is increases, the interactions become stronger and ξ(Q) increases in absolute terms. As Q reaches as certain value, Q * , a phase transition from gas to a solid occurs. The exact position of this transition will be commented later. For ever stronger interactions the potential energy becomes much larger than the kinetic one and in the limit Q → ∞ the energy of a classical crystal (40) is recovered.
While ξ(Q) experiences a kink at the critical value, Q = Q * , the discontinuity in the first derivative can be hardly perceived in the curve shown in Fig. 1. We find that the following Padé approximant accurately describes the thermodynamic dependence of ξ where the coefficients b 0 = 0.026; b 1 = 0.587; b 2 = 0.682; b 3 = 0.0345; a 2 = −2.557 are obtained from the fit and conditions a 1 = −C BMF b 0 , a 3 = E C b 3 are imposed to reproduce the BMF expansion (20) and the energy of a classical crystal (40) in the corresponding limits.

Quantum phase transition
The second main prediction of our work is the position of the gas-solid phase transition. In order to find its location we calculate the energy in the gas E gas and in the solid E solid phases using guiding wave functions (31) of appropriate translational symmetries. For each of the energies the thermodynamic value is obtained by using a quadratic fit in powers of 1/N for the energy with the jellium contribution subtracted. The difference between the energies in the gas and the solid phases is shown in Fig. 2. For small values of Q the gas phase is energetically preferable while for large Q the crystal phase is the ground state. The point where the difference is equal to zero corresponds to the phase transition point. We estimate its location as Q * = 94(5). It is important to note that no Maxwell double-tangent construction should be used in the present case as the phase transition is not caused by a change in the density or pressure but rather from change of an external parameter, Q. Indeed, by changing the density it not possible to provoke the phase transition and the phase always remains the same due to the scaling property. Thus, the zero-temperature phase transition is generated by changing parameter Q in Hamiltonian (2) rather than keeping the Hamiltonian and changing the volume or the number of particles. To a certain extent this is reminiscent of a situation where the phase transition is caused by a magnetic field which is an external parameter.

Coherence and structural properties.
Bose-Einstein condensation happens at low temperature in the gas phase. In order to quantify it and validate the use of the Bogoliubov theory we calculate the condensate fraction N 0 /N. We calculate the large distance asymptote of the one-body density matrix and use an extrapolation procedure from the variational and diffusion Monte Carlo estimator to minimize a possible residual dependence on the guiding wave function. The dependence of N 0 /N on the interaction parameter Q in the thermodynamic limit is reported in Fig. 3. The condensate fraction is close to unity for small values of Q showing that the condensation is complete. The departure from this value is correctly captured by the perturbative Bogoliubov theory shown with a dashed line in Fig. 3. The condensate fraction decreases monotonically as Q is increased and becomes very small close to the phase transition point. In order to study how the spatial correlations evolve with Q we measure the pair distribution function g 2 (r, r ) = n(r)n(r ) which quantifies the density-density correlations in the system. In the gas phase it is isotropic and depends on the absolute value of the distance between two points, |r − r |. A number of characteristic examples for the gas phase are reported in Fig. 4. When the distance between two particles is small, r − r → 0, the diverging interaction potential defines the wave function f (r) according to Eq. (30). For λ > 0 the wave function vanishes at the contact resulting in g 2 (r = 0) = 0. We verify that the g 2 (r) ∝ f 2 (r) ∝ |r| 2λ for small distances and arbitrary interaction strength Q. Generally, this polytropic behavior is non-analytic unless the interaction parameter corresponds to the even powers of λ. A similar situation happens in a one-dimensional Calogero-Sutherland model [71].
The pair distribution function shows a smooth monotonic behavior in the Bogoliubov limit of small Q, typical for weakly interacting Bose gases. As Q is increased, correlations becomes stronger and peaks start being formed. Close to the transition point, the height of the peak is around 1.5 in units of the average density.
The pair correlations in momentum space can be quantified by the static structure factor S(k) which is related to the Fourier transform of g 2 (r). The results for the gas phase are presented in Fig. 5. In the weakly interacting regime, Q 1, the static structure factor is a monotonous function which grows from zero at k → 0 to its asymptotic large momentum value S(k) = 1. Instead, in the regime of strong correlations, Q 1, a peak forms. The height of the peak grows as Q is increased which can be viewed as a precursor of the crystallization happening at the critical point Q * .
Importantly, the small momentum part of the static structure factor is not linear, S(k) ∝ k, as it happens in systems with a "usual" sound nor quadratic, S(k) ∝ k 2 , as in systems with a gap in the excitation spectrum [72]. Instead, the plasmonic k 3/2 dependence is found. This plasmonic behavior is observed for different values of Q. In the Bogoliubov regime it is possible to calculate the corresponding prefactor, In the solid phase the static structure factor has a number of macroscopic peaks located at the corresponding momenta of the crystal lattice. The height of the peaks increases linearly with the number of particles signaling presence of the diagonal long-range order.   (30)). Note that the point with Q = 100 corresponds to a metastable state, as the ground state in this regime is a solid.

Excitation spectrum and plasmons
A peculiarity of the model described by Hamiltonian (1) is that the excitation spectrum is gapless but the low-lying excitations are not described by linear phonons but rather by plasmons with a square root dispersion relation. The Bogoliubov excitation spectrum (12) provides an explicit expression in the limit of weak interactions.
In order to verify that the excitation spectrum follows the unusual plasmonic dispersion relation, we use the Feynman relation between the static structure factor S(k) and the excitation spectrum E(k), Estimation (48) provides a rigorous upper bound for the lower boundary of the excitation spectrum. This expression becomes exact when the excitation spectrum is exhausted by a single type of excitation as happens in the limit of plasmons for k → 0. We present the excitation spectrum as calculated from the Feynman relation in Fig. 6 for a number of characteristic values of the interaction parameter Q. The results are presented on a double logarithmic scale so that any power-law dependence appears as a straight line. We observe that for small momentum the excitation spectrum corresponds to plasmons with √ k dispersion relation. For the smallest considered interaction strength, Q = 0.1, there is an excellent agreement with the Bogoliubov excitation spectrum given by Eq. (12). For stronger interactions, the plasmonic behavior can be seen for small momenta, E(k) ∝ √ k. The coefficient of proportionality grows as Q is increased reflecting stronger interactions in agreement both with the Bogoliubov theory and prediction Eq. (45) obtained within the harmonic crystal approximation. According to the Landau criterion for superfluidity, a condensate moving with a group velocity V which is smaller than min k E(k)/k (i.e. the speed of sound c for linear spectrum with no roton) remains energetically stable. The "speed of sound" calculated as a first derivative of the dispersion relation diverges, c → ∞, which means that according to the Landau criterion the homogeneous gas state is superfluid. For large momentum the excitation spectrum given by Eq. (48) tends to the free particle spectrum, E(k) =h 2 k 2 /(2m), shown in Fig. 6 with a thin solid line. Close to the transition point, a roton minimum is formed at momenta similar to these where the Bragg crystal peak will be formed in the solid phase.

Critical parameters
It is illustrative to compare the properties of the system at the transition point to those of other quantum systems in order to check a possible universality. In particular the Lindemann parameter, the height of the peak in the static structure factor and the condensate fraction usually are of the same order as compared to other three-dimensional systems.
The Lindemann ratio, γ = r 2 /d, quantifies the strength of the mean square fluctuations of a particle close to the lattice site, r 2 , compared to the nearest neighbor separation in the lattice (d = a 0 / √ 2 for fcc lattice with elementary unit cell length a 0 ). The gas-solid phase transition is produced when the Lindemann ratio reaches a critical value which depends on the dimensionality and statistics of the system and to a lesser extent on the interaction potential.  (18); thin solid line, free particleh 2 k 2 /2m energy. Units ofh 2 n 2/3 /m are used for the energy and n 1/3 for the momentum. Table 1. Literature overview for the typical values of the Lindemann parameter γ (solid phase), height of the peak of the static structure factor S peak (k) and condensate fraction N 0 /N (gas phase) at the critical point of zero-temperature gas(liquid)-solid phase transition.

Reference
Authors The Lindemann ratio is approximately constant along the transition line and is relatively independent of the types of interaction potential and the crystal packing. For example, the phase diagram for the Yukawa potential is governed by two parameters and the explicit calculation of the transition line [85] is very close to the prediction based on a constant value of γ = 0.28(2) [76]. Thus, the value of the Lindemann ratio provides important information on the phase transition location and its value can be compared with what is observed in other systems. The Lindemann ratio is also approximately constant in classical systems at the transition, although its typical values 0.1 γ 0.15 are smaller compared to the quantum systems [86] Table 1 summarizes the value of the Lindemann parameter at the zero-temperature phase transition point in a number of different systems in three (3D) and two (2D) dimensions. It can be noted that even if the interaction potentials might be very different, including long-range ones, the actual value of the Lindemann ratio is limited to a rather narrow range, 0.23 < γ < 0.29. We find that at the transition the Lindemann parameter is equal to γ = 0.24 (1). In other words the present system falls into the same class of quantum phase transitions. Instead, classical systems have typically a smaller value of the Lindemann ratio at gas-solid transition [87].
It might be noted that also the gas phase has same parameters which are approximately constant at the critical point, which are the height of the peak in the static structure factor and the condensate fraction. We find that the condensate fraction is 0.8(2)% at the transition point and the height of the peak of the static structure factor is 1.63(5).

Universal scaling properties
A number of system properties are universal in that they can be mapped to the properties of an ideal Fermi gas and are defined by a single function, which can be formally introduced as the dependence of the Bertsch parameter ξ on the interaction parameter Q. The total energy per particle can be generally written as where the energy unit, apart from a numerical factor, coincides with the Fermi energy and we explicitly add the energy of the jellium background. The chemical potential can be calculated from the total energy (49) according to µ = ∂E/∂N resulting in The Bertsch parameter reported in Fig. 1 is obtained by subtracting the diverging jellium contribution from the energy per particle. The same is true for the chemical potential, once the jellium contribution is subtracted, the remaining term follows the ideal Fermi gas scaling. While the mean-field (jellium) contribution can be conveniently subtracted from the energy and the chemical potential as the potential energy can be calculated with a certain offset, the jellium contribution to the compressibility becomes crucial to the small-momentum properties. The pressure P = − ∂E/∂V| N is the sum of the "Fermi gas" contribution 2 3 ξ(Q)h 2 n 5/3 /m and the diverging "jellium" contribution ∝ N 1/3 Qh 2 n 5/3 /m. The latter term has a dramatic effect on the compressibility κ which at zero temperature can be calculated as κ −1 = (V/m 2 ) ∂ 2 E/∂N 2 V . Due to the jellium contribution its value vanishes in the thermodynamic limit, κ → 0. The relation between compressibility and the speed of sound, c 2 = mn/κ, results in a diverging value of the speed of sound c, seen in Fig. 6 as an infinite slope of the plasmonic excitation spectrum for small momentum. A related effect was observed in classical simulations of two-dimensional charges in a large trap [88] where long-range interactions resulted in vanishing compressibility leading to a larger concentration of charges at the border of the trap.
The potential energy E pot can be obtained by using the Hellmann -Feynman theorem by noting that Q enters in the Hamiltonian (1) only in the interaction energy, so that the potential energy per particle is obtained by differentiating the Hamiltonian with respect to Q and exchanging the order of the derivative and averaging The potential energy per particle (51) diverges in the thermodynamic limit (N → ∞ taken with N/V − const) due to the long-range nature of interactions. At the same time, the kinetic energy per particle remains finite and scales ash 2/3 n 2/3 /m, similar to the kinetic energy of an ideal Fermi gas. The kinetic energy appears due to a non-linear dependence of the Bertsch parameter on Q. In the limiting case of a classical crystal, Q → ∞, the Bertsch parameter becomes linear according to Eq. (40) so that the leading contribution to the energy of a classical crystal comes from the potential energy. The beyond-mean field energy per particle is given by Eq. (20), consequently for small Q, E kin /|E BMF pot | = 2/5.

Considerations for experimental realization
Although no explicit realization of the system under study is known to us, a number of closely related systems already exist or can be experimentally realized in the near future.
There is a close analogy between the statistical properties of zero-temperature one-dimensional quantum gases and terraces on crystal surfaces [89][90][91][92]. Typically, the crystal surface is covered by molecular layers of the same height (terraces). Different terraces are separated by steps at which the elevation is changed by the height of a single elementary cell. The border of a step seen from above draws a trajectory on a two-dimensional (x, y) surface and can be interpeted as a word line of a quantum 1D line in the (x, τ) plane where imaginary time τ plays the role of another dimension. Energetically it is not favorable to have a step of double height so that the steps (and equivalently the world lines) do not cross each other, creating an analogy with the Pauli exclusion principle. The mass for quantum particles is than mapped to the step stiffness β while thermal energy k B T replaces the Planck's constanth [91]. The elastic repulsion when meandering is modest (variation in x is small compared to the average distance between steps ) leads to energy AT 2 /(β 2 ) where A is a constant [90,92] fixed by the material. The mapping to the quantum system results in particles interacting via an inverse square potential. The dimensionless parameter of the problem,Ã = Aβ/(k B T) 2 is equivalent to Q in our terminology and it changes in a wide range, 0.1 A 100 [91] in real materials. The mapping with the Calogero-Sutherland model turns out to be useful for comparison of the energy and density-density (step-step) correlation function [89][90][91][92].
Another class of systems where the inverse square interaction potential can be experimentally realized are ions with controllable spin-spin interactions [93] relevant for the creation of gates needed in a quantum computer. It was experimentally demonstrated that power-law interactions, V(r) ∝ 1/r α , with 0 ≤ α ≤ 3 can be induced in a two-dimensional triangular crystal lattice of hundreds of particles [94], including "monopole-dipole" α = 2 case corresponding exactly to the inverse square interactions. One-dimensional ion chains with power-law interactions with 0.75 ≤ α ≤ 1.75 were experimentally realized in Ref. [95] and observed by different behavior of the light-cone as a function of α [96]. It is rather probable that the required interaction potential will be also created in three dimensional ion lattices.
The inverse-square interaction is relevant for Rydberg atoms [97] and polar molecules [98] as well as polymer physics [99]. The renormalization-group theory was used to predict properties in arbitrary number of dimensions [100][101][102]. Also the relation between the Laughlin state and the wave function of Calogero-Sutherland model was shown in Refs. [55,103,104].

Conclusions
We have studied the ground state properties of a three-dimensional quantum system with particles interacting via an inverse squared pair potential of strength Q. Its intrinsic property is that it is scale-free with the density being the only parameter providing a length scale. The system properties can be divided into two categories. The first one corresponds to mean-field quantities which are sensitive to the long-range nature of the interaction potential and which diverge in the thermodynamic limit, including total energy, potential energy, speed of sound, etc. The second category describes beyond-mean field properties which remain finite in the thermodynamic limit, including the kinetic energy, excitation spectrum, condensate fraction, Lindemann ratio, etc. The properties belonging to the second category can be naturally expressed in terms of the Fermi momentum and the Fermi energy. Diffusion Monte Carlo method is used to calculate numerically the ground-state properties for a wide range of parameters, 0.1 < Q < 1000. The guiding wave function is constructed from a power-law two body solution at short distances and a "plasmonic" long range Jastrow tail. We demonstrate that this system possesses a gas-solid phase transition. The energies of the gas and solid phases are separately calculated. We estimate the critical value of the interaction parameter as Q * = 94 (5). Notably, the density dependence of the solid phase is still that of an ideal Fermi gas.
For weak interactions, Q → 0, we develop a perturbative approach based on Bogoliubov theory, complemented with a jellium model which is used to remove the mean-field divergence of the energy in the thermodynamic limit. For strong interactions, Q → ∞, we use the harmonic crystal theory. A peculiarity emerging from the long-range nature of the interactions is that the low-lying excitations are gapless plasmons with a square root dispersion relation which is demonstrated within Bogoliubov theory and the harmonic crystal approach. According to the Landau criterion the homogeneous gas is superfluid at zero temperature.
We validate the results derived within the Bogoliubov theory for small Q and the harmonic approximation for large Q. In particular, we show that the energy correction is well reproduced by beyond-mean-field terms for small values of Q and approaches the energy of a classical crystal when Q → ∞. We verify that condensate depletion is caused by quantum fluctuations in the limit of weak interactions and becomes very small in the gas phase close to the transition point. A number of characteristic examples of the pair-distribution function g 2 (r)are reported across the gas phase. We show that for small distances it follows a non-analytic law, g 2 (r) ∝ |r| 2λ . The correlations in the system are quantified by the static structure factor. The excitation spectrum is approximated by using the Feynman relation and shows the plasmonic dispersion for small momenta. In the limit of weak interactions it coincides with the Bogoliubov excitation spectrum.
The unitary scaling with the density permits us to find intrinsic relations between different thermodynamic quantities. In particular, by using the Bertsch parameter ξ(Q) and its derivatives, it is possible to predict the energy, chemical potential and pressure. One of the main results is the prediction for the Bertsch parameter ξ(Q) as calculated from the ground state energy, E/N = ξ(Q)h 2 n 2 /(2m), for which we provide a Padé approximant. By using Helmann-Feynman theorem we also provide explicit expressions for the potential and kinetic energy. We verify such predictions by using Monte Carlo data which demonstrates the internal consistency of the obtained results.
As a consequence of the scaling properties, the dynamics is equivalent to that of an ideal Fermi gas and this model can be interpreted as a realization of a unitary regime in a Bose system. Importantly, in our case the particles stay in a genuine ground state and not in a metastable state, as instead happens in experiments with short-ranged Bose gases.