Finite Temperature Properties of a Supersolid: a RPA Approach

We study in random-phase approximation the newly discovered supersolid phase of ${}^4$He and present in detail its finite temperature properties. ${}^4$He is described within a hard-core quantum lattice gas model, with nearest and next-nearest neighbour interactions taken into account. We rigorously calculate all pair correlation functions in a cumulant decoupling scheme. Our results support the importance of the vacancies in the supersolid phase. We show that in a supersolid the net vacancy density remains constant as function of temperature, contrary to the thermal activation theory. We also analyzed in detail the thermodynamic properties of a supersolid, calculated the jump in the specific heat which compares well to the recent experiments.


Introduction
One of the biggest accomplishments of theoretical condensed matter physics is the ability to classify various phases and phase transitions by their mathematical order. These mathematical orders are usually expressed by order parameters reflecting certain limiting behaviours of two particle correlation functions. This concept makes it possible to describe the properties of physically very different systems in a common language and establish a link between them. Without this concept, the counter-intuitive idea of a supersolid, i.e. a solid that also exhibits a superflow, would not have been conceivable. In this language supersolidity, firstly proposed by Andreev and Lifshitz [1] in 1969 and picked up by Leggett and Chester in 1970 [3,2] is the classification of systems that simultaneously exhibit diagonal and off-diagonal long range order. Yet the idea of a supersolid was reluctantly received by the scientific community, because early experiments failed to detect any such effect in Helium-4. Apart from John Goodkind [4] who had previously seen suspicious signals in the ultrasound signals of solid helium, physicists were surprised when in 2004 Kim and Chan [5,6] announced the discovery of supersolid helium. Although equipped with a head start of more than 30 years the theoretical understanding of the supersolid state lags vastly behind, not least because supersolidity as we observe it today is nothing like what the pioneers in early 1970's had anticipated. It is now evident that defectons and impurities play a crucial role in the formation of the supersolid state. However, the data and the results of various experiments draw a rather complex picture of the supersolid phase. The normal solid to supersolid transition is not of the usual first or second order type transition but bears remarkable resemblance to the Kosterlitz Thouless transition well-known in two dimensional systems. Furthermore annealing experiments show that 3 He impurities play a significant role but the data remain all but conclusive, as the measured superfluid density varies by order of magnitudes among the different groups. Recently a change in the shear modulus of solid helium was found and the measured signal almost identically mimics the superfluid density measured by Kim and Chan [5,6]. Some popular theories give plausible explanations of some aspects of the matter. That, such as the vortex liquid theory suggested by Phil Anderson [9], is in good agreement with properties of the phase transition; theories based on defection networks are capable of describing the change in the shear modulus. However, we feel that a satisfying and comprehensive theory is still missing. Part of the problem stems from the complexity of the models. Many models are only at sufficient accuracy solvable with numerical Monte Carlo methods. The results, doubtlessly useful, are seldom intuitive and the lack of analytical results do not meet our notion of the understanding of a phenomenon. Other approaches on the other hand are limited in their analytical significance and rather represent a phenomenological ansatz. In this work we attempt to fill a part of this gap in the theoretical understanding of the supersolid phase. We present a theory of supersolidity in a Quantum-Lattice Gas model (QLG) beyond simple mean-field approaches. Following the approach of K.S. Liu and M.E. Fisher [12] we map the QLG model to the anisotropic Heisenberg model (aHM).
The method of Green's functions proves to be very successful in the description of ferromagnetic and antiferromagnetic phases and we use this method to investigate the supersolid phase which corresponds to a canted antiferromagnetic phase. Applying the Random-phase Approximation (RPA), we broke down using cumulant decoupling, the higher order Green's functions to obtain a closed set of solvable equations. This method gives a fully quantum mechanical and analytical solution of a supersolid phase. We will see that quantum fluctuations have a significant effect on the net vacancy density of the supersolid and we will also see that the superfluid state is unstable against zero-energy quasi-particle excitations. We also derive important thermodynamic properties of this model and derive formulas for interesting properties such as the power law exponents and the jump of the specific heat across the normal solid to supersolid line. The paper is organized as follows: In Section 2 we introduce the generic Hamiltonian of a bosonic many body system and discretize it to a quantum lattice gas model. This model is equivalent to the anisotropic Heisenberg model in an external field and we will identify the corresponding phases in Section 3. In Section 4 we re-derive the classical mean-field solution and discuss briefly their significance before we in Section 5 recapitulate the Green's functions for the anisotropic Heisenberg model in the random-phase approximation at zero temperature. The ground state of the system is obtained by solving the corresponding selfconsistency equation. In the following two sections we derive basic thermodynamic properties as well as ordinary differential equations to calculate the first and second order phase transition lines. In Section 8 we analyse the quasi-particle energy spectrum and in the last section we calculate phase diagrams for various parameter sets and analyse the properties of the supersolid phase and the corresponding transitions.

Generic Hamiltonian and Anisotropic Heisenberg Model
If we neglect possible 3 He impurities supersolid Helium-4 is a purely bosonic system whose dynamics and structure are governed by a generic bosonic Hamiltonian: Here ψ † (x) are the particle creation operators and ψ(x) destruction operators and obey the usual bosonic commutation relations. Hamiltonians such as in Eq. (1) are not, due to the vast size of the many body Hilbert-space, diagonalisable even for simple potentials V (x). The accomplishment of any successful theory is to find an approximation that sustains the crucial mechanism while reducing the mathematical complexity to a traceable level. Here we follow the method of Matsubara and Matsuda [14] who successfully introduced the Quantum Lattice Gas (QLG) model to describe superfluid Helium. We believe that the quantum lattice gas model is particularly useful for supersolids as the spatial discretization of this model serves as a natural frame for the crystal lattice of (super)-solid helium and a bipartite lattice elegantly simplifies the problem of breaking translational invariance symmetry for states that exhibit diagonal long range order. According to Matsubara and Tsuneto [13] the generic Hamiltonian Eq. (1) in the discrete lattice model reads: Here u ij are solely finite for nearest neighbor and next nearest neighbor hopping and otherwise zero. The values of u nn and u nnn are such that the kinetic energy is isotropic up to the 4th order. As atoms do not penetrate each other there can exist only one atom at a time on a lattice site and consequently a † and a are the creation and annihilation operators of a hard core boson commuting on different lattice sites, but obey the anti-commutator relations on identical sites. Equation (2) is the Hubbard model in 3 dimensions for hard core bosons. Neither fully bosonic nor fermionic, the lack of Wicks theorem inhibits the application of pertubative field theory though hard core bosonic systems are algebraic identical to spin systems. A simple transformation, where the bosonic operators are substituted by spin-1/2 operators generating the SU(2) Lie algebra, maps the present QLG model onto the anisotropic Heisenberg model: Here the correspondence between the QLG parameters and the spin coupling constants is given by:

Phases
The self-consistency equations in the random-phase approximation as we will derive in a later chapter are very lengthy and therefore it is our primary goal to keep the present model as simple as possible while still being able to describe the crucial physics. For this reason we will define the anisotropic Heisenberg model on a bipartite BCC lattice which consists of two interpenetrating SC sub-lattices as Figure 1 shows. In this lattice geometry the perfectly solid phase is composed of a fully occupied (on-site) sublattice A while sublattice B refers to the empty interstitial and is consequently vacant. As there is no spatial density A B Fig. 1. The BCC lattice consists of two interpenetrating SC sub-lattices. In the perfect solid phase one sub-lattice (i.e. sublattice A) serves as on-site centers and is occupied while sublattice B represents the empty interstitial and is vacant (For simplicity, we only present the two dimensional case). variation in the liquid phases both sublattices are equally occupied, and the mean occupation number simply corresponds to the particle density. We are aware that the above choice of SC sublattices does not reflect the true 4 He crystal structure which is hexagonal closed-packed (HCP). Nevertheless we believe that crucial physical properties such as phase transition and critical constants do not depend on the specific geometry as long as no other effects such as frustration are evoked. Table 3 charts the various magnetic phases and identifies the corresponding phases of the 4 He system. According to their spin configurations we identify the four magnetic phases to be of ferromagnetic (FE), canted ferromagnetic (CFE), canted anti-ferromagnetic (CAF) and antiferromagntic (AF) orders. The off-diagonal long range order parameter m 1 and the diagonal long range order m 2 are: which readily identify the corresponding phases of the Helium system.

Mean-Field Solution
In our previous work [10] we have already re-derived the classical mean-field solution of the anisotropic Heisenberg model at zero temperature. As this model provides an easy and intuitive access to the model we extend the approximation to finite temperatures as was done by K.S Liu and M.E. Fisher [12] and briefly discuss its solution and phase diagram. The mean-field Hamiltonian is obtained by substituting the spin-1 2 operators with their expectation values: Here J 1 = −q 1 J i∈A,j∈B , J 2 = −q 2 J i∈A,j∈A , J ⊤ 1 = −q 1 J ⊤ i∈A,j∈B and J ⊤ 2 = −q 1 J ⊤ i∈A,j∈A where q 1 = 6 and q 2 = 8 are the number of nearest and next nearest neighbors. The mean value of S y drops out as the randomly broken symmetry S x ↔ S y (ODLRO) allows for S y = 0. The standard method of deriving the corresponding selfconsistency equations is to minimize the Helmholtz's Free energy F = H − T S. The entropy S is given by the pseudo spin entropy of the system: where S A = S zA 2 + S xA 2 and S B = S zB 2 + S xB 2 . In the canted anti-ferromagnetic and the canted ferromagnetic states there are 4 self-consistency equations in the ferromagnetic and anti-ferromagnetic phases where S x = S y = 0 they are reduced by two. These equations are readily obtained by differentiating the free energy with respect to S z and S x respectively. The resulting equations can be rearranged to yield: (9) where The two equations (Eq. (9)) are dismissed in the ferromagnetic and the anti-ferromagnetic phases as the transversal fields S x A and S x B become zero. At zero temperature T = 0 the free energy and H coincide. This allows us to deduce the phases at absolute zero in the limiting cases (limits of h z ) from Equation (6). In the limit of h z → ∞ the Hamiltonian reduces to an effective one particle model: Consequently the system will be in the energetically favorable ferromagnetic phase. In the opposite limit, h z → 0 and with antiferromagnetic coupling J 1 < 0 the term of nearest neighbor interaction is the only term that significantly lowers the energy. Therefore the ground-state is the anti-ferromagnetic state. At T=0 and a suitable choice of parameters the canted ferromagnetic and the canted anti-ferromagnetic phases are realized in-between these limits as seen in Figure 2. However with increasing temperature the regions of canted ferromagnetic and canted anti-ferromagnetic phases deplete and at sufficiently high temperatures only the antiferromagnetic and ferromagnetic phases survive. As we can see from the phase diagram in Figure 2 most transitions are of second order. Only for parameter regimes where the CAF does not appear the resulting CFE-AF is of first order. The boundary lines are defined by ordinary differential equations and generally need to be calculated numerically. Nonetheless there exist analytical expressions for the locations of the transitions at absolute zero as derived by Matsuda and Tsuneto [13]. The ferromagnetic to canted ferromagnetic transition point is determined by Equations (9) if we set S zA = S zB = 1 2 and S xA = S xB = 0: Equally the canted anti-ferromagnetic to anti-ferromagnetic transition is defined by S zA = − S zB = 1 2 and S xA = S xB = 0: The canted ferromagnetic and the canted anti-ferromagnetic phases coexist where the order parameter of the DLRO ,m 1 = S zA − S zB approaches zero. We replace S zA and S zB in Equation (9) with m 1 and m 2 = S zA + S zB and retain only linear terms of m 1 . Subtracting and adding both equations respectively yields: We used that S zA 2 + S xA 2 = 1 4 at T=0. The solution of these two equations determine the corresponding transition point which is given by:

Green's Functions
Although the classical mean-field theory is quite insightful and gives an accurate description of the variously ordered phases its fails to take quantum fluctuations and quasi-particle excitations into account. Hence, in order to overcome these shortcomings and to derive a fully quantum mechanical approximation we solve the anisotropic Heisenberg model in the random-phase approximation which is based on the Green's function technique. At finite temperature the retarded and advanced Tyablikov [15,16] commutator Green's function defined in real time are: The average involves the usual quantum mechanical as well as thermal averages. The most successful technique of solving many body Green's function involves the method of equation of motion which is given by: The commutator [S x i , H] can be eliminated by using the Heisenberg equation of motion giving rise to higher, third order Green's functions on the RHS. In order to obtain a closed set of equations we apply the cumulant decoupling procedure which splits up the third order differential equation into products of single operator expectation values and two spin Green's functions. The cumulant decoupling [17] is based on the assumption that the last term of the following equality is negligible: The set of differential equations is not explicitly dependent of the temperature, i.e. temperature dependence solely comes with thermal averaging of the single operator expectation values. Therefore the solution is formally identical to the zero temperature solution and the detailed derivation of the Green's function in their full form can be found in Ref. [11]. The averages of the spin operator, appearing in the cumulant decoupling scheme, have to be determined selfconsistently. Two self-consistent equations can be derived from correlation functions corresponding to the Green's functions. The self-consistency equations of the canted ferromagnetic (superfluid) and canted antiferromagnetic (supersolid) phases are structurally different from the ferromagnetic (normal solid) and antiferromagnetic (normal fluid) phases, as the off-diagonal long range order increases the number of degrees of freedom by two and hence two additional conditions, resulting from the analytical properties of the commutator Green's functions, apply. Thus, the self-consistency equations of the canted phases, can be written as a function of the temperature, the external magnetic field and the spins in x-direction: similar the self-consistency equations for the ferromagnetic and anti-ferromagnetic phases: These equations involve a 3 dimensional integral over the first Brillouin zone. As explained in the work [10] on the zero temperature formalism this integral can be reduced to a two dimensional integral by introducing a generalized density of state (DOS) and gives the model a wider applicability.

Thermodynamic Properties
The relation between then QLG and the anisotropic Heisenberg model is such that the chemical potential µ corresponds to the external field h z , i.e. the grand canonical partition function in the QLG corresponds to the canonical partition function in the anisotropic Heisenberg model where the number of spins is fixed. Consequently, thermodynamic potentials of the two models are related as: Here Θ refers to an arbitrary thermodynamic potential.
In the same way as the ground state minimizes the internal energy U = H at absolute zero, the free energy F = U − T S is minimized at finite temperatures. We wish to stress that the internal energy in the present approximation cannot be derived accurately from the expectation value of the Hamiltonian in the following way: The cumulant decoupling, though a good approximation to the anisotropic Heisenberg model, is also an exact solution of an unknown effective Hamiltonian H ef f . Therefore thermodynamically consistent results are only obtained if the effective Hamiltonian is substituted in the equation above, i.e. U = H ef f . Here, as we do not know the explicit form of this effective model, we have to integrate the free energy from thermodynamic relations: This equation allows us to select the ground state in regions where two or more phases exist self-consistently according to the random-phase approximation. The entropy of the spin system is given by log( 1 1 + exp(−βω(k)) ) This formula is of purely combinatorial origin and reflects the fact that the hard-core boson system is equivalent to a fermionic system given by the Jordan-Wigner transformation. The ω(k) terms refer to the energies of the spin-wave excitations. In the solid phases both branches have to be taken into account. The entropy of the anisotropic Heisenberg model given for a fixed number of spins corresponds to the entropy of the QLG model at a constant volume. Therefore, in order to obtain the usual configurational entropy of the QLG, we have to divide by the number of particles per unit cell: (S conf = 2S nA+nB ). Here n A = 1/2 − S z A and n B = 1/2 − S z B are the particle occupation numbers on lattice sites A and B. As the nature of the normal solid to supersolid phase transition is not yet satisfactorily understood recent experiments [19] have focused on the behavior of the specific heat across the transition line in the hope of shedding light on the matter. The specific heat at constant temperature and constant pressure respectively are given by: Although the external magnetic field h z in the spin model is an observable the corresponding quantity in the QLG model, namely the chemical potential is not. Therefore we are also interested in attaining a formula for the pressure associated with a certain chemical potential. The relationship is most easily derived from the following Maxwell relation: where ǫ := S z A + S z B . Note that, using this equation in order to obtain the pressure at any specific temperature we need a reference point, i.e. a chemical potential where the corresponding pressure is known. This point is given by µ → ∞ which corresponds to N → 0 and, hence, P → 0. Consequently, in order to obtain the pressure for a specific chemical potential µ ′ we have to integrate over the interval [∞, µ ′ ].

Boundary Lines
The state of the system at any point in T and h z is given by the self-consistency equations and the free energy as can be derived from Eq. (24). Nevertheless the resulting computations come at high computational cost and therefore it seems most feasible to derive ODEs which determine the first and second order transition lines. First we will derive the ordinary differential equations which define the more abundant second order transitions. The normal fluid (FE) and the normal solid (CFE) phases are determined by Equations (21) and the supersolid (CAF) and superfluid (CFE) are defined by Equations (20) and condition Equation (9). Consequently on the SS-NS and SF-NF transition line, where both the normal (FE and CFE) and the super (CFE and CAF) phases coexists following equations hold: On the SS-NS boundary line the quotient S xA / S zA is not known a priori and therefore we eliminate it in the equation above yielding: We introduce a variable s which parametrizes the boundary curve. If we for instance choose ds=dT we get the following set of ordinary differential equations, defining the NF-SF and the NS-SS transition lines:  Upon crossing the SF-SS transition line, coming from the superfluid phase the set of possible solutions branches off into two phases, the supersolid and a non-physical (complex valued) superfluid phase. Therefore any matrix of ordinary differential equations will render a singularity and consequently we have to approach the transition line in a limiting process: The resulting ODE is which defines the superfluid to supersolid transition. As mentioned previously, there are certain parameter regimes where not all four possible phases are appearing and consequently a first order transition (mostly superfluid to supersolid) will occur. In the previous chapter we have seen that such a transition line is difficult to locate. However, a tricritical point frequently appears in the phase diagram and at this point the first order transition evolves into a second order transition. This tricritical point can be taken as a initial value for a differential equation defining the corresponding first order transition line. The relevant ODE may be derived from a Clausius Clapeyron type equation. On the transition line both phases have equal free energy. Hence: where S refers to the spin entropy as derived in the previous section (Eq. (25)) and ∆ refers to the difference of either entropy or spin mean-field of the superfluid and the normal solid phases. This equation together with ds = dT and the total derivative of the two self-consistency equations for the normal solid and one equation for the superfluid phase form a set of 5 ODEs determining S x in the superfluid phase and S zA and S zB in the normal solid phase along the boundary line in the T − h z plane:

Excitation spectrum
The superfluid phase features, due to spontaneously broken U(1) symmetry, the well know gapless Goldstone bosons, i.e. linear phonons. The supersolid phase additionally exhibits a second, gapped branch which is due to the break down of discrete translational symmetry. Figure 4 reveals a zero frequency mode in the supersolid phase at [100] of the first Brillouin zone and consequently the superfluid to supersolid phase transition is characterized by a collapsing roton minimum at [100]. The dispersion relation in the superfluid (CFE) phase is given by: From this equation we can see that the energy possibly goes to zero at [100] (corresponds to γ 1 (k) = −1 and γ 2 (k) = 1) when following condition is fulfilled: Hence we obtain a further condition (supplementary to Eq. (16)) for the existence of the superfluid to supersolid transition. Equation (35) allows for the existence of a second region of the reciprocal lattice space where the dispersion relation might go soft. For γ 1 (k) = 0 and γ 2 (k) = −1 which corresponds to [111] we obtain following condition: It is also interesting to study the behavior of the excitation spectrum with increasing temperature. In a conventional superfluid the long wave-length behavior is given by: Here V(0) is the interaction potential at zero momentum and m is the particles' mass. The density of the condensate n o (T ) typically decreases with increasing temperature and for that reason we expect lower energies with increasing temperature. In Figure 3 we can see that the quasiparticle energies indeed decrease with increasing temperature. Apart from the region in the vicinity of [100] the energies at higher temperatures lie significantly lower than those ones closer to absolute zero. This is important as it will contribute to the variation of thermodynamic quantities such as the entropy or the specific heat. Figure 4 depicts the variation of the excitation spectrum with increasing temperature in the supersolid phase. In this phase the low lying phonon branch mostly depletes with increasing temperature, although there exist a region between

limit [000] and around [100]
. In comparison with the superfluid dispersion relation the excitation spectrum changes its form/shape rather than scaling down with increasing temperature as in the superfluid phase. We observe that the supersolid phase exhibits a more complex and diverse structure than the superfluid or normal solid phases alone.

On Finite Temperature Properties
In this section we will discuss finite temperature properties of the anisotropic Heisenberg model and its solution in the random-phase approximation. In order to be able to compare the temperature dependence of the model in the random-phase approximation with the classical mean-field approximation we chose the set of parameters that was extensively scrutinized by Liu and Fisher [12]: As mentioned in the previous section, Liu and Fisher [12] have chosen this set of parameters because it provides arguably the best fit to the phase diagram of Helium-4. Since we believe that the validity of the quantum lattice gas model is too limited to appropriately reproduce the behavior of Helium-4 over the whole range of temperature and pressure we do not discuss most properties in the pressure-temperature space but rather present the major part of the results in the more comprehensible magnetic field -temperature coordinates. Only where the theory can be compared to relevant experimental data, such as the heat capacity at constant pressure we work in the corresponding representation. The phase diagram of the anisotropic Heisenberg model in T and h z coordinates is given in Figure 5. The dashed lines correspond to the phase diagram of the mean-field approximation. We see that the diagrams are quantitatively quite similar but in the mean-field solution the temperature is somewhat overestimated giving an approximately 30% higher temperature for the tetra-critical point. As mentioned before the critical magnetic fields h z c , due to quantum fluctuations, are lower in the random-phase approximation. This effect is most distinctive on the supersolid to superfluid transition line as there the deletion of the spin magnitude is strongly pronounced. The net vacancy density ǫ, density of vacancies minus density of interstitials sparked interest as the question arose [18] as to whether the number of vacancies follow predictions of thermal activation theory or are due to the effects of an incommensurate crystal. Figure 6 shows the net vacancy density in the supersolid and the normal solid phase at constant pressures. In isobar curves, curves of constant pressure, the magnetic field h z is controlled through Equation (27). We see that the net vacancy density is nearly constant in the supersolid phase. Only in the normal solid phase the net vacancy density increases exponentially with increasing temperature in total agreement with classical thermal activation theory. The almost temperature independent behavior of the net vacancy density in the supersolid phase is an important finding of the quantum lattice gas model and should be observable in high resolution experiments if this effect is real. In Figure (7) we plotted the net vacancy density as a function of h z (or equivalently µ the chemical potential) across the normal solid, the supersolid and the superfluid phases for various temperatures, namely T = 0K, T = 0.4K, T = 0.7K and T = 1.0K. The term net vacancy density does not   have a physical meaning in the superfluid phase and here the quantity ǫ rather refers to the particle density of the fluid. The particle density in the superfluid phase is linear in h z and independent of the temperature T , which follows immediately from condition Eq. (9) as S x A / S x B is equal to one in the fluid phase. In the supersolid phase the dependence of the net vacancy density on the chemical potential is stronger than in the superfluid phase. The chemical potential (magnetic field h z ) is roughly inversely proportional to the pressure. Meaning that the superfluid phase exhibits a higher compressibility than the supersolid phase, which is a quite remarkable result. Interestingly the net vacancy density in the supersolid phase also increases linearly with the chemical potential. This is due to the ratio of the superfluid order parameter S x A / S x B varies as the square root of the magnetic field in the vicinity of the transition: The exponent 1 2 is typical for mean-field type approximations and appears close to all transition lines. Figure 11 shows the net vacancy density of the model on the supersolid to normal solid transition line as a function of temperature, and also reveals the mean-field type square root law dependence. At zero temperature in the solid phase the net vacancy density is equal to zero, hence the crystal is incompressible. In real systems compressibility usually occurs as a result of change in the lattice constant a, characterized by the Grueneisen parameter. The quantum lattice gas model does not take this effect into account since the lattice constant is treated as a constant. However measurements have shown that the lattice constant in the (supersolid) helium is almost a constant indicating that the net vacancy density is the crucial parameter. Also of interest are the free energy and the entropy. Figure (8) shows the free energy of the anisotropic Heisenberg model in the supersolid (CAF) phase at constant pressure. At low T the free energy usually follows a power law: The coefficient α is a universal property of the model which is constant over the entire regime of a phase. The exponent is most easily acquired in the log-log plot, shown in the inset of the plot. In this logarithmic scale α is given by the slope of the curve and the leading contribution is given by, approximately in the supersolid region. This is close in value to the usual T 4 -term attributed to the linear phonon modes, and follows from Equation (25) with ω(k) linear in k. Here the T 4 -term is solely due to the superfluid mode and in real solids an additional T 4 -term contribution, accounting for the lattice phonon-modes, will appear. The logarithmic plot also reveals the leading correction to the free energy given by α = 5.
The non-configurational entropy of the system over the entire range of temperature and magnetic field h z is given in Figure 9. All four phases are visible and as expected from a thermodynamically equilibrated system the entropy is monotonically increasing with respect to temperature. Figure 10 depicts the configurational specific heat at constant pressure in the supersolid and the normal solid phase. The jump in the specific heat at the critical temperature T c , in agreement with the second order phase transition, appears to be smeared out due to numerical inaccuracies as the specific heat is the second derivation of the free energy which had to be integrated of the interval [h z , ∞]. The jump in the specific heat may also be calculated from following formula which is an analogy to the Clausius-Clapeyron equation: As we have C P = C h − T ∂ 2 F ∂T ∂hz ( ∂hz ∂T ) P , we have for the specific heat at constant pressure: For the values corresponding to Figure 10 we obtain an estimated jump of 0.02 which is in good agreement with the curve.

First Order Boundary Lines
In parameter regimes where the supersolid phase does not appear in certain temperature regions there consequently appears a first order phase transition between the superfluid and the normal solid phase. Liu and Fisher [12] compared the free energies of the competing phases to establish the transition line. This procedure is not applicable in the random-phase approximation, as was outlined in the section on thermodynamic properties. Other than the mean-field approximation where the Hamiltonian is given by Equation (6) the effective Hamiltonian of the randomphase approximation is not known. Therefore we have to integrate the first order transition line from a Clausius Clapeyron like equation as derived in the section on Boundary lines. A set of parameters which exhibits such a first order transition at low temperatures is given by: The corresponding phase diagram is shown in Figure 12.
According to mean-field Eq. (14) and Eq. (16) the second order superfluid to supersolid transition as well as the supersolid to normal solid transition is at absolute zero at h z = 0.86603, implying that the supersolid phase does not exist at zero temperature. At higher temperatures (above T > 0.323K) the supersolid phase does exist. At T = 0.323K where the SF-SS and the SS-NS transition lines intersect there occurs a tricritical point. On this tricritical point the superfluid and normal solid phases coexist (as well as the supersolid phase) and the first order phase transition line can be calculated from the ordinary differential equation as has been derived in the section on boundary lines (Eq. (34)). The tricritical point is as such the starting point for the integration of the ODE. There also exist a regime of parameters where the supersolid phase is suppressed at higher temperatures as can be seen in Figure 13 which corresponds to: The SF-SS and the SS-NS transition lines converge before the NF-SF line is reached, hence no tetracritical point such as in Figure 5 is present. The two tricritical points are connected by a first order superfluid to normal solid transition line.

Beyond the Model
In Section 8 we have seen that the supersolid phase appears if and only if the roton dip at γ 1 = −1 and γ 2 = 1, i.e.
[001] of the first Brillouin zone collapses. Additionally we have seen that the model also allows for a collapsing minimum at [111] if condition Eq. (37) is met. A set of parameters that fulfills this condition is given by: Note that the nearest neighbor and the next nearest neighbor constants in this configuration J ⊤ 1 and J ⊤ 2 , corresponding to the kinetic energy, are relatively weak and  are both of equal strength, leading to a highly anisotropic kinetic energy. Figure 14 shows the corresponding phase diagram according to the random-phase approximation. The normal fluid to superfluid transition line starts at at absolute zero and decreases with increasing temperature. According to Eq. (14) and Eq. (16) the critical external fields h z defining the superfluid to supersolid and the supersolid to normal solid transitions are, due to negative values under the square root, imaginary and hence physically not relevant. Consequently in the classical meanfield approximation the superfluid phase extends down to h z = 0 and a first order superfluid to normal solid transition does not occur as the relatively large negative J 2 increases the free energy of a possible solid phase. The random-phase approximation however draws a slightly different picture. Analogous to the classical mean-field solution the random-phase approximation also yields a phase transition near h z = 4.5. But unlike the classical mean-fields solution, the superfluid phase here does not survive all the way down to h z = 0. Due to the particular choice of parameters the superfluid phase becomes unstable at around h z = 2; i.e the quasi-particle spectrum turns imaginary at γ 1 (K) = 0 and γ 2 (k) = −1 ([111]). Interestingly beyond this line no other stable phase exists in the present approach; there is no set of spin fields S xA , S xB , S zA and S zB that solves the self-consistency equations (20) or (21).
To understand why the present approach breaks down and how the phase below h z = 2 might look like we first investigate the physical meaning of the collapsing roton minimum at [  This density wave takes on value one on sub-lattice A and minus one on sub-lattice B, hence it reproduces the periodicity of the supersolid crystal.
In the same way a collapsing roton dip on the main diagonals [111] as given here refers to following density wave: cos(πx/a) cos(πy/a) cos(πz/a) This density wave yields zero on sub-lattice B and alternatingly one and minus one on sub-lattice A (neighbors have opposite signs). Consequently this phase refers to a (super)-solid phase exhibiting three sub-lattices A, B and C where the mean-fields take on three different values (see Fig. 15). It would be quite interesting to study the possibilities and properties of such a phase and we leave as future work the extension of the present approach to account for a three sub-lattice phase and the investigation its properties.

Conclusion
In this paper we extended the mean-field theory of the QLG model by Liu and Fisher [12] at finite temperatures and employed a random-phase approximation, where we derived Green's functions using the method of equation of motion. We applied the cumulant decoupling procedure to split up emergent third order Green's functions in the EoM. In comparison to the MF theory employed by Liu and Fisher [12] the RPA is a fully quantum mechanically solution of the QLG model and therefore takes quantum fluctuations as well as quasi-particle excitations into account. The computational results show that these quasiparticle excitations are capable to render the superfluid phase unstable and thus evoke a phase transition. Quantum fluctuations account for additional vacancies and interstitials even at zero temperature and most interestingly the net vacancy density is altered in the supersolid phase. A potential shortcoming of the RPA is the lack of knowledge of the effective Hamiltonian and therefore the internal and free energy. Usually the free energy is needed to compute first order transition lines as the ground state is given by the lowest energy state. We bypassed this obstacle by deriving a Clausius-Clapeyron like equation which defines the first order superfluid to normal solid transition line. The entropy which is an input parameter of this equation is calculated from the spin wave excitation spectrum. The jump in the specific heat across the second order superfluid to supersolid transition line reveals information about the nature of the transition. However the specific heat is the second derivative of the free energy and thus the jump is smeared out by the numerical calculations. Consequently we derived an equation which gives an optional estimation of the jump and is in good agreement with the numerical estimate. Most important our theory predicts a net vacancy density which, in the supersolid phase is significantly different from thermal activation theory. In the normal solid phase the net vacancy density roughly follows the predictions of thermal activation theory, although quantum mechanical effects give a measurable contribution. Across the phase transition, however, in the supersolid phase the net vacancy density stays rather constant as T increases.