Tethered hard spheres: A bridge between the fluid and solid phases

The hard-sphere system is a foundational thermodynamic model, providing the first example of a purely entropic transition from the disordered fluid phase to an ordered solid crystalline phase.¹ Determining the location of this transition is not trivial, as these two phases are separated by a strong first-order phase transition that cannot be avoided; as a result, straightforward thermodynamic integration cannot be used to calculate the free energy difference between the phases, due to the presence of hysteresis.² As a consequence, different methods need to be developed. Hoover and Ree³ introduced the single occupancy model to more precisely explore the hard-sphere fluid–solid transition. In this model, space is partitioned into mutually exclusive regions, such as the Voronoi polyhedra constructed from a Face-Centered Cubic (FCC) lattice, and the centers of the particles are constrained to remain in one of these regions. This model was employed by Woodcock^4,5 to finally establish that FCC lattices are marginally more stable than hexagonal close packed (HCP), a remarkable achievement at the time. Meanwhile, Speedy⁶ had proposed in 1993 a closely related model where each of the particles were “tethered” to remain within a set distance from a location in the system, rather than confined to disjoint regions of space. This approach was explored for non-spherical particles by Donev et al.⁷ in 2007 although their analysis is restricted to the close-packed limit. The tethered particle model was recently reanalyzed in 2021 from the perspective of phase space allowing the development of additional methods to calculate the free energy without considerable restriction, including even in the fluid phase.⁸ These approaches are easily generalized to arbitrary dimensions and have already been used to determine the free energy of various crystal structures of hard hyperspheres in six dimensions.⁹ It is remarkable that, despite 65 years of simulation, the hard-sphere model continues to provide insight into phase transitions, and it still remains a rich testing ground for developing new free-energy sampling techniques.

The tethered hard-sphere model has frequently only been used as a means to determine the free energy of the untethered hard-sphere system; however, this work considers the thermodynamics of the tethered hard-sphere model on its own merits to gain insight into the structure of phase space for hard-sphere systems. This additional insight is possible as the tether collision rate is directly related to the geometry of the configurational phase space⁸ and thus to the entropy and the free energy of the system. As a result of this analysis, we also aim to provide general guidelines for finding and using continuous paths between the solid and fluid phases for determining the free energy of solids.

The remainder of this paper is organized as follows. In Sec. II, the tethered hard sphere model is briefly described and the details of the molecular dynamics simulations are provided. Afterward, the data for the particle–particle and the tether collision rates from the simulations are presented and reviewed. Then, in Sec. IV, the consistency of these data is examined and then used to determine the entropy of the system. With knowledge of the entropy and the collision rates, the geometry of the accessible phase of tethered hard sphere systems examined in Sec. V. In Sec. VI, a diagram of states is presented for the tethered hard sphere model, which includes a discussion of the properties of the system in each of these states. Finally, the main findings of the paper are summarized in Sec. VII, along directions for future research.

The system studied here is comprised of N = 13 500 particles in a cubic periodic volume V, which is varied in size to control the number density ρ = N/V. The particles are identical and have a diameter σ and mass m. Each is restricted to remain within a distance r_T from a “tether” point. While the tether point can be located arbitrarily, in this study, the tether point is an FCC (face centered cubic) lattice site to generate pathways from the most-stable hard-sphere crystal state. To generate measurements of this system, event-driven molecular dynamics (EDMD) simulations are performed using DynamO.¹⁰ For each state point, three configurations of spheres are initialized in a perfect FCC lattice. The velocities are randomly assigned from a Maxwell–Boltzmann distribution and then shifted and rescaled to zero the total momentum of the systems and to set the kinetic temperature exactly equal to one. An Andersen thermostat is used to thermalize the system, as ergodicity is poor at short tether lengths,⁸ and its mean-free time is adjusted to set it to 5% of the total number of events. The Andersen thermostat is adjusted every 100 N events and is typically stable after the first adjustment.

Each configuration is first equilibrated for 10³ N events before ten production runs of 10³ N events are performed to collect averages. Measurement uncertainties are estimated by taking the standard deviation of results from the three configurations and each of their ten production runs. All results are reported in reduced units using the particle diameter σ as a length scale, particle mass m as the mass scale, and temperature to provide the time scale $mσ/(kBT)$⁠. The key dimensionless parameters of the system are the reduced density ρσ³ and the reduced tether length r_T/σ; however, it is more convenient to consider the tether length reduced by distance between nearest FCC tether points, $a/σ=(ρc/ρ)1/3$⁠, where $ρcσ3=2$ is the close-packed number density for the FCC lattice.

In this section, the rates of particle–particle collisions and tether events from the EDMD simulations are analyzed. As shown previously,⁸ the tether event rates are directly related to the geometry of the accessible configurational phase space of the hard sphere system. This is briefly summarized here.

Consider the tethered hard sphere system as a single state point moving through a 3 N-dimensional vector space, which is the configurational phase space of the system. As overlaps are not possible in this system, the configurations where pairs of spheres intersect create impenetrable walls within the configurational phase space which the system state point cannot enter. The “thickness” of these walls is directly related to the hard-sphere diameter. The state point bounces off these walls continuously as each intersection with these walls is a particle–particle collision. Ergodicity is typically assumed; thus, the state point is assumed to visit phase space evenly. This ultimately implies that the geometry of this phase space, including measures such as the surface-to-volume ratio, is directly related to the rates of particle–particle events, and thus to pressure and other thermodynamic variables. The configurational phase space is considered to be “unraveled” so that it extends infinitely in all dimensions. Thus, the periodic boundary conditions imposed on the simulations cause these particle–particle walls associated with the excluded volume interactions to repeat throughout the vector space.

With this background, it is interesting to contrast the effect of the tethers. The tethers confine the system state point to remain within a fixed region in phase space. Geometrically, this region corresponds to a hyperdimensional volume formed by a union of N three-dimensional spheres each sited on the lattice site of a particle in its coordinates only and with a radius of the tether length; thus, it is a highly localized volume in phase space. Again the tether collision rates are directly related to how much of the confining tether hypersurface is accessible to the state point particle compared to the hypervolume of the accessible space. Crucially though, the geometry of the tether volume is known; thus, measurements of the tether event rate can be used to determine the hypersurface area. As this can only be reduced by particle–particle walls intruding into the tethered region, the tether rates can indirectly tell us about that geometry too. With this established, it is interesting to examine the rate of tether collisions and how it varies in the FCC tethered hard-sphere system.

The ideal tether model is a system where the particles only collide with their tethers; there are no particle–particle collisions. This corresponds to systems with tethers such that

so that the particles are unable to interact with each other. For the ideal tether model, the tether collision rate $ṄT$ is given by⁸ 

where D = 3 is the dimensionality of the system, β = 1/(k_BT), k_B is the Boltzmann constant, and T is the absolute temperature of the system. Generally, the observed tether rate is reduced below this ideal tether model as the particle–particle walls in phase space reduce the accessible area of the tether volume; however, values above it have been observed at high densities where the accessible hypervolume has been reduced sufficiently to offset reductions in the hypersurface area as shown later.

The particle–particle collision rate $Ṅ$ is directly related^11,12 to the pressure p of the system as follows:

As a basis for comparison, the pressure of hard sphere systems in the fluid phase is well described by the Carnahan–Starling equation of state¹³ 

where y = πρσ³/6 is the volume fraction of space occupied by the spheres. For the solid phase, the equation of state from Pieprzyk et al.¹⁴ is used for comparison,

where A = 0.016 22, B = 6.151, C = 3.8437 × 10⁻⁵, D = 27.72, E = 0.495 41, and ω = ρ/ρ_c.

The tether and particle–particle event rates for tethered systems over a range of densities and tether lengths are shown in Figs. 1 and 2. Care must be taken to not label the tethered system as a fluid at any point as it always contains some long-range order; however, as it is an extension of the untethered system, it is natural to explore representative densities in the untethered fluid $(ρσ3=0.1,0.8,0.92)$⁠, solid (ρσ³ = [1.0, 1.2]), and transition (ρσ³ = 0.96) densities. Generally, tether event rates decrease with increasing tether length, which may be expected as the confining tether hypersurface area to accessible volume ratio decreases, which is the source of the 1/r_T scaling in the ideal tether rate [see Eq. (2)]. All densities in the fluid range display a dip below the ideal tether rate as particle–particle interactions begin to occur, but approach the ideal rate once the tether radius is large enough (r_T/a ≳ 1). At this point, the tethered system can access the majority of the configurations contributing to the pressure, and it can be seen in Fig. 2 that the particle–particle collision rates approach a value very closely predicted by the Carnahan–Starling equation of state for hard sphere fluids [by combining Eqs. (3) and (4)]. It can be imagined that, once particles can fully interact with their neighbors at r_T/a ≈ 1, the remaining configurations are either translations of already accessible configurations or large density fluctuations that are improbable in this system, and thus phase space begins to takes on a regular structure, which is analyzed in more depth later (e.g., see Fig. 7).

For solid densities, the tether event rates in Fig. 1 again decrease as tether length is increased; however, unlike in the fluid density range, there is no increase back to the ideal tether rate line, where the event rates appear to level off at a value close to zero temporarily before dropping off the figure to zero. The system appears to be localized in phase space by the particle–particle “walls” and unable to reach the tether surface, which becomes a redundant confinement for long tethers. It is clear in Fig. 1 that particles rarely travel further than r_T/a ≈ 0.75. The particle collision event rates in Fig. 2 for solid densities also show an increase initially with increasing tether length as in fluid densities, before quickly leveling off at a value very closely predicted by the Pieprzyk equation of state for hard sphere fluids [by combining Eqs. (3) and (5)], denoted by the horizontal dashed lines in the figures.

Finally, the density ρσ³ = 0.96 is considered as it lies within the fluid–solid transition density for the untethered system, from the onset of freezing ρ_fσ³ ≈ 0.938 90 to the melting density ρ_sσ³ ≈ 1.037 15.⁸ This density, in particular, is chosen as there is a clear minimum in the tether collision rate with tether length for densities up to ρσ³ ≈ 0.96 before they “melt” at longer tether lengths (r_T/a ≳ 0.5) to the ideal tether model suggesting that the tethered “fluid” is recovered. At higher densities, this minimum continues to deepen; however, the collision rate no-longer increases back up to the ideal tether collision line, as for the density ρσ³ ≈ 1.0 in Fig. 1. This can be interpreted to mean that the system remains in the solid phase and does not quite “escape” from the solid phase to transition to a liquid, although this cannot be precisely defined as the system may yet still melt at longer tether lengths or simulation times than studied here.

Consider again the particle collision rates in Fig. 2 and the result for the same transition density ρσ³ ≈ 0.96 and the adjacent densities. For ρσ³ < 0.91, the particle collision rates clearly trend toward and level off at the values predicted by the Carnahan–Starling equation of state, indicating that these systems closely resemble fluids. For ρσ³ > 0.97, the particle collision rates clearly trend toward and level off at the values predicted by the Pieprzyk equation of state, indicating that these systems closely resemble solid hard sphere systems; however, at values of 0.91 < ρσ³ < 0.97, the trend toward the equation of state predictions is far less clear and the highlighted curve at ρσ³ = 0.96 shows this clearly, where the initial trend with increasing tether looks to be toward the solid equation of state line, before there is a sharp increase in collision rate up to the fluid equation of state prediction, and then another “jump” back and forth, presumably between the two phases. The tether lengths that these jumps happen at match up with those where the minimum in tether collision rates are observed (see Fig. 1). This further emphasizes that the system is seemingly becoming jammed and unjammed as the tether length is made longer through densities around this transition region. This suggests there is “jumping” between phases as the system is forced into different configurations that each allow more or less freedom of movement. It is clear that a transition region exists in the tether space as well as the density space thus a phase diagram might be considered provided the free energy can be calculated, and this is discussed in Sec. IV.

As discussed earlier and in our previous work,⁸ the event rates in this system are directly related to ratio of the exposed surface area associated with the event and the volume of accessible phase space. The entropy S is directly related to the tether collision rate and the particle–particle collision rate,

where Ω_T is the density of states and D is the dimensionality. This relation provides at least two routes to calculate the entropy of the tether system.⁸ The first is to integrate the particle–particle collision rate with respect to the density and hard-sphere diameter σ, while holding the tether length constant,

The excess entropy calculated using this expression $ΔST(σ)/kB$ corresponds to standard thermodynamic integration.

The other path is to integrate the tether collision rate with respect to the tether length, while holding the density and hard-sphere diameter constant. The entropy difference between a tethered system and an ideal tethered system at the same tether length and particle density can be determined from the tether event rate $ṄT$ as

The excess entropy calculated using this expression is known as the tether integration route. The difference between the entropy of the tethered hard sphere system and that of an ideal tethered system at the same density and tether length is deemed the “excess” entropy [i.e., $ΔST(ρ,rT,σ)≡ST(ρ,rT,σ)−STid(rT)$]. The excess entropy estimated using the two different methods presented above is shown in Fig. 3. If the entropy were an analytic function for all densities and tether lengths, the two methods of calculating the entropy should lead to the same results. This is generally true at low densities and tether lengths; however, the two routes diverge in the range of densities between the freezing and melting densities at longer tether lengths.

To better understand this difference, it is illuminating to compare the calculated excess entropies against theoretical predictions for the residual entropy of the untethered hard sphere system.⁸ The residual entropy is defined as the difference between the entropy of a system and the entropy of an ideal gas with the same number of particles, volume, and temperature. The excess tethered entropy is expected to approach the residual untethered entropy in the limit of long tethers,

At densities below the hard sphere freezing density, the excess entropy of the tethered systems, calculated using both methods, lies below the residual entropies of the hard sphere system. This is expected as the tethered system is close to the maximally spaced-out state, thus reducing the possibility of interactions over an unconstrained state. The excess entropies gradually increase with the tether length and closely approach the residual entropy of the untethered hard-sphere fluid. Above the freezing density, calculated $ΔST(σ)/kB$ values lie below the residual free energy of the hard sphere FCC solid (black dashed–dotted line); as tether length increases, the excess entropy approaches the residual entropy, just as before; however, calculated $ΔST(rT)/kB$ values lie above the residual free energy of the untethered hard-sphere solid.

Above the melting density, the tether collision rate drops rapidly with tether length (see Fig. 1). This is an indication that the system becomes trapped in the region of phase space around the perfect FCC lattice by the excluded volume interactions. In this situation, it could be expected that the excess entropy will continue to decrease without bound as the tether length increases. The volume of phase space accessible by the system remains constant as the tether length increases, while the phase space volume accessible to the ideal tether system continues to grow. In fact, if $ΔST(rT)/kB$ is shifted by the difference in the entropy between the ideal gas and ideal tether systems,

the results are once again in close agreement with the residual entropy of the untethered systems, as shown in Fig. 4.

In contrast, when $ΔST(σ)$ is shifted by $ln⁡ΩTid(rT)/ΩTig(ρ)$⁠, the results move further from the residual entropies of the untethered systems for systems with longer tether lengths. This does not make physical sense, as the longer the tether length, the closer the system should resemble an untethered system. For the shorter tether lengths, these values, much like the shifted values from the tether integrations, also are in close agreement with the residuals of the untethered system.

Another manner to check the consistency of the data is to compare the pressure, as determined from the free energy obtained from tether event rates. The particle–particle collision rate for tethered systems is shown in Fig. 5. The crosses are directly from the EDMD simulations. The circles are calculated from derivatives of spline fits to the entropies calculated from tether integrations.

The collision rates derived from the integration of the tether event rates are in good agreement with the simulation data, with the exception of systems with sufficiently long tethers in and near the fluid–solid coexistence region. The fact that the pressures agree outside this region implies that the free energies from tether collision rate integration and the particle–particle collision rate integration are the same with respect to density, but are only shifted by a constant.

Figure 6 shows the difference between the excess entropies calculated via the two different methods. The difference is only non-zero at densities lying in the solid phase for tether lengths of r_T/a > 0.7. There is a critical tether length r_T/a ≈ 0.55 below which the entropies calculated using both methods are consistent across all densities. One key thing to note is that the difference between the entropies is dependent only on the tether length and not on the density. This observation is consistent with the fact that the pressures from the two routes are the same outside the transition region.

The event rates allow us to more carefully explore the geometry of phase space. For a tethered system, the maximum possible accessible volume is that of an ideal system [i.e., $Ωid=4πrT3/3N$]. The inclusion of excluded volume between the particles will decrease the accessible volume of phase space and increase the bounding surface area. As the hard-sphere diameter increases, the thickness of these bodies will increase. There are two types of bounding surfaces: One corresponding to the tethers and the other corresponding to the particle excluded volumes (the constraint that particles cannot overlap). The tether collision rates are directly related to the area Σ of bounding surface due to the tether constraint to the volume of phase space accessible by the system⁸ as given,

Let us first consider systems with a fixed tether length r_T/a and observe what happens with increasing the value of σ/a. In this case, the maximum available surface area will be constant and equal to the surface area in the ideal tether model. The fraction of the surface area in the ideal tether model that is accessible by the system can be determined by

In Fig. 7, the fraction of the tether boundary that is accessible to the system as a function of tether length and density is presented. At fluid-like/low densities and large tether lengths, the fraction of the accessible surface area approaches a constant equal to the ratio of the volume of accessible phase space to the system to that of an ideal tether system [i.e., $ΩT(ρ,rT,σ)/ΩTid(rT)$]. This can be rationalized by considering the phase space of the tethered hard sphere system to be a porous material, with the allowed regions corresponding to the pores of the system and the regions of phase space disallowed by the excluded volume and tether constraint corresponding to the pore walls. For a sufficiently large sample of a uniform porous material, the pore area fraction on an exposed surface is expected to be the same as the pore volume fraction in the bulk of the material, regardless of dimensionality.

At densities well above the freezing point, the ratio of the accessible surface area in Fig. 7 no longer approaches the ratio of the accessible phase space volume at large tether lengths, instead it falls to values far below it. At these solid densities, the system is locked into a small section of phase space by the particle–particle excluded volume interactions, and there is almost no possibility of reaching the tether bound phase space boundary. Finally, consider the fraction of accessible tether bound surface area for the density in the transition region. There is a clear minimum observed in the data, similar to that observed in the tether collision rate data for lower densities; however, the area does not recover to the untethered limit and appears to show considerable instability over the tether lengths considered.

The data for tether collisions and accessible surface area suggest that as the tether length initially increases, particle–particle collisions block the system from reaching the tether surface, which means that the pathways for the system to escape the FCC region of phase space become restricted. As the tether length is increased beyond a certain value, the fraction of accessible tether bound surface area increases, which means that pathways out of the FCC configuration widen. This decrease and subsequent increase of Σ/Σ^id with tether length suggests that there is a “bottleneck” in the pathways from the region of phase space near the perfect FCC configuration and that of the disordered fluid. At these transition densities, the particle–particle interactions seem to form a distinct barrier that distinguishes these regions of phase space from each other.

Figure 8 shows the variation of the accessible ratio of tether bound phase space area with increasing density. This can be pictured as a study of how a fixed region of phase space bound by the tether constraints becomes increasingly inaccessible by the system as the hard sphere diameter is increased. The ratio Σ/Σ^id = 1 when ρσ³ = 0 and decreases monotonically with increasing for the ratio of accessible surface area as particle diameter for all tether volumes used. The rate of decrease with density increases as the tether length increases. At very short tether lengths, the presence of excluded volume interactions do not significantly decrease the accessible surface area until the density becomes fairly large; this is, in part, due to the fact that the FCC configuration is the most efficiently packed. For very long tether lengths, the rate of drop off of Σ/Σ^id with density eventually coincides with the residual free energy of the untethered hard sphere fluid (black dashed line in Fig. 8). This corresponds to the expectation that the fraction of pores on a surface of a porous material is approximately equal to the pore volume fraction when the same size is sufficiently large. For intermediate values of the tether length, the value of Σ/Σ^id falls below the residual free energy of the fluid at higher densities. These tether lengths correspond to those where there is a bottleneck in passing from the FCC to the disordered fluid regions of phase space (e.g., see ρσ³ = 0.96 in Fig. 7).

For all tether lengths, Σ/Σ^id lie below the residual free energy of the untethered fluid at densities beyond the freezing density, eventually falling to zero. This is an indication that the system becomes trapped in the region of phase space near the FCC configuration.

If the entropy was an analytic function of the tether length r_T/a and the density ρσ³ for all state points, the value of a line integral along any closed path should be zero. This implies that the entropy calculated by any path should lead to the same value. However, from the analysis of the entropy derived from the collision rates in Sec. IV, it is observed that this is not the case for the tethered hard sphere system, as the entropy determined by integrating the tether collision rate with respect to the tether length, at constant density [see Eq. (8)], can lead to a different result to that obtained by integrating the particle–particle collision rate with respect to the density (or hard sphere diameter), at a fixed tether length [see Eq. (7)]. This implies that there must be regions in the states of the tethered hard sphere system (i.e., particular values of r_T/a and ρσ³) where the entropy is not analytic.

It is well established that thermodynamic integration cannot be performed through a first-order phase transition. In particular, the integration of the pressure cannot be used directly through the freezing transition to determine the free energy of the solid phase. This has led to the development of a variety of alternate methods to determine the free energy of the solid phase, such as the Einstein crystal method,^2,15 self-referential methods,^16–18 as well as the use of the single occupancy cell model,^3–6 along with other methods. This non-analytic behavior of the untethered hard sphere system in the fluid–solid transition region is expected to extend to the tethered systems. In fact, indications of this are given by the sharp drop in the pressure with density for sufficiently large tether lengths (see Fig. 5).

The difference between the entropy calculated from tether integration and density integration is shown in Fig. 6. In Fig. 6(a), it is observed that the differences between the two methods vanish for densities below about the freezing density. At low tether lengths (where r_T/a ≲ 0.55), the two methods lead to the same result across all densities. The consistency between the two methods suggests that the entropy is analytic in these regions. In fact, it might be expected that the entropy is analytic everywhere the pressures are consistent between the two methods (see Fig. 5).

As can be seen in Fig. 6, the entropies between the tether integration and the density integration methods differ at densities above the freezing density and r_T/a ≳ 0.55. In this region, difference between the two methods is independent of density, but depends on the tether length, becoming larger with increasing r_T/a. This observation is in-line with the fact that the two methods lead to the same value pressure in this region (see Fig. 5). This implies that for a set tether length, the density integration entropies can be shifted by an amount independent of the density to align the results with the entropy calculated by tether integration. In Fig. 6(b), the difference between $ΔST(rT)/kB$ and $ΔST(σ)/kB$ at a fixed density is shown as a function of tether length. The dashed black line is the entropy difference between the ideal tether model and an ideal gas [see Eq. (10)]. As the density of the system increases, the difference between the two calculated entropies appear to approach a common curve that is independent of density; this curve differs from zero at a tether length r_T/a ≈ 0.55 and monotonically increases, asymptotically growing as 3 ln(r_T/a) for very large tether lengths.

A schematic diagram of states for a hard sphere system tethered to an FCC lattice is shown in Fig. 9. The shaded, light blue area on the left side of the diagram denotes the ideal tether state, where the particles are too far apart to interact with each other. The dark gray shaded area gives an indication of the region where we presume that the entropy is not analytic, based on where the entropy and pressure calculations depend on the method of calculation (see Figs. 3 and 5). For an untethered hard sphere system, which is related to the limit r_T → ∞, there is a fluid–solid coexistence region with a fluid density of ρ_fσ³ ≈ 0.938 90 and the melting density ρ_sσ³ ≈ 1.037 15.⁸ The non-analytic behavior of the entropy might be expected to be related to this transition and to persist as the tether length decreases. The gray region appears to vanish below a critical tether length of about r_T/a ≈ 0.55. Outside this gray area, the entropy should be analytic. This implies that, if we take a line integral around any simply closed path in this region (which excludes enclosing any portion of the gray shaded area), the results should be zero. For example, the line integral of the entropy along either of the blue curves would be zero.

If the line integral of the entropy encloses a portion of the gray shaded region, its value would be non-zero. Examples of this include any of the red curves in Fig. 9 would be non-zero. Based on the observed differences in the excess entropies from tether integration and density integration, the value of the line integral would be positive, if the path is taken in a counterclockwise manner and negative if it is taken in a clockwise manner. The more of the gray region that is enclosed by the path, the larger the magnitude of the line integral; so the magnitude of the line integral for curve 4 would be larger than that for curve 3.

To attempt relate the non-analytical region of the entropy to structural changes of the system, the Steinhardt q₆ order parameter¹⁹ is calculated using the software package Freud,²⁰ and the result is given in Fig. 10.

The Steinhardt order parameter characterizes the local organization of particles in the system. The directions of the “bond” vectors b formed between a tagged particle and its nearest neighbors are projected along a set of spherical harmonic functions Y_lm of a particular order l,

where b is a “bond” vector. This quantity is then averaged over all bonds in the system to give ⟨Q_lm⟩. The Steinhardt q_l order parameter is defined¹⁹ as a rotationally invariant combination of these coefficients,

In this work, the value l = 6 is used as it is a good discriminator between FCC and random/fluid ordering.

The non-analytic region appears to correlate with where the order parameter rapidly changes, indicating the transition from a slightly less ordered structure held together by the tether at lower densities to a slightly more order structure at higher densities. Interestingly, there is an angled shape of the non-analytical region to below even the fluid transition density for the untethered system. It appears that staying below a tether length of r_T/a ≈ 0.55 allows the avoidance of the non-analytic region entirely while performing thermodynamic integrations with respect to density. Previous uses of the tether model to calculate the free energy have been largely limited to integration with respect to the tether length, mainly to avoid the need to pass through the fluid–solid phase transition.^7,21,22 Consequently, they did not consider the possibility of a continuous path between the fluid and solid phases, which avoids the non-analytical region. Earlier works with the single occupancy model for hard sphere, where the particles are confined to Voronoi cells of the lattice sites, have used density integration to determine the free energy;⁴ however, closer examination of the pressure variation with density indicated that these simulations did pass through a non-analytical region, although its effect was thought to be negligible.

In this work, the properties of systems of hard spheres that are tethered to an FCC lattice are examined. The tethered hard sphere model gives us detailed information of the structure of phase space, and more specifically the behavior of the entropy, around the point where the configuration is tethered. The phase space of the system can be considered to be a porous material, with the pores corresponding to the regions accessible to the system, and the pore walls corresponding to the disallowed regions due to overlap of spheres or violation of the tether constraint. The “pore volume” of phase space is directly related to the entropy of the system, while its “porosity” is directly related to the excess entropy of the system. The collision rates of the particles with their tethers give an indication of the pathways away from the perfect FCC configuration.

At low to moderate fluid densities, the porosity of phase space appears to be fairly uniform moving away from the perfect FCC configuration to the regions of phase space corresponding to the more disordered fluid phase. The tether collision rates gradually decrease to eventually reach the value for the ideal tether mode, while the particle–particle collision rates gradually increase from zero to the bulk fluid value.

At the higher solid densities, the FCC configuration becomes isolated from the rest of the phase space by excluded volume interactions. The system becomes “trapped” in the FCC configuration. Moving away from the perfect FCC lattice (by increasing the tether length), the particle–particle collision rate increases from zero to the values for the bulk, untethered system; however, the tether collision rate decreases exponentially with the tether length and then suddenly drops to zero.

At intermediate densities, there is a “bottleneck” separating the phase space corresponding to the FCC configuration and that of the disordered, fluid phase. The tether collision rate drops rapidly to low values as the tether length increases, before jumping up to the ideal tether rate. The particle–particle collision rate first increases to a value corresponding to solid phase, but then later switches to the fluid phase value, sometimes transitioning back and forth between the two before settling to the fluid value.

This change in the pathways out of the FCC configuration leads to the entropy developing a path dependence, and this can be represented as a region in r_T/a and ρσ³ where the entropy is non-analytic. This region extends at extremely large tether lengths from the solid–fluid coexistence region of the untethered hard sphere system and ends at a tether length r_T/a ≈ 0.55 and a density slightly lower than the density at the onset of freezing.

The direct geometric interpretation of the tether and particle–particle collision rates can be used to determine the entropy of the system: by integrating the particle–particle collision rate with respect to the system density or by integrating the tether collision rate with respect to the tether length. As long as both these paths do not cross this non-analytic region, they will lead to consistent results. If either path does cross the non-analytic region, the discrepancy between the two methods will depend on what portion of the non-analytic region is enclosed by the path.

The results of this work can be extended in many ways. An FCC lattice is explored here; however, the particular shape of the non-analytic region is expected to depend on the lattice used for the tethers. For example, if the spheres were tethered to a BCC (Body-Centered Cubic) lattice, the location of the non-analytic region would be expect to change due to the difference in the fluid–solid coexistence densities as compared to the FCC crystal structure. A more interesting case is for a system that is tethered to a “fluid” configuration, rather than to a regular lattice. This would allow the behavior of particle systems that are trapped in glassy states to be explored. In fact, Speedy⁶ used this approach to calculate entropy of hard sphere systems that were jammed in various disordered configurations in order to determine the “statistical” entropy of the glassy system. It would be interesting to determine if the non-analytic region appears for these systems, if so, where is it located, how does this depend on the particular arrangement of the tether sites, and how does this relate to the presence of a glass transition?

Also of interest is the relationship of the tethered particle model to the random pinning glass model,^23,24 where a fraction of particles are forced to maintain a fixed position. The “pinning” of these particles imposes a quenched disorder in the system, which can be analyzed within the framework of the random first order transition theory,^25–27 and a glass transition is expected to occur, signaled by the vanishing of the configurational entropy, as the fraction of pinned particles increases above a critical value.²⁷ Can the quenched disorder imposed by the tether constraint be analyzed in a similar manner?

The solid–fluid phase transition in the three-dimensional hard sphere system is well established and understood. It would be interesting to examine the analytic structure of the entropy in such two-dimensional hard disk system^28,29 and the parallel hard cube model,^30–32 where the transition has been more difficult to examine. After 65 years, there is still a rich variety of phenomena to explore in the phase transitions of hard systems, and the tether model opens a complementary line of investigation and interest.

The authors have no conflicts to disclose.

James MacKinnon: Data curation (equal); Investigation (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). Marcus N. Bannerman: Conceptualization (equal); Data curation (equal); Investigation (equal); Methodology (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). Leo Lue: Conceptualization (equal); Investigation (equal); Methodology (equal); Project administration (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.