Ultracold Bosons on a Regular Spherical Mesh

Here, the zero-temperature phase behavior of bosonic particles living on the nodes of a regular spherical mesh (“Platonic mesh”) and interacting through an extended Bose-Hubbard Hamiltonian has been studied. Only the hard-core version of the model for two instances of Platonic mesh is considered here. Using the mean-field decoupling approximation, it is shown that the system may exist in various ground states, which can be regarded as analogs of gas, solid, supersolid, and superfluid. For one mesh, by comparing the theoretical results with the outcome of numerical diagonalization, I manage to uncover the signatures of diagonal and off-diagonal spatial orders in a finite quantum system.


Introduction
Gases of ultracold bosonic atoms loaded in an optical lattice provide the unique opportunity to study quantum many-body effects under controlled conditions [1,2]. To a very good approximation, the atoms can be described by a Bose-Hubbard (BH) Hamiltonian [3] whose parameters can be tuned by laser light [4,5]. By changing the configuration of the lasers, many lattice geometries can be explored, making optical lattices a powerful and versatile tool.
In a system at zero temperature (T = 0), thermal fluctuations are frozen out and quantum fluctuations prevail. These microscopic fluctuations can induce phase transitions in the ground state of a many-body system, driven by a non-thermal control parameter, such as chemical potential, magnetic field, or chemical composition. As a concrete example, consider a dilute gas of bosons at temperatures low enough that a Bose-Einstein condensate is formed. The condensate is described by a wave function consisting of every particle in the same state spread over the entire volume of the system, and typically exhibits superfluidity. An interesting situation appears when the condensate is subject to a lattice potential in which particles can move from one lattice site to the next only by tunneling. If the lattice potential is increased smoothly, the system remains in the condensed phase as long as the repulsion between atoms is small compared to the tunnel coupling (assuming, by the way, that the range of the repulsion is much smaller than the lattice spacing). In this regime, where the tunneling term dominates the Hamiltonian, a delocalized wave function still minimizes the total energy of the many-body system. When the strength of the repulsion becomes large compared to the tunnel coupling, the total energy is made minimum when each lattice site is occupied by the same number of atoms; as a result, phase coherence is lost and the system becomes insulating. The addition of a longer-range repulsion will make the phase behavior richer, with the possibility of a non-superfluid density wave and a supersolid ground state where crystalline order coexists with superfluid behavior (see, e.g., Ref. [6]). Experimentally, a way to prepare ultracold gases of long-range interacting bosons is to use atoms (such as chromium [7] or dysprosium [8]) and molecules having a large magnetic or electric dipole moment. The usual BH model predicts a T = 0 phase transition from a superfluid phase to a Mott insulator phase as the ratio of the hopping matrix element between adjacent sites (t, in absolute terms) to the on-site interaction (U) is reduced. The overall number density of particles is controlled via a chemical-potential parameter µ. As µ grows at fixed t and U, the lattice becomes increasingly filled with particles, but this can only occur outside the Mott regions since the insulator phase is incompressible [3]. The BH model has been studied in many lattice geometries and with several techniques (mean-field theory [9][10][11], perturbation theory [12][13][14][15][16], and quantum Monte Carlo simulation [17,18], to name but the most commonly employed). When a further repulsion V is introduced between nearest-neighbor atoms ("extended BH model"), new phases may arise, in primis a supersolid phase [19][20][21][22][23][24][25][26][27]. Another variant of the BH model is hard-core bosons, where site occupancy is restricted to zero or one, corresponding to taking the U → ∞ limit [6,[28][29][30][31][32].
I hereafter present the results of yet another investigation of the extended BH model, now choosing a finite graph as hosting space for bosons. Even though clearcut phase transitions (i.e., thermodynamic singularities) cannot occur in a few-particle system, a convenient choice of boundary conditions may alleviate the difference with an infinite system, making the study of a finite quantum system valuable anyway. A practical solution is to use spherical boundary conditions (SBCs), which have often been exploited in the past to discourage long-range triangular ordering at high density [33][34][35][36][37][38]. On the other hand, SBCs make it possible to observe novel forms of ordering, viz. into regular polyhedral structures, that are simply unknown to Euclidean space-see Refs. [39][40][41]. An added value of a spherical mesh of points is the possibility to vary the site coordination while keeping the overall geometry strictly two-dimensional. Bosons confined in spherical (bubble) traps have been produced experimentally [42,43] and are going to be studied soon under microgravity conditions [44,45]. In the present study, the extended BH model is considered on a finite mesh of points homogeneously distributed on the unit sphere, i.e., coincident with the vertices of a regular polyhedron inscribed in the sphere (we may call it a Platonic mesh). Despite consisting of a finite number of nodes, a regular spherical mesh shares an important feature in common with an ordinary lattice, namely all sites are equivalent; for this reason, the phase behavior of particles living on a Platonic mesh would not deviate much from that of an infinite system (this intuition will be checked in a particular case).
The plan of this paper is the following. After introducing the models in Section 2, the choice of the underlying mesh is discussed in detail, giving priority to those regular grids where a subset of nodes form a mesh that is also regular. Then, in Section 3, a mean-field (MF) analysis of the ground state is carried out for all the models considered; in one case, the indications of theory are validated against the results of exact diagonalization (Section 4). From this comparison, we find the artifacts of the MF approximation as applied to a finite system. Finally, some concluding remarks are given in Section 5.

Classical Models
To illustrate the main idea, let it initially be considered a system of classical point particles defined on the sites of a cubic mesh stretched over the surface of the unit sphere. Each of the eight nodes of the mesh has three nearest neighbors (NN, at chord distance r = 2/ √ 3) and three next-nearest neighbors (NNN, at chord distance 2 √ 2/3). For simplicity, assume that the occupancy n i of site i can only be 0 or 1 (i = 1, . . . , 8). Finally, choose the system Hamiltonian to be H[n] = V ∑ i,j n i n j with V > 0 (each NN pair in the sum is counted only once). The grand potential Ω(T, µ) of this system is where β is the inverse temperature. At T = 0, the formula is simpler: Any of the points of absolute minimum in (2) represent the actual equilibrium state of the system for that µ. For particles living in a continuous space, minimization of H − µN is carried out among a selection of crystalline states that are thought to be relevant based on symmetry considerations (see, e.g., Ref. [46]). In the present case, where the total number of microstates is small (2 8 = 256), the T = 0 grand potential can be determined exactly for each µ by a scrutiny of all possible energies. While the mesh is empty for µ < 0 and completely filled with particles for µ > 3V, in the interval from 0 to 3V only half of the nodes are occupied, and these fall at the vertices of a regular tetrahedron. This twofold-degenerate state with checkerboard order can be viewed as the finite-size counterpart of a crystalline phase. The rationale behind the choice of a cubic mesh is now clear: by introducing a repulsion between occupied NN sites, we promote the occurrence of a Platonic "crystal", i.e., the regular tetrahedron ("CT model", see Figure 1 left). There is only one other possibility to obtain a non-frustrated crystal-like state at T = 0, which is using a regular dodecahedral mesh. We shall see that (i) by discouraging the occupancy of NN sites, a cubic "crystal" is stabilized for sufficiently small µ > 0 ("DC model",  For a system of classical hard-core particles defined on a regular dodecahedral mesh, the total number of microstates is 2 20 = 1,048,576. By direct inspection, we see that the minimum-Ω "phase" for the DC model is the empty mesh for µ < 0, a cube for 0 < µ < 3V/2 (a fivefold-degenerate state, since there are 5 ways to form a cube with the vertices of a regular dodecahedron), and the complement of a cube ("co-cube" in the following) for 3V/2 < µ < 3V; finally, for µ > 3V the mesh is completely filled. For the DT model, the mesh is empty for µ < 0, filled with particles located at the vertices of a regular tetrahedron for 0 < µ < 3V/2 (a tenfold degenerate state), and completely filled for µ > 9V. For µ between 3V/2 and 9V, a different ground state exists for each even number of occupied sites in the range 6 to 16, none of which corresponds to a simple geometric arrangement.
I have plotted in Figure 2 the evolution with temperature of a few properties of the CT and DC models. In addition to the grand potential Ω, the figure shows the total number of occupied sites (N), the total energy (E), and the order parameters for tetrahedral (S t ) and cubic order (S c ). The latter quantities are defined so as to discriminate the Platonic "phase" from the other T = 0 states. A proper order parameter should be insensitive to the orientation of the polyhedron, hence it can only depend on the relative angles between the vertices vector radii departing from the sphere center: for a regular tetrahedron all these angles are equal to α = arccos{−1/3} (with sin α = 2 √ 2/3), while they are π − α, α, and π for a cube. With the idea to penalize configurations that do not match the wanted structure, my choice of the OPs is the following: and where the constants k t and k c are chosen, following the advice in Ref. [47], so that S vanishes for a random distribution of angles: Looking at Figure  It is worth considering how the T = 0 phase diagram of the CT and DC models gets modified when the occupancy of a node is allowed to take any value. For simplicity, as relevant T = 0 states only the "cluster crystals" [48][49][50] originated from the previously identified ground states are considered. Call U > 0 the on-site energy and V the NN repulsion. The grand Hamiltonian now reads H[n] = U/2 ∑ i n i (n i − 1) + V ∑ i,j n i n j − µ ∑ i n i . For the CT model, the grand potential of a tetrahedral "cluster crystal" with n particles per site is Ω (4) n = 2n(n − 1)U − 4nµ, while the grand potential of the "cluster crystal" with n particles per site is Ω For a given µ, the most stable "phase" corresponds to the minimum Ω. Up to µ = 0 the stable "phase" is still the empty mesh. As µ grows further, the site occupancy increases monotonically within each of the families Ω (4) n and Ω (8) n , but the exact sequence of stable "phases" depends on U/V. I show three cases in Figure 3. For U = 5V (right panel), the behavior for small µ > 0 recalls that found in the hard-core limit; however, as µ grows, each site becomes increasingly populated, and a whole sequence of cluster states is found. The opposite occurs for U = 2V (left panel), where only the cluster states with checkerboard order are stabilized for µ > 0. A curious situation occurs for U = 3V, where the competition between different cluster states becomes so stringent that "crystals" of the two families coexist at regular intervals along the µ axis (center panel). As a matter of fact, hybrid ground states with features intermediate between those of the tetrahedral and cubic cluster states may also occur. For example, I have checked for U = 5V that a microstate with two overlapping particles at the vertices of one tetrahedron and one particle at the vertices of the other tetrahedron is the most stable ground state in the µ interval from 8V to 11V.
Moving to the DC model, we have three different families of "cluster phases": the cubic family with grand potential Ω (8) n = 4n(n − 1)U − 8nµ; the co-cubic family with grand potential Ω (12) n = 6n(n − 1)U + 6n 2 V − 12nµ; and the family of those cluster states where every site is filled with the same number of particles, with grand potential Ω (20) n = 10n(n − 1)U + 30n 2 V − 20nµ. Without delving into details, the phase behavior as a function of U is similar to the CT model, with multiple coexistence between all three types of cluster crystals now occurring for U = 2V.

Quantum Models
The previous analysis has served to set the stage for the forthcoming study of the extended BH model on a regular spherical mesh; only the quantum versions of the CT and DC models (denoted QCT and QDC, respectively) are considered below.
The system is defined by the following grand Hamiltonian: In the present investigation, H describes a system of bosons on a spherical mesh of M sites, a i , a † i are bosonic field operators (i = 1, . . . , M), and n i = a † i a i is a number operator. Moreover, t ≥ 0 is the hopping amplitude between NN sites, U > 0 is the on-site repulsion, and V > 0 is the strength of the NN repulsion. In the hard-core limit U → ∞, the site occupancy is restricted to zero or one and the U term in (6) vanishes; it is only this limit that is treated hereafter.
In the BH model on a standard lattice, at T = 0 we observe an insulator-superfluid transition with increasing t for every positive µ. The addition of V may stabilize (depending on the lattice) a density wave at low t, as well as a supersolid phase at the boundary between the insulator and superfluid phases. These features are also present in the hard-core limit, as reported, e.g., in Refs. [30,32].

MF Investigation
MF theory is the method of choice when a new many-body problem is attacked; it has been frequently applied for continuous quantum systems as well, as an effective means to identify the ground states and quantum transitions between them (cf. Refs. [51][52][53]). Various versions of MF theory exist in discrete space; here the so-called decoupling approximation [11,29,32] is employed, which gives the advantage of a fully analytic treatment of the problem. Clearly, the accuracy of MF theory would be questionable when applied for a finite quantum system, even one lacking a boundary surface. This will make urgent an assessment of MF theory against exact results, which however is delayed until the next Section.
In the decoupling approximation, the two-site terms in the Hamiltonian (6) are linearized so that the eigenvalue problem becomes effectively one-site. This is accomplished by first writing the exact identity with δa i = a i − a i (for i = 1, . . . , M), and then neglecting the last term. Similarly, n i n j ≈ n i n j + n i n j − n i n j .
The averages a i ≡ φ i and n i ≡ ρ i (either ground-state averages or thermal averages) are to be determined self-consistently; φ i and ρ i (with 0 ≤ ρ i ≤ 1) represent the superfluid order parameter and the local density for site i, respectively. For hard-core bosons, we readily obtain H ≈ H MF with In Equation (9), F i = ∑ j∈NN i φ j and R i = ∑ j∈NN i ρ j are sums over the nearest neighbors of the i-th site. For a bipartite mesh, sites are either of type A or B, hence the unknown parameters are four, namely The simplification is obvious: rather than working on a Fock space of dimensionality 2 M , for a bipartite mesh the basis states are just four, namely |0, 0 , |0, 1 , |1, 0 , and |1, 1 , corresponding to the possible occupancies of a pair of A and B sites.

QCT Model
For the QCT model, the mesh is bipartite and formed by two tetrahedral sub-meshes with four points each. A point of grid A has three neighbors, all belonging to grid B, and vice versa. Hence, Using a subscript A (B) for operators relative to a single A (B) site, the MF Hamiltonian reads: with On the 4-vector basis B = {|0, 0 , |0, 1 , |1, 0 , |1, 1 }, the Hamiltonian (11) is represented by a 4 × 4 Hermitian matrix: Non-superfluid phases have φ A = φ B = 0. In this case, the matrix becomes diagonal, meaning that the basis vectors are all (grand-)energy eigenstates. From the self-consistency equations ρ A = n A and ρ B = n B , it readily follows that ρ A = ρ B = 0 for |0, 0 (empty mesh), ρ A = 0 and ρ B = 1 for |0, 1 (tetrahedral-B "crystal"), ρ A = 1 and ρ B = 0 for |1, 0 (tetrahedral-A "crystal"), and ρ A = ρ B = 1 for |1, 1 (fully occupied mesh). The eigenvalue gives the grand potential Ω, which is zero for the empty mesh, −4µ for the twofold-degenerate tetrahedral "crystal", and 12V − 8µ for the filled mesh. These grand potentials are exactly the same as in the CT model. Superfluid and supersolid "phases" have non-zero, possibly distinct, complex values of φ A and φ B . For sure, φ A and φ B have equal phases since only the magnitude of the order parameter can be spatially modulated; hence, the arbitrary phase can be taken as zero. Without loss of generality, we may assume φ A and φ B to be positive quantities. We shall see in the next section that the exact eigenstates of H for t > 0 do not distinguish between A and B; hence, our search can be restricted to homogeneous (superfluid) solutions: φ A = φ B = φ and ρ A = ρ B = ρ. We are thus led to the following characteristic equation (with E 0 = 24tφ 2 − 12Vρ 2 ): Doing the simplifications, we arrive at with a = E 0 , b = E 0 + 4(3Vρ − µ), and u = −12tφ. The minimum root of (15) is A real eigenvector of E is or, in explicit terms, Now imposing the conditions we arrive at the two coupled equations which are easily solved to give The above φ solution only exists provided that −3t < µ < 3V + 3t, which is a necessary condition for the existence of the superfluid. Plugging these ρ and φ in Equation (16), we finally obtain the superfluid grand potential: This outcome can also be obtained from the general relation (see Equations (6)- (8)) For the superfluid, we have φ i = φ and ρ i = ρ, with φ and ρ given by Equation (21); hence, By comparing the grand potentials of all the "phases", we arrive at the ground-state diagram in Figure 4. Along the straight lines µ = −3t and µ = 3V + 3t, separating the superfluid from the insulator "phases" at small and large µ, the µ derivative of the grand potential (=− N ) is continuous. Instead, along the line separating the "crystalline" lobe from the superfluid, the average number of occupied sites shows a jump discontinuity. Only at the lobe vertex (V/2, 3V/2) the tetrahedral-superfluid transition is continuous. We will see in Section 4 to what extent the indications of MF theory are accurate, considering that the cubic mesh consists of 8 sites only.

QDC Model
The decoupling approximation for the QDC model works similarly as for the QCT model. Again, the starting point is Equation (9) and the mesh is partitioned into a cube (A, 8 sites) and a co-cube (B, 12 sites). We have: The MF Hamiltonian then reads: with On the B basis, the Hamiltonian (27) is represented by the matrix: For phases with φ A = φ B = 0, the matrix is diagonal and, like in the QCT model, the basis vectors are all (grand-)energy eigenstates. From the self-consistency conditions ρ A = n A and ρ B = n B it follows that ρ A = ρ B = 0 for |0, 0 (empty mesh), ρ A = 0 and ρ B = 1 for |0, 1 (co-cubic "crystal"), ρ A = 1 and ρ B = 0 for |1, 0 (cubic "crystal"), and ρ A = ρ B = 1 for |1, 1 (filled mesh). The grand potential Ω is zero for the vacuum, −8µ for the cubic "crystal", 6V − 12µ for the co-cubic "crystal", and 30V − 20µ for the filled mesh. These values are the same as for the DC model, hence the same sequence of "phases" as a function of µ is observed in the QDC model for t = 0. Moving to phases with φ A , φ B = 0, I first consider the possibility of a superfluid (φ A = φ B = φ > 0 and ρ A = ρ B = ρ). In this case, the characteristic equation takes the form: Observing that a + d = b + c ≡ s, the minimum eigenvalue of H MF turns out to be (see Appendix A): Notice that this energy is exactly 5/2 times that of the QCT model (see Equation (16)). Using p = (b − c)/2 = (d − a)/10 = 2(3Vρ − µ) and q = u/3 = v/2 = −12tφ, the coordinates (x, y, z, w) of an eigenvector |ψ E of E satisfy the following linear system: (32) implying that |ψ E is also an eigenvector of the simpler matrix with eigenvalue −5 p 2 + q 2 . For t > 0, one such vector is . (35) This is identical to the state (18) describing the superfluid in the QCT model. The reason of this equivalence is that, in both cubic and dodecahedral meshes, every site has three neighbors; hence, in the superfluid, the quantities F A , F B , R A , R B are the same. No surprise, then, if also the expressions of ρ and φ turn out to be equal for the two models (the consistency equations are identical). As already commented, the energy of the QDC superfluid is instead a factor 5/2 larger (in absolute terms) than in the QCT model: The main novelty with respect to the QCT model is the existence of a stable supersolid "phase" in the QDC model. To seek for MF solutions having the character of a supersolid, I have first derived the exact equations for the MF parameters ρ A , ρ B , φ A , φ B , without a priori assuming them to be equal in pairs (this is done in Appendix B). Then, I have solved these equations numerically so as to find the region where the supersolid is more stable than the other phases. This task is made simpler by noting that, thanks to a symmetry property of the equations, from a supersolid solution with µ > 3V/2 it is possible to obtain another solution with µ < 3V/2. Indeed, it easily follows from Equation Moreover, for µ = 3V/2 the two solutions share the same grand potential. It turns out that the (t, µ) region where the supersolid is stable is symmetric about the µ = 3V/2 axis, lying across the boundary between the solid "phases" and the superfluid (see Figure 5). To all evidence, the boundary between the supersolid and the superfluid lies at t = 1/3. In fact, the supersolid consists of two distinct "phases", SS1 below µ = 3V/2 (where ρ B < 1/2 < ρ A ) and SS2 above µ = 3V/2 (where ρ A < 1/2 < ρ B ), with φ A < φ B < 1/2 in both. In SS1, cubic sites are more occupied, on average, than co-cubic sites, whereas the opposite occurs in SS2. The two supersolids coexist for µ = 3V/2, for all t in the range 0.25 to 1/3. Hard-core bosons on the triangular lattice have a similar phase diagram [32], but in that case the supersolid extends down to t = 0 and the solid lobes do not overlap each other.
The full T = 0 phase diagram of the QDC model according to the decoupling approximation is reported in Figure 5. Along the straight lines µ = −3t and µ = 3V + 3t, the µ derivative of the grand potential is continuous. Instead, along the curves which respectively separate the cubic and co-cubic regions from the superfluid region, the total number of particles jumps discontinuously. Only at the lobe vertices (V/3, V) and (V/3, 2V) the two "crystal"-superfluid transitions are continuous. The line µ = 3V/2 separating the two "crystalline" regions ends at the point (3V/10, 3V/2) where the two boundary curves (37) cross each other. However, this triple point is only metastable since, in the whole region bounded by the blue circuit of Figure 5, the stable phase is actually supersolid.

Assessment of MF Theory
Let us again reconsider the QCT model at T = 0. Due to the relatively small dimensionality of its Hilbert space (=256), the model can be solved numerically, determining (among others) the exact ground state |g and its eigenvalue (i.e., the grand potential) as a function of t and µ. This is obtained by representing the grand Hamiltonian H on the Fock basis {|x 1 , x 2 , . . . , x 8 } (with x i = 0 or 1 for all i) and diagonalizing the ensuing matrix. An extensive mapping of a few quantum averages enables us to clarify the nature of the "phases" present.
As usual, the sites of the mesh are classified according to what tetrahedral sub-mesh, A or B, they belong to. Then, I compute the average occupancies of A and B sites (corresponding to the MF parameters ρ A and ρ B , respectively), the average values of a A and a B (corresponding to the MF parameters φ A and φ B , respectively), and the superfluid density ρ SF (see, e.g., Ref. [13,22]). In the present case, the latter quantity reads: where a 0 = (1/ √ 8) ∑ 8 i=1 a i is the zero-momentum field operator. Notice that, in a large lattice of M sites, a † 0 a 0 = N 0 is the average number of condensate particles, hence ρ SF = N 0 /M is indeed the condensate density.
As far as the occupancies are concerned, exact diagonalization shows that they are always equal for an A and a B site, with the only exception of t = 0 and 0 < µ < 3V where the occupancies are as in MF theory (i.e., either n A = 1 and n B = 0 or vice versa). In particular, for t > 0 the equivalence n A = n B holds in the whole (t, µ) region pertaining, according to MF theory, to the tetrahedral "phase". Hence, no spontaneous symmetry breaking does really occur in the QCT model, except for t = 0. Indeed, for t > 0 the Fourier coefficients of |g relative to any pair of Fock states equal by A-B inversion are the same. I show in Figure 6 the average site occupancy in the whole space of parameters. We see that the (t, µ) plane is divided in zones where the site occupancy, which overall grows with µ, takes the same constant value. The possible values are k/8, for k = 0, 1, . . . , 8. This outcome is not entirely unexpected considering that the Hamiltonian commutes with the N operator, implying that the ground state (which is non-degenerate for t > 0) should also be an eigenstate of N. Like in MF theory, the occupancy is zero below µ = −3t and 1 above µ = 3V + 3t; in a whole region around µ = 3V/2 the occupancy is 0.5 for all t, with no abrupt transition from "crystal" to superfluid values. Another difference with MF theory concerns the ground-state averages of a A and a B : these are identically zero for all t and µ values, which may seem in stark contrast to the behavior of the MF parameter φ. In fact, it is not these averages that should be monitored in an exact treatment but rather the ρ SF quantity defined at Equation (38). In the top panel of Figure 7, I have reported the ρ SF values computed along a few iso-t lines; each jump discontinuity of ρ SF occurs in coincidence with one of the site occupancy. For comparison, in the bottom panel of Figure 7, I show the MF values of φ 2 along the same lines. We see a clear correlation between the two behaviors, which demonstrates the existence of signatures of superfluidity in a finite quantum system. However, intriguing differences also exist: (i) first observe that ρ SF = 0.125 in the filled mesh, which is an oddity for a perfectly insulating state. In fact, 0.125 is 1/M for M = 8, suggesting that this is an artifact of the finite system size. (ii) Another finite-size effect is the non-zero values of ρ SF in the purported tetrahedral region. However, ρ SF is here sufficiently small so that its peaks in the µ gaps between the tetrahedral and insulator regions remain well visible-in MF theory these gaps fall in the superfluid region. Finally, in Figure 8, I plot the QCT grand potential as a function of µ for a number of t values in the range 0 to 0.5. In the MF curves, cusp singularities are associated with the crossing of first-order transition lines. We see that MF theory systematically overestimates exact values, the more so the larger t is.

Conclusions
I have worked out the zero-temperature phase diagram of two systems of hard-core bosons defined on the nodes of a regular spherical mesh. The interaction is of Bose-Hubbard type, with a further repulsion between neighboring particles. Choosing a suitable mesh, bosons are pushed to form a Platonic "crystal" in a range of chemical potentials. In the QCT model, the mesh is cubic and a tetrahedral "crystal" is formed; in the QDC model, the mesh is dodecahedral and a cubic "crystal" is formed instead.
Using a mean-field approximation, I have obtained fully analytic results for the thermodynamic properties of the two models. In addition to a number of insulating phases, both systems also exhibit a superfluid ground state. In the QDC model, triple coexistence between two solids and the superfluid is superseded by the occurrence of a more stable supersolid phase. Clearly, while the predictions of mean-field theory are generally accurate for an infinite Bose-Hubbard system, deviations will unavoidably be observed in a small system, where no true singularity can occur. However, discrepancies are probably less strong for bosons on a regular spherical mesh, which has equivalent sites and is devoid of natural boundaries.
To check this expectation for the QCT model, I have diagonalized its Hamiltonian exactly, finding the ground state and a number of ground-state averages. Overall, the exact T = 0 behavior of the finite system does not depart much from theory, though obviously phase-transition lines are now reduced to simple crossovers; however, clear marks of superfluidity have been detected in the same region of parameters where MF theory predicts it to occur. In conclusion, I hope that the present work may stimulate research aimed at the realization of a new experimental platform for ultracold bosons that is akin to a regular spherical mesh.
Funding: This research received no external funding.

Conflicts of Interest:
The author declares no conflict of interest.

Appendix A. On the Solutions to an Algebraic Equation
In this Appendix, I show how to determine the eigenvalues of a matrix like (29), namely the four real solutions λ to the characteristic Equation (30) with a + d = b + c ≡ s.
Upon doing the algebra, the characteristic equation turns out to be: Let us now reshuffle the lhs of (A1) term by term: The MF Hamiltonian is Equation (29) with real φ A and φ B . Therefore, the eigenvalue equation is as in (30), with a = E 0 = 48tφ A φ B + 12tφ 2 B − 24Vρ A ρ B − 6Vρ 2 B , b = E 0 + 12(2Vρ A + Vρ B − µ), c = E 0 + 8(3Vρ B − µ), and d = E 0 + 4(6Vρ A + 9Vρ B − 5µ); moreover, u = −12t(2φ A + φ B ) and v = −24tφ B . The minimum eigenvalue is like in the first line of Equation (31), but its explicit form in terms of MF parameters is now Probably, the simpler method to solve these four coupled non-linear equations in four unknowns is to minimize a suitable non-negative function constructed in such a way as to vanish when the (A15) are all fulfilled. This is easy to do numerically with a computer. By this method, I have drawn the supersolid boundaries in Figure 5.