Classical and Quantum Gases on a Semiregular Mesh

The main objective of a statistical mechanical calculation is drawing the phase diagram of a many-body system. In this respect, discrete systems offer the clear advantage over continuum systems of an easier enumeration of microstates, though at the cost of added abstraction. With this in mind, we examine a system of particles living on the vertices of the (biscribed) pentakis dodecahedron, using different couplings for first and second neighbor particles to induce a competition between icosahedral and dodecahedral orders. After working out the phases of the model at zero temperature, we carry out Metropolis Monte Carlo simulations at finite temperature, highlighting the existence of smooth transitions between distinct"phases", The sharpest of these crossovers are characterized by hysteretic behavior near zero temperature, which reveals a bottleneck issue for Metropolis dynamics in state space. Next, we introduce the quantum (Bose-Hubbard) counterpart of the previous model and calculate its phase diagram at zero and finite temperatures using the decoupling approximation. We thus uncover, in addition to Mott insulating"solids", also the existence of supersolid"phases"which progressively shrink as the system is heated up. We argue that a quantum system of the kind described here can be realized with programmable holographic optical tweezers.


I. INTRODUCTION
Investigating the behavior of a many-particle system has an undeniable charm: despite microscopic interactions are undirected, various forms of self-organization ("order") can develop at the macroscale. In the last century, countless examples of emergent order have been described, each with its own practical realization, and many more can be devised by exploring through theory physical situations that are somehow atypical. These indications can stimulate new experimental work or simply be aimed to clarify and expand the scope of the theory itself.
A way to produce novel, unconventional phase behaviors is to consider many-body systems under geometric constraints, since local interactions are frustrated and unusual ground states then appear. A classic example is a (finite) system of hard particles confined in the surface of a sphere [1]. The sphere topology forces an excess of fivefold coordinated particles over sevenfold ones, leading to high-density packings with defects [2][3][4][5][6][7][8]. We note that bosonic atoms confined in thin spherical shells [9] have already been realized [10,11] and are currently studied in microgravity [12,13]. In other cases, frustration is directly embodied in the interaction law -like in spin glasses or in antiferromagnets on a triangular lattice [14,15].
In this paper, we consider a discrete system of particles ("lattice gas") on a spherical mesh of points, which is chosen such that a rich interplay arises between distinct "phases" having the symmetries of a Platonic solid. Clearly, on a finite mesh well-definite phases only exist at zero temperature (T = 0), since for T > 0 any phase transition will be smeared out, i.e., replaced by a smooth crossover region. Using a finite mesh, we greatly reduce the computational effort without however making the phase behavior trivial. If a toy model of classical particles on a finite mesh may look somewhat artificial and hardly corresponds to a real-world system, its quantum counterpart might be different.
The last decades have witnessed a considerable progress in the manipulation of quantum atoms at low temperature, opening the way to a systematic study of correlation effects in many-body systems [16][17][18][19]. While optical lattices [20] are routinely employed in numerous laboratories worldwide as a tool for confinement of quantum atoms, in the last few years a laser technology has been invented, based on the use of optical tweezers [21,22], which allows virtually any type of structure (not necessarily a lattice) to be realized with cold atoms. We are thus encouraged to consider the quantum (Bose-Hubbard) counterpart of the lattice gas on a spherical mesh, with the explicit purpose to compare their thermal behaviors. In particular, we devise a quantum variational theory that predicts supersolid phases and returns the results of a classical mean-field theory when quantum tunneling is precluded.
The rest of the paper is organized as follows. In Sec. 2 we introduce our model and the methods used to investigate its phase behavior. Next, we present our results, first at zero temperature (Sec. 2.1) and then at finite temperatures (Sec. 2.2). In Sec. 3 we deal with the quantum extension of the model in Sec. 2. Using the decoupling approximation, we not only work out the ground-state diagram (Sec. 3.1) but also a few finite-temperature properties (Sec. 3.2). Lastly, we give our concluding remarks in Sec. 4.

II. LATTICE-GAS MODELS ON A SPHERICAL MESH
As anticipated in the Introduction, we hereafter explore the possibility of unusual orderings in a system of particles occupying the nodes of a spherical mesh, chosen to be sufficiently regular that some polyhedral (Platonic) arrangement can occur.
We focus our attention on the pentakis dodecahedron [23] (PD, see Fig. 1), a Catalan solid with 32 vertices obtained by augmenting the dodecahedron with 12 right pyramids on its pentagonal faces, in such a way that the resulting polyhedron is dual to the truncated icosahedron (clearly, the pyramidal apices form the vertices of an icosahedron). Even though the PD is not inscribable, implying that no spherical mesh can be drawn from its vertices, it is straightforward to obtain a biscribed solid by a small distortion of the PD that preserves its connectivity properties and its full icosahedral symmetry. A biscribed solid is any convex polyhedron that has concentric circumscribed and inscribed spheres, where the sphere center is also the centroid of the vertices. The five Platonic solids are biscribed solids, but none of the Archimedean or Catalan solids are. As for the PD, it suffices to adjust the height of the pentagonal pyramids only slightly to force all the vertices to be on the same sphere. The outcome of this construction is the biscribed form of the PD. It is this variant of the PD that is considered hereafter.
The PD has 60 faces (isosceles triangles) and 90 edges (60 short and 30 long). We call PD mesh the skeleton of the PD, i.e., the mesh formed by its edges. The PD mesh is the finite analog of a lattice; the nodes of the mesh (i.e., the PD vertices) are its "sites".
While five edges depart from an icosahedral vertex/site, the number of edges departing from a dodecahedral vertex/site is six (in this sense, icosahedral and dodecahedral sites are "inequivalent"). Setting the circumscribed radius equal to 1, the short-edge length (i.e., the shortest distance in the mesh) is 1 = 30 15 − 15(5 + 2 √ 5) /15 0.64085 . . ., whereas the long-edge length (the second shortest distance in the mesh) The phase behavior of a lattice-gas model on the PD mesh is better studied in the grandcanonical ensemble. Upon increasing the chemical potential µ at fixed T , the mesh becomes increasingly populated, with the possibility of "transitions" between qualitatively distinct arrangements. Denoting c i = 0, 1 the occupation number of site i, we call first (second) neighbors any two sites/particles that are linked by a short (long) edge. According to our nomenclature, an icosahedral site has five first-neighbor sites and no second-neighbor site, whereas a dodecahedral site has three first neighbors and three second neighbors (see Fig. 1).
With these specifications, the grand Hamiltonian of our model reads: where V 1 > 0 (V 2 ) is the coupling between two first (second) neighbor particles. This model can mimic a system of particles adsorbed on a "substrate" sculpted like a PD mesh (think of, e.g., the interstices between atoms in a C 60 molecule), and interacting via a sphericallysymmetric potential with a hard core followed, at larger distances, by a soft short-range repulsion. By suitably tuning the ratio γ = V 2 /V 1 between the couplings, we anticipate the existence of a competition between icosahedral order and dodecahedral order at low temperature, with the further possibility of arrangements with intermediate order.
In the following, we use V 1 as the unit of energy; in turn, this defines a reduced temperature, T * = k B T /V 1 (k B being the Boltzmann constant), and a reduced chemical potential,

A. Zero-temperature phases
At zero temperature, the stable phase at fixed µ is the one minimizing the grand potential Ω. We expect the absolute minimum Ω to be reached in one of a few microstates/configurations, chosen among those exhibiting a homogeneous occupancy of equivalent sites. To be clear, there are five ways to select -out of 20 dodecahedral sites -eight sites forming a cube [24] (similarly, each cube is the union of two tetrahedra). Then, the PD vertices are naturally grouped in three sets of equivalent nodes: we call A the set of icosahedral sites, B any set of cubic sites, and C the set comprising the remaining dodecahedral sites Ω ICO+CUB = 24V 1 − 20µ, and Ω ICO+COC = 36V 1 + 6V 2 − 24µ.
By a lengthy (though elementary) calculation, we arrive at the T = 0 phase diagram depicted in Fig. 2. Here, PQ, PR, and RS are straight lines with equations µ * = 3 + 3γ, µ * = 9 − 6γ, and µ * = 3 + (3/2)γ, respectively; the two half-lines departing from the origin, namely OT and OU, have equations µ * = (15/4)γ and µ * = (3/2)γ, respectively. All the phase boundaries are first-order (since the number of particles changes discontinuously from one phase to the other), except for RT (since both the particle number N and the energy E change continuously through it). Looking at Fig. 2, we note that: i) the ICO phase is only stable for γ > 0, in a region delimited by the µ * = 0 and µ * = 3 lines; ii) the DOD phase is stable in a wide region of the γ-µ plane, bounded from the right by γ = 4/5 and from the above by µ * = 5; iii) ICO+CUB and ICO+COC are each stable in an unbounded set of positive γ values, and coexist along the RS line; iv) for each γ, "empty" and "full" are respectively stable for all values of µ that are sufficiently small or sufficiently large. We monitor a number of equilibrium properties as a function of µ: the number of particles (N = i c i ) and the energy (E = H + µN ), together with their self-and crosscorrelations; the reduced isothermal compressibility, and two specific heats, namely (S being the entropy), expressed in terms of grand-canonical averages through Eqs. (A6) and (A15) of Appendix A. In addition, we also determine a number of order parameters (OPs, see Appendix B for a definition), in order to establish the nature of the system "phase" and the crossover behavior at its boundaries.
For the sake of illustration, take γ = 1/2. For T = 0 the succession of phases is We explore the phase behavior of this model at five temperatures, T * = 0.1, 0.2, . . . 0.5, and in the µ * range from −2 to 7. To account for the possibility of hysteresis, we carry out our simulation runs in sequence: for a given value of µ * , the run is started from the last configuration generated in the previous run at a slightly larger or smaller µ * . Our results are collected in Figs at sufficiently low temperature order can be so robust that we observe hysteresis -the most characteristic feature of a first-order transition. This is the clue to insufficient sampling of the equilibrium distribution, which can be cured by either substantially increasing the length of the MC trajectory or changing the MC algorithm (see more below).
Curiously enough, hysteresis is only found for one of the three transitions present at T = 0. To clarify this point, it is instructive to compute the acceptance probability of the Metropolis move(s) that initially drive the system out of one phase into another. The crucial quantity to look at is the ensuing variation in H, i.e., ∆H = ∆E − µ∆N , with µ given by the transition-point value.
Starting from the empty mesh at µ = 0, the repeated addition of particles in icosahedral sites occurs with probability 1, since at each step ∆E = 0 and ∆N = 1. If we start from ICO, the repeated removal of particles again occurs with probability 1, since at each step ∆E = 0 and ∆N = −1. This means that for µ = 0 the system has equal probability to be in either of the two phases. Hence, no hysteresis would be observed at the empty-ICO transition, as indeed found.
Now consider the ICO-DOD transition at µ * = 15/8. Starting from ICO, the cost to annihilate a particle is ∆H * = 15/8 (since ∆E = 0); the next step of creating a particle in a dodecahedral site has a minimum cost of ∆H * = 31/8 (since ∆E * ≥ 2). Thus, for µ * = 15/8 there is a free-energy barrier for the transition to DOD, and the system then remains for long in ICO notwithstanding DOD is more stable. If we instead start from DOD, we should first annihilate a particle (∆H * = 3/8) and then create another particle in an icosahedral site (∆H * ≥ 17/8). Again, the transition from DOD to ICO is discouraged at low T (though apparently less so than the opposite transition), meaning that for µ * = 15/8 the system is more probably found in DOD than in ICO. As a result, hysteresis will be observed at the ICO-DOD transition.
Lastly, we consider the DOD-full transition at µ * = 5. If we start from DOD, the first step towards "full" is creating a particle in an icosahedral site (∆E * = 5) and the probability for this move is one. If we instead start from "full", the cost for annihilating a particle in an icosahedral site is zero again, since ∆E * = −5. Indeed, no hysteresis is detected at the DOD-full transition.
In Fig. 4 we plot a few response functions for γ = 1/2. At the lowest temperatures all these quantities exhibit a distinct peak near each T = 0 transition point, which is where the energy and particle number are subject to the sharpest variations. In the "empty" phase the reduced isothermal compressibility is close to 1 (the ideal-gas value); upon heating, every asperity in its profile becomes smoothened until ρK T becomes a monotonously decreasing function of µ. As T grows, both specific heats develop a broad maximum inside the ICO and DOD regions. Admittedly, these are the locations where, in the moderately hot system, the fluctuations of energy and particle number are stronger. Eventually, both maxima gradually deflate, the DOD bump being the last to disappear. Finally notice that inside the "empty" region the two specific heats have different behaviors at low temperature: while the covariances involving energy are practically zero, the mean square fluctuation of N is of order N (see Eq. (2)). Looking at Eqs. (A6) and (A15), we thus conclude that C µ /k B ≈ (βµ) 2 1 and C N ≈ 0.
As a second example, consider γ = 1 and T * = 0.05 close to µ * = 9/2 -which is where, at T = 0, ICO+CUB transforms into ICO+COC. By the same argument put forward before we would conclude that this transition is accompanied by hysteresis at low temperature.
However, when plotting N as a function of µ, in a small neighborhood of µ * = 9/2 we see a narrow plateau in the middle between the N levels in the two phases (top-left panel of all icosahedral sites occupied, thus confirming that ICO+CUB (ICO+COC) is the stable phase for µ * < 9/2 (µ * > 9/2); instead, exactly for µ * = 9/2 we have counted as many as 240 distinct microstates, all with N = 22 and E * = 33, having the same grand potential as ICO+CUB and ICO+COC. It is the existence of such configurations that makes the transition between ICO+CUB and ICO+COC smoother than expected. We care to stress that this occurrence is rather exceptional; usually, phases are well separated in free energy and compete with each other only in pairs.
As anticipated, hysteresis is an annoying problem due to the inadequacy of Metropolis dynamics to overcome free-energy barriers. To solve this issue we have employed the Wang-Landau algorithm [25], which directly computes the density of states and is thus particularly suited for a free-energy landscape with multiple minima. Compared to the original algorithm, the refinement parameter ln f was reduced at a slower rate (at regular intervals, ln f is divided by 1.1 rather than 2), which greatly reduces the (already small) saturation error [26]. We report results for γ = 1/2 in Figs Fig. 7 the density of states g as a function of N and U (i.e., the sum of the first two terms on the r.h.s. of (1)) and the probability density cycles performed in each run, which is 5 × 10 6 and 5 × 10 7 , respectively. Bottom panels: OPs for ICO+CUB and ICO+COC, for the simulation with 5 × 10 7 equilibrium cycles per run. We see a narrow interval of µ * values around 9/2 where the order is neither ICO+CUB nor ICO+COC.
of the particle number, P (N ). In the ICO region, well before the transition at µ * = 15/8, a second peak builds up in P (N ), which, as µ is increased, is gradually shifted to larger and larger N values until becoming centered at N = 20.
Finally, it is useful to compare MC results with the outcome of a mean-field (MF) theory.
The simplest approach is to estimate the grand potential of (1) using the Gibbs-Bogoliubov (GB) inequality with a trial probability density π[c] given as an uncorrelated product of one-site terms: In the previous equation, with 0 ≤ ρ A ≤ 1, and similarly for B and C. The rationale behind Eq. (6) is that the average occupancy takes a possibly different value in each set of equivalent nodes, being ρ A for the icosahedral set, ρ B for the cubic set, and ρ C for the co-cubic set.
An upper bound to the exact grand potential Ω is the GB grand potential Ω * , where Then, it is a simple matter to show that Observe that the value of Ω * in the putative ground states listed in Sec. 2.1 exactly reproduces their respective grand potentials. The stationary values of (9) fulfill the coupled equations , and These equations are solved numerically, seeking the solution that provides the absolute minimum Ω * for the given T and µ.
To have a flavor of how MF theory works, we draw in Fig. 8 the theoretical phase diagram on the µ-T plane for γ = 2/3, corresponding to a lattice gas where the icosahedral phase covers a µ range as wide as that of the dodecahedral phase, see Fig. 2. Owing to a symmetry property of Eqs. (10), the phase diagram is symmetric around µ = 5/2 (see also Sec. 3.2).
At T = 0 the phase boundaries in Fig. 8 are exact. As T grows, the theory predicts a gradual weakening of ICO and DOD orders, as witnessed by the decrease of |ρ A − ρ B | on heating, eventually resulting in an abrupt (first-order) transition to either "empty" (ρ A = ρ B = ρ C < 0.5) or "full" (ρ A = ρ B = ρ C > 0.5), according to whether µ < 5/2 or µ > 5/2. Clearly, this singularity is an artifact of MF theory, since no sharp transition is present in our system for T > 0. Another unphysical prediction of the theory concerns the behavior of the model in a narrow strip of temperatures and chemical potentials near µ * = 0 and µ * = 5. In

HUBBARD MODEL
The PD mesh is regular enough that we can study the quantum analog of the lattice-gas model in relatively simple terms. The obvious bosonic counterpart of (1) is the hard-core limit of the extended Bose-Hubbard (BH) model where a i , a † i are bosonic field operators and n i = a † i a i is a number operator. Moreover, t ≥ 0 is the hopping amplitude, taken for simplicity to be the same for first-and second-neighbor pairs, whereas U > 0 is the on-site repulsion. In the hard-core limit U → +∞, the site occupancies are effectively restricted to zero or one and the U term can thus be discarded.
For hard-core bosons, creation and annihilation operators at different sites commute, whereas a i and a † i are anticommuting operators as a result of the dynamical suppression of Fock states with two or more particles in the same site [27].
In the original BH model [28][29][30], where the second, fourth, and fifth terms in H are absent, the tunneling term (kinetic energy) is minimized by a condensed state spread over the entire volume of the system, whereas the potential energy favors particle localization. As a result, the T = 0 system exists as either a superfluid (large t/U ) or a Mott insulating ground state (small t/U ), separated by a quantum transition. The Bose-Hubbard Hamiltonian can be derived starting from the second-quantized Hamiltonian describing a gas of ultracold bosonic atoms subject to an optical-lattice potential [18].
The phase scenario becomes richer when the range of interaction is increased: depending on the lattice, other Mott insulating ground states (density waves) may appear; moreover, crystalline order may coexist with superfluidity (supersolids) [31][32][33][34][35][36]. Supersolidity is a fascinating property of quantum matter, which has only recently been experimentally detected in a gas of dipolar atoms [37][38][39]. In a supersolid, atoms can simultaneously support frictionless flow and form a crystal. As suggested by Leggett, a rotating supersolid should have a moment of inertia that is reduced with respect to its classical value [40]. This phenomenon is called "nonclassical rotational inertia" and its first observation is reported in a paper published this year [41].
Due to the semiregular character of the PD mesh, in our system the superfluid phase would be discouraged in favor of less-symmetric condensed phases, and a supersolid region will then occur at low temperature. To check this expectation, we employ a mean-field theory, an approach known to give accurate results in the continuum [42][43][44].

A. Decoupling approximation
As done in Refs. [24,45], we analyze the phase diagram of the extended BH model in the hard-core limit using the decoupling approximation (DA) [46,47]. The latter approach consists in linearizing the hopping and interaction terms in (11) as where the thermal averages a i ≡ φ i and n i ≡ ρ i are determined self-consistently; φ i and ρ i represent the superfluid OP and the average occupancy in the i-th site, respectively (the condensed fraction is |φ i | 2 , see e.g. [45]). The DA Hamiltonian is a sum of one-site terms, given by the first and second neighbors of i, respectively). While referring to [45] for a full justification of DA, it is worth to underline that the self-consistency equations for φ i and ρ i are also the conditions under which the grand potential of (13) is stationary. If more stationary solutions are found, we must select the one providing the minimum grand potential.
As discussed before, the 32 nodes of the PD mesh are naturally classified as icosahedral (A), cubic (B), or co-cubic (C), implying that the number of variational parameters in (13) is reduced to six. A phase with φ A = φ B = φ C = 0 is a Mott insulator, whereas a homogeneous occupancy together with φ A = φ B = φ C = 0 defines a superfluid. Any unbalance between φ A , φ B , and φ C corresponds to a supersolid.
Looking at Fig. 1 we soon realize that in such a way that the DA Hamiltonian becomes 16) and In the hard-core limit, the eigenvalues of each partial Hamiltonian in (16) follow from the diagonalization of a 2 × 2 matrix. We easily obtain: Therefore, for T = 0 the grand potential of (13) reads Ω = 12λ with E 0 = 12E For T > 0, the partition function of (13) reads: yielding the grand potential In seeking the stationary solutions of (22), it can be assumed -without loss of generality -that φ A , φ B , φ C are real and positive [24,45]. Putting the derivative of (22) with respect to each free parameter equal to zero, and suitably rearranging the formulae, we arrive at the coupled equations: The above equations can be solved numerically by, e.g, the method described in [24]. Upon combining the six equations (23) together, the following identities are easily derived: indicating that the value of each superfluid OP is comprised between 0 and 1/2.
Before moving to numerical results, we show that for t = 0 the DA theory is equivalent to the MF theory described in Section 2.2. Indeed, for t = 0 the grand potential (22) becomes with Upon differentiating (25) with respect to each density parameter and putting the result equal to zero, the same equations (10) are eventually obtained. If these equations are substituted back into (25), then the grand potential (9) is obtained, indicating that the DA phase diagram of the t = 0 quantum system is exactly identical to the phase diagram of the lattice-gas system in the MF approximation.

B. Numerical results
Until now, the value of V 2 was arbitrary. For the sake of example, we henceforth take γ = V 2 /V 1 = 2/3. We note that the case γ = 1 was considered in [45].
We have first solved Eqs. (23) numerically for T = 0 and various µ values, being careful that the minimum Ω solution is picked out in each case. The resulting phase diagram is shown in Fig. 9. In addition to the "empty" phase (ρ A = ρ B = ρ C = 0) and the "full" phase (ρ A = ρ B = ρ C = 1), we observe an icosahedral phase (ρ A = 1, ρ B = ρ C = 0) and a dodecahedral phase (ρ A = 0, ρ B = ρ C = 1). All these phases are insulating (φ A = φ B = φ C = 0) and incompressible (i.e., the density is constant throughout the phase). Upon increasing t at fixed µ a supersolid phase invariably appears, characterized by ρ A = ρ B = ρ C . We also note that our phase diagram is symmetric around µ = 5/2. Indeed, we see from Eqs. (23) that ρ A,B (µ) = 1 − ρ A,B (5 − µ) and φ A,B (µ) = φ A,B (5 − µ), and an analogous symmetry property holds for (22). It is worth stressing the similarities and differences between Fig. 9 and the phase diagram of hard-core bosons on a triangular lattice [33,47,48]: The overall structure is the same, but the nature of the condensed phase at large t is different, being herein supersolid rather than superfluid -owing to the frustration effect associated with the existence of inequivalent nodes. We actually distinguish four different supersolid phases (see Fig. 10, where the OPs are plotted for two representative values of t). While φ B is always slightly larger than φ A , in the region between the ICO and DOD lobes (where, in particular, t < 0.266) we find ρ A > ρ B for µ < 5/2 (supersolid 1b) and ρ A < ρ B for µ > 5/2 (supersolid 2b). Outside the lobes, we instead find ρ A < ρ B for µ < 5/2 (supersolid 1) and ρ A > ρ B for µ > 5/2 (supersolid 2), namely the densities are in reverse order with respect to the reference "solid" phase. For t 0.278, the values of ρ A and ρ B jump discontinuously at µ = 5/2, signaling that the phase transitions along this line are first-order. A further first-order line runs vertically near t = 0.266, separating the supersolid phases 1b and 2b from the supersolid phases 1 and 2, respectively. Finally, there are four second-order transition lines: the two lines separating "empty" and "full" from the adjacent supersolid, the descending part of the boundary between ICO and supersolid 1b, and the ascending part of the boundary between supersolid 2b and DOD.
Imposing B-C symmetry, we may simplify Eqs. (23) and then determine the equations for the continuous-transition loci, following the same procedure as illustrated in Ref. [45]. We eventually find: t ("empty"-supersolid boundary) ; In fact, looking at Fig. 9 we see that the lower branch of the ICO-supersolid locus is only virtual, since this transition is preempted by a first-order phase transition. A similar comment applies for the upper branch of the DOD-supersolid locus.
For T > 0 the phase diagram evolves in the way illustrated in Fig. 11. Clearly, the indications of DA for non-zero temperatures are less accurate; moreover, the prediction of sharp phase boundaries is an artifact of the approximation, the transitions being actually smooth crossovers. Already for T * = 0.5 we observe a retreat of every supersolid phase with respect to T = 0. The smallest t value for which the system can be supersolid is now slightly larger than 0.20. ICO and DOD too lose ground in favor of "empty" and "full", respectively, a trend that will become more marked on increasing T further. An effect of finite temperature is that the occupation unbalance between A and B/C is no longer sharp in ICO and DOD, and is moreover µ-dependent; but, similarly to T = 0, the site occupancies are independent of t. Notice that the ICO and DOD densities as well as the ranges of stability are exactly the same as predicted by the MF theory of Sec. 2.2 (see Fig. 8). For any T > 0 the densities are µ-dependent and independent of t also in the former "empty" and "full" phases, though still homogeneous throughout the mesh. For T * = 0.7 the supersolid 1b and 2b phases are nearly disappeared; furthermore, the µ extent of ICO and DOD is more than halved with respect to T = 0. Finally, in the phase diagram for T * = 1 not only 1b and 2b but also ICO and DOD are no longer present. Moreover, the supersolid sector lies to the right of t = 0.35; for t < 0.35, the occupancies evolve continuously through µ * = 5/2.
In spite of the elementary character of the DA, the main traits of the thermal evolution outlined above are correct, being in line with other theoretical studies [49][50][51] and simulations [52][53][54]. To recapitulate, the quantum phases are weakened by the thermal fluctuations associated with finite temperatures. Thus, for T > 0 a normal fluid appears in the system.
Here, the superfluid OP is zero and the density at each site of the mesh becomes non-integer. This is to be contrasted with the incompressibility of insulating quantum phases, which have integer occupancy at each site. Thermal fluctuations also undermine the supersolid phases, which are shifted towards higher and higher hopping amplitudes as T is progressively increased.

IV. CONCLUSIONS
Using a combination of analytic calculations and numerical experiments we have worked out the phase behavior of a lattice-gas model on the skeleton of the PD. Only particles connected by a PD edge are allowed to mutually interact, with different couplings for short and long edges. Depending on the ratio γ between these couplings, various ordered phases are observed at T = 0 (in addition to "empty" and "full"): icosahedral, dodecahedral, icosahedral+cubic, and icosahedral+co-cubic. For T > 0 we study the phase behavior of the lattice gas by MC simulation. The total number of sites (32) is sufficiently small that the system is quickly equilibrated at any temperature, with negligible uncertainties on the thermal properties. Despite the absence of sharp phase transitions in a finite system, at low temperature we observe strong hysteresis at some of these boundaries. We have shown that hysteresis, which would occur with any MC algorithm with local updates, can be cured by making an entropic sampling. In this respect, it would be intriguing to examine whether anything similar to the concept of nucleation barrier [55,56] applies for this model (but we leave this for future work).
A variation on the theme of the present model is one where the occupancy of sites is unrestricted. In this case, for high chemical potentials we expect to observe the formation of cluster phases, as in [7]. A mean-field theory similar to that formulated in [57] would probably suffice to obtain accurate predictions for the phase behavior of this system, thus making the simulation unnecessary.
Finally, we have considered the quantum analog of the lattice-gas model and solved it using the decoupling approximation. This theory predicts various Mott insulating ground states, each being the counterpart of a phase of the classical model, as well as a number of supersolids for higher hopping amplitudes. Admittedly, it is the frustration effect due to the semiregular character of the PD mesh that causes the superfluid phase to be superseded by a supersolid phase. The take-home message is that confining bosonic particles in a semiregular mesh is an easy way to stabilize a supersolid in an ultracold quantum gas. Upon heating, the extent of all the T = 0 phases gets progressively reduced, leaving room to normal-fluid behavior.

Acknowledgments
We express our gratitude to one of the Referee for suggesting entropic sampling as a cure to the hysteresis apparent in Figs. 3 and 4.
Appendix A: Calculation of specific heats In this Appendix, we derive the statistical-mechanical formulae for the specific heats in Eq. (3).
Working in the grand-canonical ensemble, the partition function (a sum over microstates) with T, V , and µ as control parameters (V is the system volume). The grand potential Ω, that is the thermodynamic potential in the T, V, µ representation, is given by Ω = −k B T ln Ξ.
Next, we focus on a different specific heat, calculated by keeping N (and V ) fixed: To enforce N = const., the chemical potential must be in a suitable relation with T and V : N (T, V, µ) = const. =⇒ µ = µ(T, V ) . Hence and In turn, the µ derivative of S is given by a Maxwell relation: Putting Eqs. (A7)-(A11) together, we end up with: Now, using (A3) we obtain: ∂N ∂µ T,V = β δN 2 (A13) and Finally, substituting Eqs. (A6), (A13), and (A14) into (A12), we arrive at µ. Clearly, there is no unique way to define these OPs -our proposal below is just one possibility. Measuring the degree of ICO+CUB order is more subtle, since in the phase region where this "phase" is stable the simulated system circulates, even at the lowest temperatures, between five different basins of microstates. As a result, in a long simulation the average occupancy will be the same for all dodecahedral sites. Here, the crucial observation is that the cosine of the angle θ ij formed by the vector radii relative to two distinct cubic sites i and