Real Space Triplets in Quantum Condensed Matter: Numerical Experiments Using Path Integrals, Closures, and Hard Spheres

Path integral Monte Carlo and closure computations are utilized to study real space triplet correlations in the quantum hard-sphere system. The conditions cover from the normal fluid phase to the solid phases face-centered cubic (FCC) and cI16 (de Broglie wavelengths 0.2≤λB*<2, densities 0.1≤ρN*≤0.925). The focus is on the equilateral and isosceles features of the path-integral centroid and instantaneous structures. Complementary calculations of the associated pair structures are also carried out to strengthen structural identifications and facilitate closure evaluations. The three closures employed are Kirkwood superposition, Jackson–Feenberg convolution, and their average (AV3). A large quantity of new data are reported, and conclusions are drawn regarding (i) the remarkable performance of AV3 for the centroid and instantaneous correlations, (ii) the correspondences between the fluid and FCC salient features on the coexistence line, and (iii) the most conspicuous differences between FCC and cI16 at the pair and the triplet levels at moderately high densities (ρN*=0.9, 0.925). This research is expected to provide low-temperature insights useful for the future related studies of properties of real systems (e.g., helium, alkali metals, and general colloidal systems).


Introduction
The study of equilibrium triplet structures in 3D N-particle systems with quantum behavior remains a pending task in condensed matter research at low temperatures. Apart from a number of early developments focused mainly on the proposal and indirect testing of the so-called closures [1][2][3][4][5][6][7][8][9][10][11] or on the use of alternative order parameters [12], just a few computational works based on Feynman's path integrals (PI) [13] can be found in the recent literature on this field [14][15][16][17][18]. Not much is known about the behavior of quantum triplets, hence the interest in undertaking this task. This is not only a logical step further in current statistical mechanics, allowing one to formulate thermodynamic properties beyond the pairwise approach [19], but also it is central to outstanding condensed matter properties. Among the latter, one can mention the following: phonon-phonon interactions in helium-II [4], the N-particle interpretation of fluid entropies [20][21][22], multiple scattering [23], theories of phase transitions [24,25], and glassy dynamics [26][27][28]. Although the whole PI quantum triplet task is computationally daunting at the present time, one can always seek to identify the main triplet features that may serve as a guide for the necessary future work on this topic.
The triplet topic encompasses both real-space {r}-triplets and reciprocal (Fourier)-space {k}-triplets. Nevertheless, there is no direct experimental determination of triplet functions [29,30]. Thus, one must resort to theoretical computations for extracting this sort of information, which makes these computations the only "experimental" method of solving the triplet problem. In the quantum In stark contrast, the classical system only shows one category of a physically significant structure [56]. Therefore, it is easy to understand the much higher complexity and computational cost of the quantum problem when treated in depth.
The aim of the present article is to expand the study of the PI triplet features in many-body quantum systems [14][15][16][17][18]. The system selected is that composed of quantum hard spheres (QHS system hereafter). Hard spheres are known to be a very useful reference. They have been used to model classical and quantum systems, ranging from helium [33,36,66,67,[71][72][73] to complex colloids [74,75], and under very different conditions (i.e., fluid, boson superfluid, superlattices, solids, and suspensions). This work concentrates on the real space instantaneous and centroid triplets, leaving aside the total thermalized-continuous linear response case to keep the related computations affordable [17]. As stressed earlier, in understanding triplets through closures, a thorough consideration of the pair structures is needed. This benefits the triplet-closure computations and the analysis of the correspondences between the salient pair and triplet features (in different phases or within the same phase). Therefore, the necessary attention is also paid to the pair prerequisite.
The exact computational method chosen is PIMC, which avoids the PIMD difficulties linked to the discontinuity of the hard-sphere potential, and it is complemented with the closures: KS3, JF3, and their average AV3 = (KS3 + JF3)/2. A wide range of QHS fluid and solid conditions, within the purely diffraction regime, is studied: reduced de Broglie wavelengths 0.2 ≤ λ * B < 2 and reduced number densities 0.1 ≤ ρ * N ≤ 0.925, where λ * B = h/(2πMk B Tσ 2 ) 1/2 , ρ * N = Nσ 3 /V, σ and M being the hard-sphere diameter and mass, respectively. The specific {r}-space targets pursued are the following: (i) Analyzing in detail the potential usefulness of AV3 for quantum fluid triplet studies. This is prompted by the encouraging results obtained in [17] for liquid para-hydrogen and liquid neon. (ii) Gaining triplet structural insights from the comparison, in the short-medium range of distances, between the coexisting fluid and FCC (face-centered cubic) solid [66,67]. (iii) Comparing the salient triplet features of the cubic solids FCC and cI16 at moderately high densities, (ρ * N = 0.925, λ * B = 0.2) and (ρ * N = 0.9, λ * B = 0.8). The so-called cI16 lattice is a distorted superstructure of the body-centered cubic (BCC) lattice [76,77]. (Hard-sphere BCC lattices are known to be mechanically unstable in both classical and quantum applications [78,79]). One also notes that cI16 phases have been reported for Li and Na at very high pressures [76,77], hence the interest of this point. (iv) In connection with (iii), there is the question of establishing a clear identification of the QHS bcc-qIII phases, observed in [67,78], as genuine cI16 phases. (The bcc-qIII phases arise from the PIMC evolution of initially perfect BCC lattices). The cI16 lattice has been identified recently for classical hard spheres in the insightful simulation work reported in [80], and the patterns of the related results suggest that bcc-qIII is in fact cI16. Proof of this is given in this article, which adds more value to the QHS system for modeling quantum solid-solid changes of phase [67] and enhances the meaning of the related triplet calculations.
It is hoped that the reported "experimental" results will form a useful basis for comparison when extensive studies of triplets in real quantum systems are undertaken. The outline of this article is as follows. Section 2 contains the theoretical background. Section 3 describes the relevant computational details. Section 4 gives the results and their discussion, and Section 5 collates the main conclusions of this work.

Path Integral Monte Carlo (PIMC)
The canonical partition function Z(N, V, T) of a quantum monatomic system (number density ρ N = N/V), under conditions where exchange interactions can be neglected, can be approximated by the accurate PI formula (Tr = trace) [31][32][33] where H N is the Hamiltonian, assumed normally to be composed of one-and two-particle terms, M is the particle mass, β = 1/k B T, P is the discretization in beads of the necklace representing and actual particle j, r t j denotes the real space coordinates of bead t (t = 1, 2, . . . , P) belonging to necklace j, and W NP is the "effective potential" ruling the whole set of N × P beads (hereafter all of them equivalent: X = P). In what follows, the optimal P will be assumed. In addition, it is worthwhile to note that (a) consecutive beads, t and t + 1, in a necklace are separated by β /P in Euclidean imaginary time β ; (b) then, a given bead labeled t is associated with the imaginary time β t/P; and (c) the cyclic property t + 1 = P + 1 ≡ 1 is satisfied.
In the case of the QHS system, an appropriate choice for W NP is based on Cao-Berne's CBHSP pair action [35], and it can be written as [14,78] W CBHSP In Equation (3a), one recognizes the superposition of the free-particle behaviors [13]. Equation (3b) shows Cao-Berne's correction, where the adjacent-bead effects are to be noted. Equation (3c) is the expression of the singular hard-sphere potential extended over all the pairs of necklaces, which interact in an equal "t-time" bead-to-bead fashion (ET). The symbols P * in the sum and product above mark the t−cyclic property already mentioned. For the specific thermodynamic property formulas that can be derived from CBHSP, the reader is referred to [67,78]. At this point, it is important to give the definition of the CBHSP centroid of a given necklace j This quantity plays an important role in PI theoretical developments [13,47], in particular in (a) the appealing centroid approaches to quantum dynamics [81][82][83][84][85]; (b) the exact formulation of the equation of state of quantum fluids [39,86]; and (c) the characterization of quantum solid phases [67].
A key feature of the QHS system is that its state points can be uniquely characterized by giving two parameters, namely the reduced number density ρ * N and the reduced de Broglie thermal wavelength λ * B . This fact was early realized at the level of semiclassical treatments based on -expansions (see, for instance, [87][88][89][90]). Within PI, the same fact is just a property contained in the mathematical structure of the QHS partition function, regardless of the particular proper form that W NP may take (see [47] for a discussion of QHS propagators). Accordingly, the QHS system properties can be expressed in reduced units, thereby being independent of the actual parameters M, σ, T, and ρ N employed [66,67]. Therefore, for example, internal energies E can be given as E * = E/RT, and pressures p can be given as p * = pMσ 5 / 2 . For the pressure, an indication to guide the interested reader will suffice: when using different sets of parameters to define the state points 1 and 2, if ρ * N , λ * , then necessarily (PV/RT) 1 = (PV/RT) 2 , and also p * 1 = p * 2 . The same general type of argument applies to the real space structures g n q 1 , q 2 , . . . , q n , for which, when reporting QHS system results, it is most useful to do it using interparticle distances in reduced form: r * 12 = q 1 − q 2 /σ. In order not to burden the notation, the formulation of the structural concepts below will utilize the distances and related quantities in their non-reduced forms, as in Equations (1) to (4).
Another technical point seems worth placing here. It is related to the three-(and higher-order) particle contributions to the quantum Hamiltonian H N of the system, which may yield more complete and elaborate forms for the propagators and W NP . While this is a question that can be tackled in various ways when continuous interparticle potentials are involved [91][92][93], to this author's knowledge, no QHS system PI actions beyond the pair level are available, and such an extension remains intractable for now. Nevertheless, because of the strong similarity between helium atoms and quantum hard spheres, the related effects on the QHS system can be expected to become significant at very high solid-phase densities (and sufficiently low temperatures) [33]. Furthermore, owing to the QHS infinite repulsion at the hard core, Equation (3c), the wave functions of the QHS system must vanish for interparticle distances r ≤ σ (i.e., there can be no tunneling); hence, quantum hard spheres repel one another before "classical contact" can occur [89,90]. (Within PI, this means that the related forbidden region brought about by Equation (3c) arises only for the "equal-time" bead g n correlations). Therefore, given the lack of any attraction, the "preemptive" QHS repulsions can be expected to cause a strong impediment to the coming together of triplets of quantum hard spheres. Using the quantum diffraction parameter γ = ρ * N λ * 3 B , the latter triplet effects should not play any significant role unless γ becomes truly high. The largest value of γ in this work is 2.8, which is compatible with the QHS pair modeling of normal fluid and solid helium-4 [66]. Therefore, the pair-level CBHSP approach can be deemed adequate to compute structures under the fluid and solid conditions investigated in this work.

PI Triplet Structures
Within PI-CBHSP, the centroid (CM3) and the instantaneous (ET3) three-point number densities can be cast as the ensemble averages [17] ρ (3) where one notices that (i) the multi-index summations run over the whole set of permutations of N particles taken three at a time; (ii) the instantaneous case contains a further P average involving "equal-time" beads in different necklaces; and (iii) these definitions are completely general, since they depend on the position vectors of the representative set of three particles and can be applied to the statistical description of monatomic systems, which are either fluid or solid. Due to the high computational cost, no attempt is made in this work at calculating total thermalized-continuous linear response triplets [14,17]. For homogeneous and isotropic fluids, one finds simpler formulas [17] ρ (3) CM3 q 1 , q 2 , q 3 = ρ 3 N g CM3 q 1 − q 3 , q 2 − q 3 = ρ 3 N g CM3 (r 12 , s 13 , u 23 ), where the triplet correlation functions g CM3 and g ET3 depend only on the three generic interparticle distances: r 12 = q 1 − q 2 , s 13 = q 1 − q 3 , and u 23 = q 2 − q 3 . This exact reduction from nine to three independent variables makes the intricate triplet problem more accessible for the study of monatomic fluid state points [14][15][16][17]. Rigorously speaking, the related exact framework for a monatomic solid is contained in Equations (5) and (6). Nevertheless, affordable approximations to this even more costly problem can be obtained by applying Equations (7) and (8). Actually, such an approach is consistent with the same idea, already exploited successfully, at the pair level in the study of regular solid phases, since the g CM2 (r) and g ET2 (r) pair radial functions retain many significant traits of the underlying solid structure [66,67,80]. Furthermore, as a first step, the use of Equations (7) and (8) facilitates the comparison of the global salient triplet features appearing in different solid phases.

Additional Pair Structural Quantities
To supplement the PIMC triplet calculations in the canonical ensemble, the following quantities can also be computed: (a) The pair radial functions for the centroid (CM2) and the instantaneous (ET2) correlations, in both the fluid and the solid phases [47]. Their PI ensemble averages can be cast as where r 12 = q 1 − q 2 .
(b) The pair static structure factors S  ET (k) associated with the foregoing pair radial structures in the fluid phase [47] S (2) where h 2 = g 2 − 1 stands for the corresponding pair total correlation function, and c (2) (k) stands for the corresponding pair direct correlation function in Fourier space (k = |k|). These structure factors can be fixed with great accuracy, at a very low cost and for every k ≥ 0 wave number [48], via the Ornstein-Zernike framework [94][95][96] developed by this author [47,86,97,98]. Apart from their intrinsic usefulness, they are decisive in achieving a number of significant improvements in the study of fluids with quantum behavior [39,[47][48][49]. In particular, S CM (k) and S ET (k) can be utilized for (i) extending the ranges of the simulated g CM2 (r 12 ) and g ET2 (r 12 ) [47], which serves to perform triplet closure computations; and (ii) gaining insight into their associated triplet structure factors S  ET (k 1 , k 2 ) [17,18]. (c) In simulation work using cubic boxes, the PI sample size is composed of N S necklaces, each with P beads, enclosed in a volume V S = L 3 . To characterize solid phases, one can employ the normalized-to-unity solid-phase configurational structure factors at the centroid and instantaneous pair levels [67,99,100]. They can be written as and are taken at their maximal values arising from the simulation runs [67,78]. In these simulation conditions, the wave vectors k to be analyzed must be commensurate with the box, which means k = 2πL −1 k x , k y , k z , where the components k x , k y , k z take integer values [56]. In connection with this, notice that cubic-based perfect lattices can be associated with sets of three commensurate wave vectors, {k w } n = {k 1 , k 2 , k 3 } n , which are maximal in that: (i) For the perfect FCC and BCC lattices, one can single out sets {k w } n such that they reach the maximum value, S 2 (k max ) = 1, (w = 1, 2, 3). For a perfect cI16 lattice, which is not so highly regular, one obtains S 2 (k max ) < 1, (w = 1, 2, 3), as will be shown later on. (ii) The following result holds Therefore, comparison of the above standard perfect-lattice results with those arising from the simulated cubic solid phase allows one to identify its type and relative order. Obviously, the values of the simulated configurational structure factors are lower than the perfect reference values; they appear associated with each of the three maximal wave vectors and are close to one another, but, as a rule, they are not equal: one of them can be singled out as the maximum, whereas the other two remain slightly below [67,78]. As a guide for quantum work [99], the following centroid estimates are worth quoting: CM2 (k max ) for partially crystalline solids, while typically S CM2 (k max ) should be between the two foregoing limits). It is important to stress that although somewhat expensive to calculate, the quantities S ET2 (k max ) are global for the simulation sample. Therefore, in this context, these quantities seem more complete than local-order parameters (e.g., the rotationally invariant Q l ) [67,80,101].
Before going any further, it is convenient to consider the general issue of the simulation sample size N S for the solid phases, thus allowing one to introduce cI16 basic details. The conditions for FCC and BCC are well-known, and for N S > 100, one finds: (i) N S (FCC) = 4n 3 , with n = 3, 4, 5, . . .; and (ii) N S (BCC) = 2n 3 , with n = 4, 5, 6, . . .. However, the case of cI16 is not so standard. cI16 is a distortion of BCC and is characterized by the so-called fractional displacement parameter, which is usually denoted by x [76,77], so as to have the particles occupying the 16c Wyckoff site (x, x, x) of the space group I43d. This means that its body-centered unit cell does contain 16 particles. Consequently, there are some extra restrictions that may make the N S (cI16) values different from those of BCC. Thus, again for N S > 100, one finds (iii) N S (cI16) = 16n 3 , with n = 2, 3, 4, . . .. The reader is referred to [76,77,80] for specific details.

Closures for Fluid Triplets
The two basic closures analyzed in this work are Kirkwood superposition KS3 and Jackson-Feenberg convolution JF3. Both can be applied to the fluid centroid (CM3) and instantaneous (ET3) triplet correlations. Their expressions can be written as follows [1,4]: (18) where v j4 = q j − q 4 , h 2 = g 2 − 1, and g 2 = g CM2 or g ET2 . Although explicitly stated in Equation (18), it is important to remark that JF3 lacks the triplet-product term h 2 (r 12 )h 2 (s 13 )h 2 (u 23 ), which should appear in an h 2 −expansion. This absence has deep consequences as will be shown in this article. An easy and direct way to recover such contribution (half of it), while at the same time keeping the convolution integral (half of it) contained in Equation (18), is via the average closure AV3 that reads as As regards the properties of these closures, suffice it to say that (i) KS3, JF3, and AV3, satisfy Equations (9a), (9c), and (9d); and (ii) only KS3, as induced by its construction, satisfies Equation (9b), which is a special case of the general triplet behavior g 3 → 0 when two particles approach closely each other [14][15][16][17].

Computational Details
The main target of this work is the determination of QHS equilateral and isosceles triplet correlations (centroid and instantaneous), namely the types of functions g 3 (r, r, r) (or g EQ 3 for brevity when necessary) and g 3 (r, s, s). For the sake of interpretation, they are complemented with the additional structural properties discussed in Section 2.3. The state points studied are shown in Table 1. They span a wide range of conditions, from the normal fluid phase to the distinct solid phases FCC and cI16. Special attention is paid to the two sides of the fluid-FCC coexistence line, as determined in [67] λ * B ≤ 0.8 and [66] λ * B > 0.8 . Moreover, the study is extended to (i) fluid state points under very strong diffraction effects λ * B ≈ 2 , with a view to establishing a meaningful correlation of triplet structures when going toward the change of phase by increasing ρ * N at constant temperature, and (ii) the lattices FCC and cI16 at (ρ * which are conditions that are significantly far from the very high-density regions.
The necklace normal mode algorithm [62,63] is used to generate the collective P movements of a given necklace, with a Metropolis acceptance ratio of 50%. (As in previous works, the actual hard-sphere parameters are M = 28.0134 amu and σ = 3.5 Å; 1 Å = 10 −10 m). The necklace sample sizes N S are 864 for the fluid and the FCC solid phases, and 1024 for the cI16 solid phases. The quantum P convergence for the results is studied as shown below (12 ≤ P ≤ 40). One kpass is defined as a set of 10 3 N S × P attempted bead moves, and one Mpass is then 10 3 kpasses. After equilibration, most of the simulation runs are arranged into 40 blocks for the g 2 calculations and 30 blocks for the g 3 calculations.
The respective block sizes are (i) 92.6 kpasses for the fluid simulations; (ii) 92.6 kpasses for the FCC simulations; and (iii) 78.125 kpasses for the cI16 simulations. Therefore, the run lengths associated with the g 2 and g 3 calculations are in between 2.34 Mpasses and 3.7 Mpasses. (The extra simulations using P = 36 and 40 have lengths of about 1 Mpass). Block subaverages for g 2 and g 3 are obtained by gathering statistics every 5000 (g 2 )/7000-8000 (g 3 ) configurations generated. The configurational structure factors given by Equations (14) and (15) are analyzed four times per block, at equally spaced intervals, by recording the ten largest values for the final assessment. To do so, triplets of integers k x , k y , k z are monitored in the mesh 25 ≤ k 2 x + k 2 y + k 2 z ≤ 200, with the components in −10 ≤ k ν ≤ 10 (symmetry properties allow one to reduce the calculations). Given that the information provided by the correlation functions, complemented with that arising from the structure factors, is more than sufficient to characterize the current solid structural results, the Q l order parameters [101] are not evaluated, thereby alleviating the considerable computational effort involved in this work.
The pair and triplet sructures g 2 and g 3 are fixed in the established ways using histograms. The case of g 2 is straightforward and well-known [56], and the simulations are utilized as the reduced width of the bins ∆ * = 1/35 (or σ/35 =0.1 Å). However, the case of g 3 includes a good number of subtleties [57,58]. The detailed description of the related procedure can be found in [14]. For the current purposes, suffice it to say that for a triplet of distances (R, S, U), the basic g 3 expression is given by ; ET3 and CM3 (20) where (∆n T ) is the number of times mutual distance triplets lie within the ranges R − ∆ < r 12 ≤ R, S − ∆ < s 13 ≤ S, and U − ∆ < u 23 ≤ U, and (∆V) RSU stands for the appropriate volume element [58]. Once again, in these calculations, ∆ * = 1/35. The histogramming of triplets extends up to distances r 12 , s 13 , and u 23 , which are < L/4. Statistical errors (one-standard deviation) for the average structures computed with PIMC are fixed with the corresponding block subaverages. For example, for the first peaks heights of g 2 and g 3 , the errors remain below 1% for most of the present calculations. In this connection, Table 2 gives some representative g 3 results (mean first peaks (FP)) in the close vicinities of the absolute maxima of the structure indicated, together with the associated errors. (More on this in the Supplementary Materials). Note that the P convergence is influenced by both λ * B and ρ * N . For the fluid and FCC state points on the coexistence line, under the most extreme quantum conditions studied herein (λ * B = 1.9832, γ 2.7 − 2.8), P = 30 is sufficient to produce practical convergences in the centroid and in the instantaneous functions. For the solid state points at densities ρ * N = 0.9, 0.925, it is worthwhile to note that there is a slowing down of this convergence with decreasing temperatures (λ * B = 0.2 → 0.8) , which becomes more noticeable (a) for the triplet centroid quantities and (b) for the cI16 lattice, its openness playing a significant role in the fixing of the final particle distributions.

Closure Calculations
The current calculations at the actual fluid state points on the coexistence line use the new PIMC information obtained with sample sizes larger than those employed in [49]. (The new and the former results are in excellent agreement). The JF3 convolution integrals involve the h 2 extension to longer ET (k). The convolutions can be obtained by employing a well-known expansion in Legendre polynomials P n [10,24] where φ is the angle between s 13 and u 23

Results and Discussion
The results reported in this section are complemented with data in the Supplementary Materials.  (Figure 1c,d) display the expected traits of FCC lattices. General comments on these pair radial functions are (i) the higher order in the solid functions that does not disappear with increasing distances; (ii) the outward shift and smoothing of the features with increasing quantum effects (on the coexistence line analyzed, one has 0.006 < γ < 2.81); and (iii) the proximity between the locations of the fluid and solid first maxima (also between the dominant second-maximum regions), revealing that the system is ready to effect the change of phase. It is also interesting to note in passing that on the fluid side, the absolute maxima of the pair structures show dependences upon γ that can be fitted in the form g 2 (Max) = aγ −b , the associated linear correlation coefficients r corr. being reasonably good:   (14) and (15): the fluid phase maximal values obtained remain ( ) < 0.1. Moreover, Table 3 contains the observed variations in the maximal values of ( ) corresponding to the FCC centroid and instantaneous cases. For the current calculations, a representative FCC-set of maximal wave vectors can be defined by their k-integer components: (−6,6,6), (−6,6, −6, ), (6,6,6) .  Table 3. Solid phase variations in the maximal values of the centroid (CM2) Equation (14) and instantaneous (ET2) Equation (15) configurational structure factors at the pair level fixed with PIMC.  Figure 2 shows the pair radial correlation functions, centroid and instantaneous, of the FCC and cI16 state points at the moderately high densities ρ * N = 0.9, 0.925. There is a sharp contrast between the FCC and cI16 structures, since the usual coordination shells existing in the highly regular FCC lattice are absent from cI16. The most characteristic trait of cI16 is, perhaps, the presence of a convoluted inner structure, with two conspicuous big dips, for distances below r * ≈ 2.5. The FCC solid structures (Figure 2a,c) are the "compression" (at constant temperature) of the corresponding FCC structures on the coexistence line. The current cI16 results (Figure 2b,c) agree feature for feature with the pair structures displayed by the bcc-qIII phases in [67]. (Differences between the first peaks are due to the B-spline smoothing carried out in Figure 9 of [67]; see the Suppplementary Material for non-smoothed data). This deserves to be highlighted, since the PIMC-QHS origins of both types of structures are not related: the former bcc-qIII phases arose from the evolution of initially perfect BCC lattices (N S = 432), while the present (delocalized) cI16 phases are just the results obtained from the evolution of initially perfect cI16 lattices (N S = 1024). To complete the foregoing information, Table 3  There is still the further question related to the characterization of cI16 phases via the fractional displacement parameter . In the quantum case, the delocalization makes this task a three-fold one, since there are three types of distinct structures. Given the current scope, only the centroid and instantaneous estimates are determined in this work. A convenient way is through the tabulation for perfect cI16 lattices of ( , ( ) ( )), which can be computationally fixed by varying . Thus, for the interval 0.025 ≤ ≤ 0.038, using Δ = 0.001, the related parabolic least-squares fitting (better than just linear) leads to ( ) ( ) = 1.0931 − 8.269 − 108.75 ; CM2 and ET2,

FCC PHASE on the Coexistence Line
which guarantees absolute errors of orders ≤ 10 in the estimated values of the reference maximal There is still the further question related to the characterization of cI16 phases via the fractional displacement parameter x. In the quantum case, the delocalization makes this task a three-fold one, since there are three types of distinct structures. Given the current scope, only the centroid and instantaneous x estimates are determined in this work. A convenient way is through the tabulation for perfect cI16 lattices of x, S interval 0.025 ≤ x ≤ 0.038, using ∆x = 0.001, the related parabolic least-squares fitting (better than just linear) leads to S (C) 2 (k max ) = 1.0931 − 8.269x − 108.75x 2 ; CM2 and ET2, (22) which guarantees absolute errors of orders ≤ 10 −4 in the estimated values of the reference maximal structure factors (for higher precision, the reader is referred to the tabulation in the Supplementary Materials). Note that the higher the x is, the lower the S (C) 2 (k max ) becomes, as expected. In this regard, note that for x = 0, which is out of the above interval, one must retrieve the perfect BCC limit S  [76,77,80]; and (c) the CM2-ET2 differences increase with the quantum effects. Another point to consider here is related to the fact that samples of classical hard spheres can be "squeezed" more than samples of quantum hard spheres, because of the latter's "preemptive" repulsions. This means that via low temperatures, one can expect QHS-cI16 phases to appear for lower densities than in the classical hard-sphere system ρ * N ≥ 1.1 [80], which is indeed the case.  Figure 3, which collects results at two state points along the lowest isotherm λ * B = 1.9832. First, as occurred on the pair level, the centroid CM3 features are far more pronounced than those of the instantaneous ET3 case. Second, and associated with the equilateral data, one notes that the first maximum and the first minimum positions of a given g 3 (r * , r * , r * ) occur in the close vicinities of the corresponding first maximum and first minimum of the associated g 2 (r * ) shown in Figure 1. Third, although the closures KS3 and JF3 fail to reproduce the exact PIMC behavior, their average AV3 shows a remarkable performance for both the centroid and the instantaneous correlations. Fourth, as the density increases along isotherms, and when going toward longer distances, AV3 loses predictive power to fit the profiles of the isosceles correlations g 3 (r * , s * , s * ). In relation to this, see Figure 3d, where s * = s * M is such that g 3 s * M , s * M , s * M absolute equilateral maximum. Finer equilateral facts follow.

Triplets in the Fluid Phase
(i) Figure 3a,b displays explicitly, at state point (ρ * N = 0.1, λ * B = 1.9832), the equilateral asymptotic behavior g 3 (r * , r * , r * ) → 1 with increasing r * for the PIMC centroid and instantaneous correlations. (ii) Figure 3c illustrates the isosceles asymptotic behavior g 3 (r * , s * , s * ) → g 2 (r * ) , when the two s * distances increase. (iii) As seen, the short-range behavior of AV3 is non-correct (due to that of JF3), whereas KS3 behaves properly. (iv) At constant temperature, there is a sharpening and shifting inwards of the structures with increasing density. For example, at λ * B = 1.9832 in the vicinities r * FP of the equilateral first maxima, the g EQ 3 = g 3 (r * , r * , r * ) behave as follows: (1) ρ * N = 0.1, r * FP = 1.9, g EQ CM3 = 2.16 and r * FP = 2, g EQ ET3 = 1.41 ; (2) ρ * N = 0.3, r * FP = 1.5571, g EQ CM3 = 12.65 and r * FP = 1.5429, g EQ ET3 = 3.54 ; and (3) ρ * N = 0.348, r * FP = 1.5, g EQ CM3 = 19.74 and r * FP = 1.4714, g EQ ET3 = 4.51 . An analogous behavior can be observed at the pair level. (Use σ = 3.5 Å and rounding-off to two decimal places to transform the foregoing r * into the actual r of the (M, σ) system utilized in the current calculations, e.g., r * = 1, 5571 → r * = 5.5 Å). minimum of the associated ( ) shown in Figure 1. Third, although the closures KS3 and JF3 fail to reproduce the exact PIMC behavior, their average AV3 shows a remarkable performance for both the centroid and the instantaneous correlations. Fourth, as the density increases along isotherms, and when going toward longer distances, AV3 loses predictive power to fit the profiles of the isosceles correlations ( * , * , * ). In relation to this, see Figure 3d, where * = * is such that ( * , * , * ) ≅ absolute equilateral maximum.  Figure 3c illustrates the isosceles asymptotic behavior ( * , * , * ) → ( * ), when the two * distances increase. (iii) As seen, the short-range behavior of (Use = 3.5Å and rounding-off to two decimal places to transform the foregoing * into the actual of the ( , ) system utilized in the current calculations, e.g., * = 1,5571 → * = 5.5Å). Figure 4 shows the equilateral correlations at three representative state points on the fluid side of the coexistence line. The aforementioned trends of KS3 and AV3, as compared to PIMC, appear again for both types of correlations CM3 (Figure 4a) and ET3 (Figure 4b). In going from higher to lower densities/temperatures on the fluid side, one observes that the larger the quantum effects, the flatter the structural triplet features become.       (Figure 4b). In going from higher to lower densities/temperatures on the fluid side, one observes that the larger the quantum effects, the flatter the structural triplet features become. Table 4 contains the absolute maxima, fixed with quadratic interpolations of the adequate PIMC data, of the fluid equilateral correlations. (See the Supplementary Materials for more related numerical data). Once more, in an attempt to connect the foregoing maximum behaviors with the quantum parameter γ, one notes that simple empirical decay fittings g EQ 3 (Max) = aγ −b can be found for the centroid and for the instantaneous cases, their associated linear correlation coefficients r corr. being reasonably good. Thus, one finds for the centroid case a 23.6702, b 0.1877, r corr. = −0.9959 and for the instantaneous case a 6.437, b 0.3449, r corr. = −0.9999. This general pattern is to be regarded as a reflection of the very same observed at the pair level. Three additional points are worthwhile to mention: (i) the quality of this type of fitting remains comparable if one tries the modification g EQ 3 (Max) = aγ −b + c; (ii) exponential decays, e.g., g EQ 3 (Max) = aexp(−bγ), give poor fittings; and (iii) the potential energy discontinuity at r * = 1 precludes one from retrieving the classical limit at λ * B = 0. Although there is no apparent physical basis for the empirical γ pattern found, this line of thought might be well worth exploring in future work. Table 4. Absolute maxima of the PIMC equilateral correlations g EQ 3 = g 3 (r * , r * , r * ) on the fluid and the FCC sides of the QHS coexistence line. Discretizations at λ * B = 0.8 and 1.9832: P = 12 and 30, respectively. r * = r/σ. Four decimals shown in g EQ 3 to avoid rounding-off errors.  Figure 5 contains a quick description of the isosceles correlations g 3 (r * , s * , s * ) at two representative fluid state points, for the centroids CM3 in panels (a)-(c) and for the instantaneous ET3 in panels (b)-(d). Three especial r * distances are selected from the g 3 (r * , r * , r * ) information obtained at each state point, namely r * FP , r * FV , and r * SP , which are positions in the close vicinities of the equilateral maxima and minima: first maximum (FP), first minimum (FV), and second maximum (SP), respectively. Apart from the expected AV3 unphysical behavior for r * ≤ 1, the good overall performance of AV3 is certainly surprising. Two weak points are to be remarked. First, the AV3 (and KS3) behavior for low s * distances, 1 < s * < 1.5, when r * increases: for example, at r * SP where the closure maxima are overestimated. (This is directly related to the AV3 trend displayed by the upper profile plot in Figure 3d). Second, Figure 6 shows a detailed image of the isosceles deterioration of the PIMC-AV3 agreement with increasing densities, the worse results for centroids being a consequence of this key fact (centroids mimic a fluid at a higher density than the actual one).   Table 4 and Figures 7 and 8 show selected results for the PIMC equilateral and isosceles correlations of FCC state points on the solid side of the fluid-FCC coexistence line, within the  Table 4 and Figures 7 and 8 show selected results for the PIMC equilateral and isosceles correlations of FCC state points on the solid side of the fluid-FCC coexistence line, within the approximations obtainable from Equations (7) and (8). For visualization purposes, the associated PIMC fluid results are also incorporated into these figures. Figure 6. QHS fluid centroid (CM3) and instantaneous (ET3) * profiles of the heights in the close vicinities of the first peaks of the isosceles correlations at two selected state points on the fluid-FCC coexistence line. (a) Fluid functions at ( * = 0.789, * = 0.2). (b) Fluid functions at ( * = 0.533, * = 0.8). * = distance in the close vicinity of the absolute maximum of the equilateral correlations. Acronyms for methods as in Figure 3. The vertical line at * = 1 marks the position of the hard core. Table 4 and Figures 7 and 8 show selected results for the PIMC equilateral and isosceles correlations of FCC state points on the solid side of the fluid-FCC coexistence line, within the approximations obtainable from Equations (7) and (8). For visualization purposes, the associated PIMC fluid results are also incorporated into these figures. A general idea can be obtained by observing Table 4. The fluid and solid absolute maximum positions are close to one another, and the structures are shifted outwards with increasing quantum effects. Moreover, higher values appear on the solid side (e.g., at * = 0.2, ≈ +39% for CM3, and ≈ +51% for ET3). This trend is far more pronounced for the centroid correlations, the ratio  Figure 9 and Table 5 contain equilateral PIMC results for the FCC and cI16 state points in the region of moderately high densities ( * = 0.9,0.925). The centroid CM3 and the instantaneous ET3 A general idea can be obtained by observing Table 4. The fluid and solid absolute maximum positions are close to one another, and the structures are shifted outwards with increasing quantum effects. Moreover, higher g 3 values appear on the solid side (e.g., at λ * B = 0.2, ≈ +39% for CM3, and ≈ +51% for ET3). This trend is far more pronounced for the centroid correlations, the ratio increasing monotonically with λ * B (e.g., at λ * B = 1.9832, ≈ +102% for CM3). However, for the instantaneous ET3 correlations, such a ratio is not monotonic; it goes through a maximum (at λ * B = 0.4, ≈ +62%) and then falls monotonically (at λ * B = 1.9832, ≈ +20%). These behaviors can be ascribed to the two effects present on the coexistence line. On the one hand, there is the decreasing density, which contributes to diminishing the structural features. On the other hand, there is the increasing delocalization with λ * B , which makes PI structures become more and more smeared out, the instantaneous case being always much more sensitive to this. As regards the question of finding a γ−fitting of the solid equilateral absolute maxima, the situation is less clear than on the fluid side (γ is slightly higher on the solid side). Although one can obtain reasonable dependences g EQ 3 (Max) = aγ −b (r corr. < −0.991), on closer inspection, these fittings cannot cope with the apparent inflection in 0.13 < γ < 0.3 (or in 0.6 < λ * B < 0.8) for centroids g EQ CM3 (Max), nor with the large discrepancies for low γ between the original and the estimated instantaneous values g EQ ET3 (Max). In Figure 7, one observes that the equilateral FCC and fluid g 3 (r * , r * , r * ) patterns are qualitatively similar within the first maximum regions. It is also noticeable that the FCC state points develop easily identifiable peak structures with increasing distances (r * 2). The main two maxima of the FCC equilateral triplets can be put into direct correspondence with the main two maxima obtained at the FCC pair level (Figure 1), since they appear located close to one another.

Triplets in the FCC and cI16 Denser Solid Structures
The FCC g 3 (r * , r * , r * ) display deep first valleys, almost at the zero-ground level, appearing for both the centroid and the instantaneous structures, e.g., for centroids and ρ * N = 0.573, the region in Figure 7 located in 1.6 r * 2.1. In general, this feature is far more pronounced in the centroid structures than in the instantaneous structures and is consistently shifted outwards with increasing quantum effects. If comparison with Figure 1c,d is made, one notes that this triplet region corresponds to the FCC pair region where the smallest maximum shows up. (Such region fades away with increasing quantum effects in the instantaneous case, Figure 1d). To get a feeling of the depth of these valleys, it seems worthwhile to quote some significant results: In addition, Figure 8 contains typical isosceles g 3 (r * , s * , s * ) behaviors of the fluid and the FCC solid at the lowest (λ * B = 0.2) and the highest (λ * B = 1.9832) de Broglie wavelengths. These graphs display significant r * −slices (i.e., at r * FP and r * SP ) of the tabulated functions in the close vicinities of the corresponding first (FP) and second (SP) maxima of the equilateral correlations. The parallels between the triplets of the solid and fluid phases coexisting at equilibrium are manifest once more. Figure 9 and Table 5 contain equilateral PIMC results for the FCC and cI16 state points in the region of moderately high densities (ρ * N = 0.9, 0.925). The centroid CM3 and the instantaneous ET3 correlation results, with P = 12 for both lattices at (ρ * N = 0.925, λ * B = 0.2), are P converged (Table 2). At (ρ * N = 0.9, λ * B = 0.8), convergence for the instantaneous correlations with P = 36 is guaranteed (practical convergence already occurs with P = 24), whereas for the centroid correlations, there is still room for further improvement. Nevertheless, the centroid results obtained with P = 36 are expected to      1 There is a cI16 small bump at r * = 2.2143, g EQ CM3 = 0.01 (P = 12), 0.18 (P = 24), 0.14 (P = 36).

Triplets in the FCC and cI16 Denser Solid Structures
In Figure 9, the equilateral correlations of the FCC and cI16 state points at (ρ * N = 0.925, λ * B = 0.2) are considered within r * < 2.5. Three well-defined features can be seen in each case, and they can be put into correspondence with the results obtained at the related distances on the pair level ( Figure 2). Thus, three separated maxima arise from the triplet FCC calculations (as occurred on the coexistence line). However, four maxima arise from the triplet cI16 calculations, with the second and third forming an overlapping structure. This reminds one of the characteristic shallow split showing up past the first maximum in the g 2 (r * ) of amorphous systems [99]. Moreover, the FCC features are more pronounced than those of cI16, as was to be expected. In addition, for 1.5 < r * < 2.5, cI16 and FCC are somewhat complementary regarding the positions of their peaks. One observes the clear quantitative differences between the centroid CM3 and the instantaneous ET3 results. The patterns of the salient features shown in Table 5 for the two density-temperature conditions are fully consistent with each other and with the underlying pair information (Figure 2).
A closer inspection of the equilateral flat regions between the first and the second maxima may be worth carrying out. The following results correspond to the discretizations: (i) P = 24 at (ρ * N = 0.925, λ * B = 0.2), although P = 12 results are not significantly different; and (ii) P = 36 at (ρ * N = 0.9, λ * B = 0.8).
(a) As regards the FCC results, these regions are related to those found on the coexistence line, but now the behavior is much more extreme: the zero-ground level is effectively reached. The solid triplet flat regions arise from the combination of the role of the QHS interactions and the unavailability of space due to the variations in ρ * N and λ * B . As a result, the solid equilateral structures analyzed turn out to be simpler than their pair radial counterparts (Figure 2), which is especially true of the cI16 lattice. Use of this fact might find application to characterizing irregular solid structures and/or monitoring their formation. (See the Supplementary Materials for more information on these structures). Another observation is related to the order shown by these two lattices. FCC appears as more ordered than cI16, which is clear from the respective regularities in their pair and triplet correlation functions and also from the maximal configurational structure factors (Table 3). Therefore, under the same (ρ * N , λ * B ) conditions, FCC entropies (and free energies) [67] must be lower than those of cI16, the determination of these properties being possible via the Einstein crystal quantum technique [66,67].

Conclusions
This article has analyzed several real space triplet correlation issues in the PI-QHS system under conditions in which quantum exchange can be neglected. Triplet PI centroid and instantaneous correlations (equilateral and isosceles) in significant fluid and FCC-solid-state points have been studied. Furthermore, the positive identification of the formerly denoted bcc-qIII solid phases [67,78] with proper quantum cI16 solid phases has been achieved by utilizing information at the pair level (radial structures and maximal structure factor values). Triplet calculations have also been carried out at two cI16 state points. The results lead to the following conclusions.
(1) Fluid phase and the use of closures.

(a)
The centroid results display far more structured triplet functions than the instantaneous results. These structures tend to be shifted outwards with increasing λ * B (delocalization) and inwards with increasing ρ * N (localization).

(b)
From the comparison between PIMC and the closure results, one concludes that the role of pair correlations in shaping triplet structures is more relevant in the quantum domain than previously thought. The combined use of KS3 (for short range) plus AV3 = (KS3+JF3)/2 (beyond short range), although not exact, is found to be a useful and simple choice to understand the related main {r} triplet features, either centroid or instantaneous, of fluids with quantum diffraction effects. (c) The AV3 success appears to be linked with the fact that this closure adopts the form of a "complete" h 2 expansion truncated to first order in the density, which includes explicitly the triple-h 2 product absent from JF3. Given that along isotherms, AV3 deteriorates with increasing distances as the fluid-solid coexistence is approached, improvements on AV3 may be of interest and should incorporate at least second-order density terms in the h 2 expansion. (d) The foregoing finding extends the previous results obtained in [17] for liquid para-hydrogen and liquid neon, since the current study has involved a purely repulsive interparticle potential. Therefore, applications of an improved AV3 (supplemented with KS3 as said above) might be expected to provide reliable pictures of what is behind triplet correlations in fluid helium over a wide range of conditions [4,15].
(a) There is a close correspondence between the positions of the main structural features at short range of both phases at equilibrium, not only at the pair level but also at the triplet level. Such a phase correspondence between triplet positions appears in both the equilateral and the isosceles correlations. These are clear signs of the system readiness to undergo the phase transition.
The triplet features are far more pronounced on the solid side. In addition, the centroid features are always sharper than those of their instantaneous counterparts (e.g., more elevated peak regions and lower valley regions for centroids).
(3) On the fluid side, the absolute maxima of the pair and the triplet-equilateral correlations, centroid and instantaneous, appear to follow empirical behaviors that depend on the quantum parameter γ = ρ * N λ * 3 B in the general form g EQ interactions dominate, this might be a further structural signature of the fluid phase on the quantum crystallization line [17,18,49], and it deserves further examination. (4) FCC and cI16 solid phases.
(a) FCC state points show a significantly higher order than their cI16 counterparts, at the same density-temperature conditions, which can be ascribed to the openness of the cI16 lattice. This is observed for the two structures, centroid and instantaneous, in all the forms computed g 2 , g 3 , S ET2 (cI16). FCC entropies must certainly be lower than their cI16 counterparts, and it is tempting to explore the relationships between the solid entropy and the values of the quantum structure factors in future work.
Within the short-medium range of distances (i.e., 1 < r * < 2.5) the equilateral functions adopt shapes simpler than the pair radial functions. This effect turns out to be much more remarkable for cI16 state points, which show quite a convoluted peak/valley behavior. Accordingly, for the purposes of monitoring the onset of crystallization and/or characterizing irregular solid phases in general, triplet centroid information may advantageously complement the usual pair level information. (c) PIMC calculations of solid centroid triplet structures converge slowly with increasing quantum effects, which contrasts with the more rapid convergence of the centroid pair calculations. This fact should be kept in mind when studying centroid triplets in high-density solid phases at low temperatures.
(5) Finally, one must dwell a little more on the (mechanically stable) QHS-cI16 phase that, as is shown in this work and [67], arises for lower densities than in the classical case. Once the question of its appearance from the PIMC evolution of perfect BCC lattices has been settled, there are no symmetry problems related to the calculations of cI16 free energies [67]. The selection of an appropriate reference system (Einstein crystal [66,67]) can be well defined now [80], and the way to computational studies of stability is open. Although there is every reason for believing that, as in the classical case [80], quantum-cI16 is metastable with respect to FCC (or to HCP = hexagonal close-packed) at low temperatures, due to the cI16 higher energies and pressures [67], the assessment of such behavior seems highly valuable. In this connection, one notes the potential QHS-cI16 relations to (i) the high-pressure solid-solid transitions in alkali metals at low temperatures and (ii) the special responses to external fields that these solid structures, which are less tight than FCC (or HCP), might exhibit.
There is work in progress to tackle the issues raised above and to identify some essential facts associated with quantum condensed matter triplets in the real and the Fourier spaces. Funding: This research received no external funding.

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