Condensation and Crystal Nucleation in a Lattice Gas with a Realistic Phase Diagram

We reconsider model II of Orban et al. (J. Chem. Phys. 1968, 49, 1778–1783), a two-dimensional lattice-gas system featuring a crystalline phase and two distinct fluid phases (liquid and vapor). In this system, a particle prevents other particles from occupying sites up to third neighbors on the square lattice, while attracting (with decreasing strength) particles sitting at fourth- or fifth-neighbor sites. To make the model more realistic, we assume a finite repulsion at third-neighbor distance, with the result that a second crystalline phase appears at higher pressures. However, the similarity with real-world substances is only partial: Upon closer inspection, the alleged liquid–vapor transition turns out to be a continuous (albeit sharp) crossover, even near the putative triple point. Closer to the standard picture is instead the freezing transition, as we show by computing the free-energy barrier relative to crystal nucleation from the “liquid”.


Introduction
An age-old question in statistical physics is to what extent particles residing on the nodes of a regular lattice can reproduce the emergent properties of a continuous many-body system, such as, e.g., those encoded in the phase diagram [1][2][3][4][5][6]. Leaving aside systems such as the Potts lattice gas [7,8] or the mixture of hard hexagons and points [9], for which their phase diagrams also recall that of a simple fluid, we restrict our discussion to onecomponent lattice gases. This class of models has been recently revitalized by a series of computational studies [10][11][12][13][14][15][16][17] aimed at establishing how the order-disorder transition of hard-core lattice particles depends on the range of the forbidden region and whether the order of the phase transition can be anticipated by symmetry considerations. Systems of asymmetric hard-core particles have also been investigated [18][19][20], but here we focus on isotropic interactions. As a general rule, the more extended the range of exclusion around a particle, the more "conventional" the melting behavior. While the solid-liquid transition is dominated by strong short-range repulsive forces, a sufficiently long-ranged attraction is needed to promote a liquid-vapor transition [5,21].
Without wanting to make an exhaustive analysis of the problem, we fix our attention at model II of [2], probably the simplest instance of a lattice-gas system with three phases (solid, liquid, and vapor). With interactions extending up to fifth-neighbor sites on the square lattice, this model experiences the miracle of a phase diagram of standard type with a limited number of ingredients. Furthermore, the model can be refined (see Section 2) in such a manner as to induce another crystalline phase at higher pressures, thus making it even more appealing. All these conclusions are drawn from transfer-matrix calculations, which, however, are only feasible for lattice strips of relatively small width. Hence, a supplement of analysis is needed to conclude that the indications of the transfer matrix are genuine, i.e., that it reflects the transition behavior of a simple fluid. To this aim, we ran grand-canonical Monte Carlo simulations across the purported liquid-vapor coexistence line, not far away from the putative triple point, in order to see whether the peak of compressibility grows with lattice size, as dictated by the theory of finite-size scaling. As for the freezing transition, the litmus test will be to find a strong hysteresis from either side of the coexistence line and/or hints that the transition is thermally activated.
The paper is organized as follows. In the next Section, we describe the model and the method employed. In Section 3, we present and discuss our transfer-matrix and simulation data (additional information that would be too cumbersome to be included in the body of the paper is placed in Appendixes A and B). The conclusions follow in Section 4.

Model and Method
We study a lattice-gas model on the square lattice, with a spherically symmetric interaction extending up to fifth-neighbor sites. Calling c i (0 or 1) the occupation number of the ith site (i = 1, . . . , N s ) and r n the nth-neighbor distance, the Hamiltonian of the system reads H = ∑ i<j u(r ij )c i c j with the following potential.
In the above equation, > 0 is an arbitrary energy unit. This model departs only slightly from model II of Orban, Van Craen, and Bellemans [2] (hence, the name "modified OVB" or MOVB model), the only difference lying in the extension of the core: The originally infinite repulsion at third-neighbor distance is replaced in the MOVB model with a finite-strength repulsion. Accordingly, denoting the lattice step with a, the maximum value of the particle-number density changes from (1/5)a −2 (for the square crystal in the left panel of Figure 1) to (1/4)a −2 . Interestingly, the square crystal with density (1/4)a −2 is not the unique configuration with this density, since a one-step shift of a line of particles would keep the density unchanged (while giving rise to a different system configuration). Among the infinite number of close-packed configurations, the one with minimum energy is the centered-rectangular (c-ret) crystal represented in the right panel of Figure 1: only this crystal ensures an optimum of four particles at a distance r 4 apart from any given particle. For lower pressures, however, the stable solid at zero temperature (T = 0) is the square crystal in Figure 1 left panel, which holds the minimum energy content among all configurations. To locate the transition between the two crystals, one simply observes that for T = 0, the energy and particle numbers of the square crystal are E = 2u(r 4 )N s /5 and N = N s /5, respectively, whereas in the c-ret crystal is E = (u(r 3 ) + 2u(r 4 ))N s /4 and N = N s /4. By comparing the grand potentials, we see that the square crystal overcomes in stability the c-ret crystal for chemical potentials µ lower than 5u(r 3 ) + 2u(r 4 ) = 4.1 , corresponding to a reduced pressure Pa 2 / = (µ − 2u(r 4 ))/(5 ) = 1.3. A similar comparison between the square crystal and the T = 0 vapor (i.e., the empty lattice) sets the corresponding transition at µ = 2u(r 4 ) = −2.4 or P = 0.
To sketch the complete phase diagram of the MOVB model, we use the transfer-matrix method (see, e.g., [1]), which computes the exact pressure as a function of T and µ for the system defined on a lattice "strip" L × ∞, being finite in the row direction and infinite in the other (with periodic conditions at the boundaries of a row). Due to the rather long range interaction, the basic lattice unit for the definition of the transfer matrix consists of two rows, implying that the size of the transfer matrix equates the number of states of 2L sites. In turn, the pressure is provided in terms of the dominant eigenvalue λ 1 of the transfer matrix as follows: where k B is the Boltzmann constant. The quantity λ 1 in (2) should be evaluated numerically.
In practice, we can take advantage of a few symmetries to reduce the size of the matrix while keeping the maximum eigenvalue unchanged (see Appendix A). Once P has been determined, number density ρ and isothermal compressibility K T are obtained from the following formulae.
The peaks of ∂ρ/∂µ or K T as a function of µ at fixed T may be taken as indication of singularities in the thermodynamic limit (see Section 3.1), but only if the peak height scales as a specific power of L.
In order that both crystalline phases of the MOVB model fit into the lattice, L should be a multiple of 10. It turns out that the only viable case is L = 10, since for L = 20, the transfer matrix is large. For the OVB model, L should be a multiple of five, and we are able to treat lattice strips with up to 20 sites in a row. To assess and, where necessary, strengthen our transfer-matrix predictions, we carry out grand-canonical Monte Carlo (MC) simulations of the MOVB model (Section 3.2) with single-site moves: At each MC step, we attempt to flip the occupancy of a randomly chosen site; then, the move is accepted or rejected in accordance with the Metropolis criterion. Each simulation run starts from a typical equilibrium configuration of the system at a nearby state point; after long equilibration, we generate a trajectory that is a few million cycles long-one MC cycle consisting of N s trial moves. By dividing the production run in large blocks, statistical errors are estimated as sample standard deviations of block averages. In addition to density and energy, we compute isothermal compressibility through the fluctuations of particle number, according to the following well-known formula: where N = ∑ i c i is the current particle number, · · · denotes a grand-canonical average, ρa 2 = N /N s is the reduced density, and δN = N − N . We also monitor the density histogram, P(ρ) = δ N ,ρV (denoting δ the Kronecker delta and V = N s a 2 the system volume), which informs on the "strength" of any transformation involving a distinct density change.
The last methodology implemented is umbrella sampling (US), which we apply to the determination of the nucleation barrier for the liquid-to-(square) crystal transition of the MOVB model so as to confirm that the onset of solid occurs by the same process of thermal activation that works in the continuum (Section 3.3). An important test will be to show that the cost of solid formation decreases with increasing liquid supersaturation. In addition to a paper on the Potts lattice gas [22], we do not know of any other investigation of crystal nucleation in a realistic lattice system. Sophisticated methods exist for the nucleation barrier of the Ising lattice gas [23], which, however, cannot be easily adapted to the solid-liquid transition of interest here.
Given a criterion to identify solid-like particles within a predominantly liquid system and choosing the size n of a solid cluster as unique reaction coordinate, the work of cluster formation (namely, the free-energy difference between the supersaturated liquid with and without a solid cluster) reads ∆Ω(n) = −k B T ln(N n /N s ) [24,25], where N n is the average number of n-clusters per configuration (see Appendix B). The maximum of ∆Ω(n) is the height of the nucleation barrier, whereas its abscissa is the critical cluster size discriminating (in probabilistic terms) between extinction and growth. However, for low to moderate supercooling/overcompression, the spontaneous occurrence of a large solid cluster in the metastable liquid is a rare event. This poses a problem of poor statistics in the estimate of N n by MC, which in a US simulation is overcome by the use of a biasing potential that couples with size n max of the largest cluster. By properly reweighing the sampled microstates, one eventually recovers the ordinary ensemble averages (more details on the technicalities of the US method in [26]). In the present study, the biasing potential is 0 in a window around the target size, while it is infinite otherwise. The main obstacle to the calculation of ∆Ω(n) is the need for identifying the largest cluster after every MC move. This complication can be largely mitigated by the use of a hybrid scheme [27,28].

Transfer-Matrix Phase Diagram of the MOVB Model
The transfer-matrix method provides an elegant as well as exact solution to the manybody problem, at least in cases where the transfer matrix is not too large. First considering the OVB model, we extend the transfer-matrix calculations in [2] to L = 20, with the purpose to validate the conclusions reached in that paper. We choose two values for inverse temperature β = (k B T) −1 , i.e., η ≡ e β = 6.5 and η = 9, the latter pretty close to the triplepoint value (η ≈ 9.5). In Figure 2, we report density ρ and reduced compressibility ρk B TK T as a function of βµ for three values of L. Of the two peaks in ρk B TK T , the one for lower µ refers to the liquid-vapor transition, whereas the other peak highlights the transition from liquid to solid. While the freezing transition is very sharp for both temperatures, the liquid-vapor transition is milder, at least for η = 6.5, where it likely corresponds to a smooth crossover.
Moving to the MOVB model, we only have data for L = 10. In Figure 3, we plot the βµ derivative of the density as a function of βµ for a number of η values between 1 and 15, along with a number of density profiles (in the inset). We see that, while the low-µ regime is nearly identical for the OVB and MOVB models, the high-µ regime is completely different, since a second crystalline phase appears in the MOVB model, heralded by the double-peak structure of the density derivative and confirmed by the density profile. The double (rather than single) peak indicates that the square crystal will not directly transform into the c-ret crystal but rather through an intermediate fluid phase.  Only results for L = 10 are available. In the main figure, the βµ derivative of the density is plotted as a function of βµ. In the inset, a few density plots are shown. Any peak of the density derivative signals a more or less steep rise in the density, which is, in turn, indicative of the possibility of a phase transition in the thermodynamic limit.
A neater picture emerges from the phase diagram. Taking the location of peaks in the density derivative as the finite-size estimate of transition points, we obtain the diagrams plotted in Figure 4 (only the purple dots in both panels were obtained through scans along constant-µ loci, see more below). In the low-µ/low-P sector, we recognize the same phase diagram of the OVB model, but it has increasingly larger deviations as µ progressively increases. Beyond a certain βµ value, the solid-liquid locus bends towards low temperatures, implying the existence of a maximum melting temperature for the square crystal. Above µ = 4.1 (or P = 1.3 a −2 ) the stable solid is the c-ret crystal. With the possible exception of very low T, the dense fluid creeps in between the two solids.  From the knowledge of the P dependence on T for fixed µ, we can derive the entropy density s and the constant-µ specific heat per unit volume c µ by the following formulae.
A sample of these quantities is shown in Figure 5. The abrupt fall of entropy on cooling occurs at the crossover from fluid to square crystal. Hence, the peak in the T-derivative of the entropy (or, alternatively, the specific-heat maximum) provides an estimate of the square crystal-to-fluid transition point. The imperfect matching in Figure 4 between the "coexistence loci" obtained from µ and T scans is a finite-size effect, which is all the more evident when peaks are not sufficiently sharp.
To sum up, the transfer-matrix treatment of the L = 10 strip is sufficient to sketch the complete phase diagram of the MOVB model. Obviously, this analysis is by no means exhaustive, since the low-temperature regime cannot be accessed for such a small L. Furthermore, nothing can be said about the exact location of the liquid-vapor critical point.

The Liquid-Vapor Transition Is in Fact a Crossover
A long-standing issue in statistical physics is how to distinguish between first-order and second-order phase transitions in a finite system. At a first-order transition, the specific heat and the compressibility exhibit delta-function singularities in the thermodynamic limit. This should be contrasted with a second-order transition, where the second-order free-energy derivatives diverge algebraically. An infinite-volume system does not anticipate a first-order transition as the transition point is approached. Finite-volume systems do instead anticipate the onset of a phase transition of any order. This feature is exploited by numerical methods, which examine the finite-size scaling (FSS) of extrema of quantities being singular in the thermodynamic limit at the transition point. In finite systems, the counterpart of these singularities are smooth peaks, the height and shape of which depend on the strength of the phase transition.
According to the theory of FSS [11,29,30], at a first-order transition, the height of the compressibility peak on a symmetric lattice increases linearly with volume V and its width at half maximum shrinks as 1/V; at a second-order transition, the peak has a slower increase in height ∼ V γ/(dν) (with γ/(dν) < 1) and a broader width ∼ V −1/(dν) (1/(dν) < 1), where γ and ν are the usual critical exponents and d is the dimensionality of space.
In light of the above arguments, we reconsider the transfer-matrix evidence for condensation and freezing in the MOVB model. A good example is provided by the data reported in Figure 2. Although actually referring to the OVB model, these data would also apply for the MOVB model-in view of the relatively low µ values. The easier case is freezing, where the delta-function-type increase in compressibility is the clear imprint of a first-order transition. Less clear is the status of the liquid-vapor transition, since both the increase in the peak and the reduction in width are rather slow, even near the triple point. To settle the question, we have carried out Monte Carlo simulations of the OVB and MOVB models for η = 9, considering L × L lattices of various sizes (up to L = 120), in order to observe how compressibility behaves near the transition from vapor to liquid as a function of L 2 . Equilibrium averages are computed over two million MC cycles. Our results, reported in Figure 6, clearly indicate that compressibility converges to a finite value with increasing L, meaning that no strict phase transition actually occurs here. A similar conclusion follows from the plot of the density distribution (see Figure 7, which refers to L = 120). We observe a smooth changeover from a vapor-like peak to a liquid-like peak, passing through a broad histogram for µ = −5.27, with no evidence for a valley between peaks, which makes us envisage the absence of a free-energy barrier also in the infinite-volume limit.
To conclude, in the (M)OVB model, liquid and vapor are only loosely separated (i.e., there is no sharp distinction between the two) even in the thermodynamic limit. This may come as a surprise, considering that even the simple Ising lattice gas has a liquid-vapor transition. However, when the range of exclusion extends beyond the central particle, the lattice-gas model actually represents a fluid-solid model. In this case, we expect the onset of sublattice order at low temperature/high pressure, which would generally be accompanied by a first-order transition (due to the symmetry breaking involved). Whether the inclusion of an attraction outside the core induces a further liquid-vapor transition probably depends on the symmetry and dimensionality of the underlying lattice and the range of the attraction. In the inset, the maximum of k B TK T is reported vs. volume on a log-log scale for both models to show that in the thermodynamic limit no phase transition is likely to occur in either of the models.

Features of Crystal Nucleation from the "Liquid"
Over the last decades, there has been considerable interest in the study of crystal nucleation from the liquid [31,32]. Describing this process in detail is of great importance for many practical applications (such as drug synthesis and production, ice formation in clouds, and the prevention of amyloid diseases), all impacted by the issue of polymorph selection during the early stages of solidification. However, nucleation is of utmost relevance primarily for fundamental reasons, since any discontinuous phase transition is triggered by the formation of a droplet of the stable phase (at least, up to precursors [33]). Here, we exploit the universality of this connection as a key to demonstrate that the freezing transition in the MOVB model is of a standard type.
In addition to the evidence provided by the compressibility, the first-order character of freezing is also evident in the hysteresis found in sequences of MC runs initiated from either side of the transition point. For instance, for η = 6.5, we have been able to overcompress the MOVB liquid up to βµ = −3.65 (well above the nominal transition point at βµ −3.87) and to expand the solid down to βµ = −4.15. Not surprisingly, we instead found no little trace of hysteresis across the liquid-vapor "transition". To characterize the resistance of the metastable liquid to conversion into solid, we calculate the free-energy cost of solid-cluster formation as a function of cluster size, using µ as a driving parameter (no attempt will be made to estimate the nucleation rate).
Preliminary to any study of crystal nucleation is the choice of a local measure of crystallinity, which could enable the distinction, in any system configuration, between solid-like and liquid-like particles. To simplify things, we consider a metastable liquid at moderate pressure so as to avoid competition between different polymorphs (in the reasonable assumption that, in the relevant range of densities, there would be no crystalline structure capable to compete in energy with the square crystal). In this case, the crystallinity criterion can be directly tuned to the target solid structure (that is, to the square crystal).
Consider a configuration of the metastable liquid. We attach the "solid-like" label to any particle forming bonds with two or more particles at distance r 4 from it, but only if at least two of these bonds are perpendicular to each other. Let 1 be a solid-like particle forming two mutually perpendicular bonds with particles 2 and 3; then, the triplet {1, 2, 3} is called a "wedge of center 1" (we may also say that particle 1 is solid-like in the given configuration if it is the center of a wedge). Two solidlike particles, say 1 and 2, will be part of the same cluster if 1 belongs to a wedge of center 2, and vice versa. With these rules established, we can (i) identify the solid-like particles present in the given configuration and (ii) enumerate the connected assemblies (clusters) of solid-like particles by the Hoshen-Kopelman algorithm [34]. A typical outcome of our clustering algorithm is illustrated in Figure 8, which refers to a metastable liquid for η = 6.5 and βµ = −3.70. We see that most solid clusters, even the smallest ones, have a distinct square-ordered structure. Occasionally, we see two particles at r 3 distance within the same cluster, which in this case can be considered as "polycrystalline". Clearly, we could have dubbed solidlike any particle at the center of two or four wedges or might have chosen a different descriptor of local order, such as a Steinhardt order parameter [35] or the like [36], but in this case, the only difference would be in the statistics of small clusters and in the critical size, with little influence on the height of the nucleation barrier [26,37].
We are now in a position to compute the free-energy cost of a solid cluster, using the US method. We take η = 6.5 and a 100 × 100 lattice. We divide the n max range in windows, i.e., [n 0 − 10, n 0 + 10] with n 0 = 10, 20, . . ., and carry out the simulations in sequence, performing 10 5 cycles for each n 0 to equilibrate the system, followed by a fivetimes longer production run (which proved sufficient to determine the nucleation barrier with enough accuracy). Then, the separate free-energy branches are vertically shifted so as to match with each other and with the cluster free-energy curve resulting from an unbiased MC simulation of the metastable liquid (in the latter, simulation only the statistics of clusters with size not larger than ≈15 turns out accurate). The final cluster free energy is shown in Figure 9 for three supersaturations (βµ = −3.75, −3.70, −3.65). For each βµ two quantities are plotted, namely β∆Ω(n) and β∆Ω * (n) + ln N s (the latter one being defined in Appendix B), which should coalesce for large n. The residual discrepancy may be ascribed to statistical uncertainties. All in all, the cluster free energy has the usual shape and dependence on the supersaturation. Only the small-n behavior is non-standard, being strongly n-dependent and in the same terms non-monotonic for all βµ. This behavior certainly reflects the peculiar definition of crystallinity adopted, which, e.g., discourages 3-clusters relative to 2-clusters and 4-clusters. It is, by the way, clear that the cluster free energy is, by its very definition, a meaningful concept only for large clusters. For large n, most solid-like particles are gathered in a single big cluster, only for entropic reasons [38]. The typical shape of the critical cluster for βµ = −3.70 can be appreciated in Figure 10. Rather than circular, the critical cluster is slightly elongated and has an irregular contour, confirming a non-trivial role in nucleation for the length/area of the cluster boundary [39].   We add a final comment on the large-n fate of the curves in Figure 9. For βµ = −3.75, we had to stop the US simulation when the extension of the biggest cluster reached the box edge so as to avoid any spurious influence from the interaction of the cluster with its periodic images. Instead, the US simulation for βµ = −3.65 was stopped when the system completely solidified, which happened abruptly through a jump in the density and the formation of a few coexisting big clusters with 300-350 solid-like particles each.

Conclusions
The phase behavior of lattice-gas systems has both similarities and differences with that of particles in the continuum. An enlightening example is provided by the OVB model [2], a lattice-gas model defined on the square lattice. In [2], this system was claimed to have the same phase diagram of an ordinary simple fluid, with a square crystal as solid phase. However, this conclusion ensued from a transfer-matrix treatment of a strip only 10 sites wide. In the analysis presented here, we extend this width to 20, but this is still insufficient to assess the nature of the transition from vapor to liquid. We have thus carried out Monte Carlo simulations of L × L lattices, with L up to 120, by which we definitely exclude that vapor and liquid are distinct phases: Their apparent coexistence locus is no more than a disorder line. The conclusion remains if we take the repulsion at third-neighbor lattice distance to be finite; in this manner, we switch from the OVB to the MOVB model. However, the latter model at least allows for the existence of a second crystalline phase at high pressure, which brings the MOVB model closer to real-life materials.
Next, we have inquired into the first stages of the transition of "liquid" into square crystal, asking whether the nucleation process would occur in the same terms as in continuous three-dimensional systems. To this aim, we introduce a notion of crystallinity tuned to the solid at hand. By computing the nucleation barrier, we find that the cluster free energy as a function of cluster size has the same shape as for ordinary supercooled liquids, except for the smallest sizes where lattice peculiarities cause a characteristic non-monotonic variation.
Our study raises the interesting question as to whether a sharp distinction between liquid and vapor in a lattice-gas model could be achieved by either increasing the range of interaction or changing the underlying lattice. This gives us the opportunity to revisit the conclusions reached in [5] for a number of triangular-lattice-gas models, for which a more careful analysis of the liquid-vapor equilibrium is planned for the next future. Other directions of research development may concern the clustering of two-dimensional lattice particles with overlapping cores (i.e., the discrete-space counterpart of the study in [40]) or the finite-size phases of particles living on the nodes of a dense polyhedral mesh (much denser than considered in [41]).

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. On the Dominant Eigenvalue of the Transfer Matrix
The calculation of the dominant eigenvalue of the transfer matrix T i j can be conducted numerically by power iteration [42]. However, the size of the transfer matrix increases exponentially with the length L of a row, eventually leading to an exceedingly long cputime-if not even to a problem of memory overflow. To partially overcome this limitation, Runnels and Combs have proposed [43,44] to make use of symmetry properties in order to replace the transfer matrix T i j with a smaller matrix τ α β without affecting the dominant eigenvalue. The idea, which we here reproduce for the reader's convenience, is to gather together in the same equivalence class α all the two-row states obtained from a given state i by a horizontal shift (that is, a cyclic permutation of the occupancies) or a reflection relative to the vertical axis of the strip. Then, the matrix will not depend on the particular state i ∈ α. Now, by noting that τ α β is a primitive matrix [43,45], it is easy to prove that any eigenvalue λ of τ α β is also an eigenvalue of T i j (while the opposite is false). Let a u vector be given such that and take v i = u α for all i ∈ α. Then, we obtain the following: Therefore, the eigenvalue λ of τ α β is also an eigenvalue of T i j and the corresponding eigenvector is v. Since all components of the dominant eigenvector of τ α β are positive (by the Perron-Frobenius theorem), also the corresponding v vector will have positive components; hence, it will be the dominant eigenvector of T i j (again, by the Perron-Frobenius theorem).
This implies that the dominant eigenvalue of τ α β coincides with the dominant eigenvalue of T i j . The number m(L) of equivalence classes (i.e., the size of τ α β ) is typically much smaller than the original number n(L) of two-row states. For the MOVB model, n(10) = 1025 and m(10) = 78. However, for L = 20, the matrix τ α β is already so large (n(20) = 1,048,577 and m(20) = 27,012) that we could not manage to store it in the memory of our computer. We have less problems with the OVB model where, due to a more extended core, the sizes of T i j and τ α β are smaller. For example, n(20) = 196,333 and m (20) = 5140 (which is still amenable).

Appendix B. Cluster Free Energy and Its Relation to the Cluster-Size Distribution
There is a simple expression for the free energy cost ∆Ω(n) of a n-sized cluster in terms of the cluster-size distribution. Assume that we have preliminarily identified, by some reasonable criterion, the solid-like particles present in the given configuration c = {c 1 , . . . , c N s }. Following Maibaum [25], let s i (c) be the size of the cluster containing the ith site (under the assumption that c i = 1 and this particle is solid-like), whereas s i (c) = 0 otherwise. Then, the number N n (c) of n-clusters reads as follows: and its thermal average is as follows.
N n = N s δ s 1 ,n n .
Assuming that there is a cluster containing site 1, the probability that it has size n is as follows.
Then, the cluster free energy (in reduced, k B T units) reads as follows: β∆Ω(n) ≡ − ln P(s 1 = n ∩ there is a solidlike particle in 1) + ln n = − ln N n N s , (A7) by taking into account the degeneracy implicit in the choice of a particular cluster particle. When n is large enough, all solid-like particles belong to a single cluster, which is also the largest cluster: s i (c) = S(c) for each site i in the cluster (denoting with S(c) the size of the largest cluster in c). Hence, ∑ N s i=1 δ s i (c),n = S(c)δ S(c),n = nδ S(c),n and N n (c) = δ S(c),n . As a result, we have the following: where β∆Ω * (n) = − ln δ S(c),n (A9) which represents the reduced free energy associated with the probability distribution of S(c). For small n, the two sides of Equation (A8) yield different numbers and, as argued by Maibaum [25], the r.h.s. may be expected to be larger than the l.h.s., due to the penalty involved in constraining S(c) to a value smaller than its typical value in the metastable liquid.