{'fraction_non_alphanumeric': 0.027315512048192773, 'fraction_numerical': 0.01120105421686747, 'mean_word_length': 4.5432536783888136, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Measuring glass entropies with population annealing

We combine population annealing Monte Carlo and several thermodynamic integration techniques to measure equilibrium vibrational and configurational entropies in the metastable fluid regime beyond the dynamic glass transition. We obtain results for a three-dimensional binary mixture hard sphere system. Our results suggest that the configurational entropy vanishes before the equilibrium pressure diverges, which implies that an underlying thermodynamic glass transition exists for this system. The computational methods are general and can be applied to a variety of glass forming systems but are restricted to small system sizes.


I. INTRODUCTION
Configurational glasses have been an active area of research for decades, however, open questions about the nature of the glass transition remain. Broadly speaking, a glass transition corresponds to a dynamical phenomenon where, upon cooling, a fluid suddenly exhibits extremely slow dynamics and effectively behaves as an amorphous solid. The glass transition is, by convention, defined as the temperature where the viscosity reaches 10 12 Poise [1]. This temperature does not correspond to a thermodynamic transition and was chosen from a practical standpoint: fluids with viscosities much larger than 10 12 Poise can not be equilibrated within reasonable lab timescales. While the conventionally defined glass transition is somewhat arbitrarily defined, there are reasons to believe that the onset of slow dynamics may be a precursor to a true thermodynamic glass transition at some lower temperature. The nature of such a transition and whether it even exists is a central question in glass physics [2][3][4][5][6][7][8].
In this paper we apply microcanonical (NVT) population annealing Monte Carlo and thermodynamic integration to a binary hard sphere fluid in the glassy regime beyond the dynamic glass transition. Binary hard sphere fluids have been studied numerically [9][10][11][12][13][14] and experimentally [13,15] in the past but, despite the relatively large body of research, a detailed understanding of the physics beyond the dynamic glass transition remains elusive. The hard sphere potential is well-suited to theoretical and numerical studies due to its simple functional form and a binary mixture can be designed to prevent homogeneous crystallization so that it displays a robust glassy regime. In hard sphere systems, temperature is not an independent control parameter and instead the inverse dimensionless pressure, 1/Z = N k B T /P V , is the relevant thermodynamic variable [16]. In most simulations, the dimensionless pressure is an observable and the density or packing fraction, φ, is the control parameter. Therefore, the picture one should have in mind for hard sphere glasses is that as the density is increased, the dimensionless pressure increases and eventually diverges once the configuration of particles can no longer be compressed at the random close packing density, φ rcp [17].
There are many different theories that attempt to describe the physics behind glassy fluids at high density. One explanation, originally proposed by Gibbs and Di-Marzio [3,18], is that there exists a dynamic ergodicity breaking transition where phase space separates into many dynamically distinct regions at some density, φ d . In this picture, for φ < φ d , the system behaves as a fluid that will ergodically sample all of phase space. After the dynamic transition, the mixing is no longer sufficient to equilibrate and phase space breaks into different regions, each of which corresponds to a different glass state. In terms of actual degrees of freedom, particles are unable to diffuse on large length-scales and, instead, remain trapped locally in cages composed of their neighbors in which they vibrate. In this picture, the phase space volume of each glass state corresponds to the local vibrational degrees of freedom of the particles and, furthermore, each glass state has a corresponding "glass" or "vibrational" entropy resulting from these local particle vibrations. In this paper we use the term "glass state" to refer to a localized region in phase space and "particle configuration" or "configuration" to refer to a specific set of particle positions within a thermodynamic state. In this terminology, a glass state is made up of many different configurations which are accessible to one another on relatively short time scales. The physics of such a system is dictated by a "configurational entropy per particle", defined as where N g is the number of glass states. Understanding how S c decreases as the density is increased and, ultimately, when S c goes to zero is paramount to understanding whether or not a thermodynamic glass transition exists. If S c goes to zero at a "Kauzmann" density, φ K < φ rcp , then this would signal that a thermodynamic transition exists. Therefore, measuring S c and extrapolating it to zero is a fundamental way to determine the high-density equilibrium behavior of supercooled fluids. In the Adam-Gibbs picture [4], the Kauzmann transition results in the structural relaxation time increasing exponentially with 1/(T S c ). As such, another tactic to estimating the value of φ K is to make fits to dynamical arXiv:2103.13837v2 [cond-mat.stat-mech] 6 Apr 2021 quantities, such as structural relaxation times or viscosities, as a function of packing fraction in order to estimate the location of a dynamic divergence corresponding to the thermodynamic transition.
Another possible scenario is that there is no thermodynamic transition at finite Z and that φ K = φ rcp . In this picture, the configurational entropy vanishes at φ rcp and, as a result, no thermodynamic transition would occur [11,16]. In a thermal glassy system, this would be analogous to a super-cooled fluid branch ending at T = 0 without a thermodynamic glass transition occurring at some T K > 0.
Much work has been done to determine the nature of the glass transition and to estimate the values of φ K and φ rcp in hard sphere systems. Past work includes fits to dynamical functions, such as the Vogel-Fulcher-Tammann law, [12,13,15], free volume fits [11,14,19], and direct measurements of the configurational entropy [9,10]. The primary obstacle in these works is in sampling equilibrium glass states beyond the dynamic transition, where known numerical methods are not efficient for binary systems and, as of yet, there is no consensus on whether a thermodynamic glass transition exists for three-dimensional binary fluids. Recent work with continuously polydisperse hard sphere systems have been able to probe unprecedentedly high densities while remaining in statistical equilibrium [16,[20][21][22][23][24]. These studies suggest that a thermodynamic transition does exist in three dimensions. However, it is not clear that the physics of continuously polydisperse systems is the same as that of binary mixtures and, additionally, there are non-trivial issues with taking the thermodynamic limit of continuously polydisperse systems that need careful consideration.
The paper is organized as follows. We begin by briefly reviewing the binary hard sphere model and the observables of interest in Sec. II. We then describe the NVT version of the population annealing algorithm and introduce two new thermodynamic integration techniques to calculate the vibrational entropy of a glass state in Sec. III. We present the results from large-scale simulations in Sec. IV and the paper closes with a discussion in Sec. V.

II. MODEL AND OBSERVABLES
The system we study is a binary hard sphere fluid with a 50:50 mixture of particles with radius ratio 1.4:1 [11,19,25]. This system is known to be a good glass former because although its high density equilibrium state in the thermodynamic limit is two monodisperse crystals separated by a domain wall, this state is inaccessible in simulations starting from a random mixture. One of the primary observables of interest in a simple fluid is the dimensionless pressure, Z, defined as where P is the pressure, V is the volume, N is the total number of particles, k B is Boltzmann's constant, and T is the temperature. It is common to use the packing fraction, φ, as a control parameter where and r 3 = (r 3 0 + r 3 1 )/2 is the average of the cubed radii of the two species. Because there is no potential energy in this system, the temperature only sets the average kinetic energy of the particles and the remaining physics depends on the dimensionless ratio Z. Thus, without loss of generality, we set k B = T = 1. We work in the NVT ensemble where N and φ are set and Z(φ) is measured. Although free energy is the thermodynamic potential for NVT, due to the triviality of the energy, all of the equilibrium physics is contained in the entropy as a function of N and φ.
The binary fluid equation of state is well-approximated throughout the fluid phase by the phenomenological Boublík-Mansoori-Carnahan-Starling-Leland (BMCSL) equation of state [26,27] where y i are constants that depend on the polydispersity. For a 50:50 mixture of 1.4:1 size particles, y 1 = 0.0513, y 2 = 0.0237, and y 3 = 0.9251. This equation of state is very accurate when compared to numerical data [19], but it must break down at high density, as is clear since it remains finite until φ = 1. For glassy systems, we are ultimately interested in the configurational entropy, S c . There are several ways of estimating S c directly, including counting inherent structures [5] and using the Franz-Parisi potential [28]. However, the most common procedure, used here, is to measure the total entropy and vibrational entropy and then use the relation, where S and S vib are the total and vibrational entropies, respectively.
To better understand the definitions and relationships of these entropies, we start with the assumption that in the glassy regime configuration space is broken into ergodically disconnected regions, each of which corresponds to a different glass state. In the NVT ensemble each glass state ν appears with probability, where the sum is over the set of all glass states, c(φ), at packing fraction φ, and Ω ν is the configuration space volume of glass state ν. The entropy of glass state ν is given by the standard definition [1], Throughout this work all entropies are defined per particle. We define S vib as the average over the entropies of the glass states, The total entropy is given by the standard definition, Using these definitions and Eq. 5 yields an equation for the configurational entropy, which is similar in form to the Shannon entropy. If we assume that there is a finite number of glass states, N g , each of which has the same statistical weight, w ν = 1/N g , then S c reduces to the standard definition given in Eq. 1, The pressure and the total entropy obey the standard thermodynamic relation which, in terms of packing fraction, is given by Using this relation, the dimensionless pressure can be integrated with respect to the packing fraction in order to obtain the entropy, see App. A for details. Numerically measuring the vibrational entropy can also be done using thermodynamic integration but is non-trivial and is discussed below in Sec. III C.

A. Microcanonical population annealing
We simulate the binary hard sphere mixture at high density using an NVT ensemble version of population annealing Monte Carlo, first described in Ref. [19]. Population annealing (PA) is a sequential Monte Carlo method [29][30][31][32][33][34] similar to nested sampling [35][36][37] that is used to simulate systems with rough free energy landscapes. The general idea of population annealing is to initialize a large ensemble of independent simulations in an easy-to-equilibrate region of parameter space and to anneal towards a difficult-to-equilibrate region. After each annealing step, the population is resampled so as to keep the distribution in equilibrium. Initially the population can be kept in equilibrium via conventional Monte Carlo schemes, however, eventually the simulation enters a regime where the system is unable to equilibrate dynamically and resampling becomes necessary to keep the population in equilibrium.
In this work, we deal with hard spheres in the NVT ensemble. All allowed configurations of hard spheres have the same energy and the roughness of the free energy landscape in the glassy regime is entirely due to a rough entropy landscape. In the NVT version of PA, each replica in the population is independently initialized as a gas of particles at low packing fraction. An equilibrating procedure such as molecular dynamics or Markov chain Monte Carlo is then applied to each member of the population. Here we use event chain Monte Carlo (ECMC) [38,39]. This step is performed to equilibrate and decorrelate the population, however, in the glass regime, it only manages to move particles within local cages. After running Monte Carlo, the population is annealed by increasing the packing fraction, φ, following an annealing schedule, {φ 0 , φ 1 , . . . , φ f }, where φ 0 is in the low density fluid phase and φ f is the highest packing fraction simulated in the glassy regime. The physical process of annealing corresponds to decreasing the box volume but, in simulations, it is computationally simpler to increase the sphere radii while keeping the volume fixed. As we will see below, these two annealing processes result in different changes in entropy so a correction must be made when performing thermodynamic integration.
After increasing the sphere radii, a fraction of the population's configurations have overlaps and are illegal at the new density. Replicas with illegal configurations are erased or "culled" and are replaced by randomly resampling the remaining legal replicas with equal weight. The fraction of the population that is culled, , is called the "culling fraction" and is an important quantity for setting the annealing schedule and integrating the entropy. After culling and resampling, this process is then repeated at the new packing fraction and annealing continues until φ f is reached.
In this work we use a "hybrid" resampling method that is similar to that in Ref. [33] in order to reduce the statistical errors associated with the resampling process. In the PA scheme of Ref. [19], described above, the new population is resampled from the final population at the end of the annealing step. The hybrid method decreases the statistical errors associated with resampling by increasing the frequency of sampling during a single annealing step. Instead of resampling the population once at the end of the equilibration process, the population is sampled sev- Figure 1. Diagram representing one hybrid microcanonical population annealing step. The packing fraction is fixed at φt and the population is equilibrated using ECMC for several intervals. After each interval, a fraction of the population is subsampled and configurations that are legal at the next packing fraction in the annealing schedule, φt+1, are saved in a reservoir. Illegal configurations (marked black) are not saved. After equilibration is completed, the reservoir is randomly sampled to produce a new population at packing fraction φt+1. eral times and additional sweeps are performed before each sampling step. As a result, the entropy of each member of the population is taken into account several different times. This procedure is shown in Fig. 1. The scheme used in this work is to perform a total of c k sweeps of ECMC on each member of the population at annealing step k, where we define the number of sweeps as the number of particle movements divided by the number of particles. The sweeps are broken into an initial burn in of c k /2 sweeps. After each subsequent sweep, a fraction of the population is subsampled and configurations that are legal at the next packing fraction are saved into a reservoir. In this work, the subsampled population consists of R/(c k /2) * 1.5 randomly chosen replicas. Finally, at the end of the annealing schedule, R replicas are chosen at random from the reservoir to represent the population at the beginning of the next annealing step.
The subsample size was chosen with consideration of the culling fraction in order to ensure that the total number of legal configurations placed in the reservoir at the end of the annealing step would be larger than the population size, R. The factor of 1.5 in the subsampling step was added as an additional safety to ensure that the reservoir is always larger than the total population. This factor acts as a tunable parameter that determines the sampling rate of the entire population. If we choose our factor so that the reservoir is nearly exactly R, then the weight of each member of the population is effectively measured once and there will be no benefit in comparison with standard PA. If we choose a large reservoir, many times the size of R, then the weight of each replica Figure 2. The outer circle corresponds to configuration space at the initial packing fraction φt and the inner circle corresponds to configuration space after annealing to the new packing fraction φt+1. The fraction of configurations that are eliminated after annealing is an estimator for the fraction of configuration space volume that has been eliminated.
will effectively be measured many times resulting in reduced statistical errors, but at the expense of using more computer memory.
Both the hybrid and the standard versions of microcanonical PA give access to the culling fraction, i , at annealing step i, which is an estimator of how much configuration space volume contracts after an annealing step. In particular, the volume of configuration space decreases by a factor of 1 − i each annealing step, as shown in Fig.  2, and the corresponding change in entropy is given by where the ratio of packing fractions corrects for the fact that we keep the system volume fixed during annealing. By summing the changes in entropy over the entire simulation, it is possible to numerically integrate the total entropy, where S(φ 0 ) is the entropy at the initial packing fraction, see App. A for details. Due to limitations in computational resources, it was necessary to carry out multiple independent runs of PA rather than one run with a very large population. Results from independent simulations can be combined using weighted averaging [19,31] to reduce both statistical and systematic errors and also to estimate the magnitude of these errors. Given M independent runs of PA, each with population size R (m) , the weighted average O of an observable O, such as the pressure, is given by, whereÕ (m) and S (m) are the estimators of the observable and the entropy, respectively, in run m. The weighted average of the entropy itself depends on a summation over the annealing schedule and is given by a different formula, For fixed population size, R (m) = R, the weighted average of an observable becomes exact in the limit of infinitely many runs, M → ∞.

B. Event chain Monte Carlo
The population annealing equilibrating procedure that we use in this work is event chain Monte Carlo (ECMC), which is particularly efficient at sampling 2D and 3D hard sphere configurations [38][39][40][41]. In the version of ECMC used here, a particle is randomly chosen and is translated in a random direction until it collides with another particle. When a collision occurs, the moving particle is stopped and the struck particle is moved in the same direction until it collides with another particle. This process is repeated until the total distance travelled by the particles, called the chain length, is equal to a predetermined length, . When this distance is reached, the current moving particle is immediately stopped. This process is then repeated by randomly choosing a new starting particle to move in a new direction. We simulate a system with periodic boundary conditions, so it is sufficient to move particles in only the positive x, y, or z directions, which violates detailed balance but preserves global balance [38].
Event chain Monte Carlo provides a way of measuring the dimensionless pressure, as shown in Ref. [39]. Consider a single chain in the x direction. When two particles, j and k, collide, the distance between their centers projected in the x-direction is x k − x j . The "lifted" distance of an event chain, x final − x initial is defined as where the sum takes place over all collisions in a single chain. If this process is repeated, then the dimensionless pressure is given by an average of the lifted distance over all of the chains, The total entropy of the fluid is then given by the thermodynamic integral of the population-averaged dimensionless pressure, where Z r is the dimensionless pressure of replica r at packing fraction φ .

C. Vibrational entropy
There are many ways to measure configurational and vibrational entropy, see Ref. [16] for a review. In this section, we focus on measuring the vibrational entropies of glass states directly using two new techniques. We call the first technique the "shell" method, which is similar to Frenkel-Ladd thermodynamic integration [42] that has been previously used to measure vibrational entropies of glassy systems [21]. We call the second method replica thermodynamic integration (RTI), which is a new technique that integrates the entropy of individual replicas from the fluid state into the glass state.

Constraining shell integration method
For the shell method, we measure the vibrational entropy of a glass state by taking the initial position of each particle as a reference. A spherical hard shell centered at each particle's reference position constrains that particle to remain within the shell. Particles are unable to Each sphere is contained within a hard shell that it cannot penetrate. Spheres only interact with their own shells or with other spheres that enter their shell. The shells begin very large (left) and, as the spheres are dynamically evolved the Metropolis algorithm (middle), the shells are decreased in size. Each decrement in size results in a culling in the population which corresponds to a loss of vibrational entropy. The shells are contracted and the vibrational entropy is numerically integrated via population annealing until shells no longer overlap (right). When the shells no longer overlap, the remaining entropy can be calculated analytically as an ideal gas.
penetrate their own shells, but are able to freely penetrate the shells of other particles. Initially the shells are much larger than the system size and, as the simulation progresses, they are gradually shrunk. Throughout this process, the entropy lost during each decrement in shell size is summed. Eventually all of the shells become sufficiently small that they no longer overlap with each other and only contain their own particles. At this point, the particles can no longer interact with each other and the remaining entropy is simply that of each particle within its own shell, as seen in Fig. 3. A single particle constrained within a hard shell is simply a particle in a box or, equivalently, an ideal gas.
The shell vibrational entropy for a given configuration of particles, S shell ( x x x), where x x x is the list of initial particle position vectors, is given by where η is a parameter that controls the shell sizes, η f is the parameter value where no shells overlap, r r r(η) is a list of the radii of all of the shells, and K[r r r(η f )] is the sum of the ideal gas entropies of each particle in its shell. During the integration, all shells shrink at the same rate. When a shell no longer overlaps with other shells, then its integration stops, its shell radius no longer shrinks, and the sphere contained within the shell no longer contributes to the numerical integral. The remaining vibrational entropy for the sphere/shell pair is calculated analytically and contributes to the constant K. The integration continues with the remaining shells that have overlaps. When no more shells overlap, the integration is complete and where i enumerates the particles and r i shell is the final radius of the i th shell. This formula is straightforward to understand, 4π 3 (r i shell − r i sphere ) 3 is simply the free volume of each shell and 3/2 is the kinetic contribution to the entropy of the ideal gas. In general, an additional factor of −3 log(λ th ) is present, where λ th is the thermal deBroglie wavelength, but, as before, we set λ th = 1.
This method of calculating S vib constrains the system so that particles are unable to switch places and, therefore, are distinguishable. This approximation becomes more accurate as structural relaxation times become very large deep within the glassy regime. This method may also be used to measure the total entropy in the fluid regime, however, since particles are no longer localized in cages, an additional term of 1 − log(N/2) must be included in order to account for indistinguishability of the particles.
The numerical integration of S shell (φ) is performed using population annealing as follows: 1. Choose a random glass sample from an equilibrium distribution at fixed φ 2. Make R shell identical copies of the glass sample 3. Initialize shells so they are infinitely large 4. Evolve population using MD or MC 5. Anneal in shell size and integrate entropy.
The value of S shell will depend on the number of sweeps per annealing step (see Fig. 9) leading to some ambiguity in the definition of S vib and S c . The details of the parameters used for the S shell integration will be described in detail in Sec. III E.

Replica thermodynamic integration method
The second method of measuring the vibrational entropy is called replica thermodynamic integration (RTI). In abstract terms, the RTI method can be described by recasting Eq. 8 into a thermodynamic integral, where Z ν (φ ) is an equilibrium trajectory of pressures of increasing density that ends in glass state ν at packing fraction φ, w ν (φ) the weight of ν at packing fraction φ, and C is a normalizing constant. This expression is different from the thermodynamic integral for the total entropy which is given by In the fluid phase, for φ < φ d , there is only one thermodynamic state and the two expressions are identical, however, in the glassy phase, the two become distinct due to the formation of ergodically separate glass states. The RTI method can be realized in a PA simulation by integrating the entropy along the history of replica r at packing fraction φ, denoted as S r (φ). Using the ECMC dimensionless pressure along the history, where κ(φ |r, φ) is the replica index of the ancestor at packing fraction φ of replica r at packing fraction φ and S(φ 0 ) is the initial entropy at the low packing fraction φ 0 . In PA, the population is an equilibrium sample of the glass states that automatically takes into account weights during resampling. This means that we can simply average S r over the population to get the average glass entropy S RTI (φ), albeit with an incorrect constant of integration set in the low density limit, In the fluid regime, where replicas are individually in equilibrium, the two quantities are equal. In the glassy regime, S RTI > S (see Fig. 6) because S RTI only averages over the surviving, high-entropy replicas whereas S averages over all replicas.
To identify the correct constant of integration requires us to use the shell method or an equivalent approach. We define S (φ ) RTI (φ) with an additive constant so that it is equal to S shell (φ ) at packing fraction φ , S (φ ) RTI (φ) provides quasi-continuous estimate of S vib (φ) if the constant of integration is set at an appropriate point φ > φ d where S shell can be reliably measured. It is worth mentioning that S RTI measures entropy in a way where particle swaps are allowed and therefore naturally includes mixing entropy.

D. Equilibration
The equilibration of our simulations can be estimated in several ways. The first and perhaps most obvious is to use one of the intrinsic equilibration metric associated with population annealing [31]. A second method, discussed at the end of this section, is specific to glasses and depends on the configurational entropy.
It can be shown [19,31] that systematic errors in PA scale as 1/R and that a prefactor of this scaling is an "equilibration population size", ρ f , where R is the population size of the simulation and S is the total entropy estimator, and the variance is measured with respect to independent runs of PA. In practice, ρ f must be estimated from multiple runs with finite but sufficiently large population size. The intuition behind relating the variance of S to systematic errors is that when the entropy estimator has large fluctuations then independent simulations sample distinct regions of configuration space and produce different results. Furthermore, these results are, on average, biased since the correct value of an observable would be obtained from an entropically weighted average over a very large number of runs.
Although results from several independent simulations can be combined using weighted averaging, ρ f is not a suitable measure of equilibration of the weighted average of many simulations. The extension to multiple runs was introduced in Ref. [19], using the quantity ρ * f , where M is the number of simulations, S is defined in Eq. 17 and is the total population size of all the simulations, with R (m) the population size of the m th simulation. It can be shown that ρ * f /R tot is proportional to the systematic errors in observables obtained from the weighted average of multiple runs, all carried out with the same annealing schedule but with, perhaps, different population sizes.
Although ρ * f was originally calculated using bootstrap resampling [19], we found that it is preferable to estimate var(S) using a weighted variance of the entropies of the runs, . (31) This method gives nearly identical results to the bootstraps method used in Ref. [19], but is easier to calculate. We can also use the configurational entropy to estimate whether the population size of the simulations is sufficient. The change in configurational entropy determines how the number of glass states decreases as the density is increased and can be used to estimate the rate of die-off of independent population members in PA. If the glass transition is suitably sharp, then it is safe to assume that the population is in equilibrium before the dynamic transition and therefore, at φ d , the number of independent glass states is approximately equal to the population size, R. Since the transition is relatively sharp, it is also safe to assume that no new glass states will be discovered after the dynamic transition and that population members remain stuck in the same glass state as their ancestor at φ d . With this in mind, the number of independent glass states in our population for φ > φ d , N g , is This means that in order to have multiple independent glass states at packing fraction φ, the population size must satisfy where ∆S c = S c (φ) − S c (φ d ). The implication of this is that PA can only go a short way beyond the dynamic transition specific to the MCMC used to equilibrate the population.

E. Simulation details
We ran two different sets of large-scale simulations. The main set of simulations produced ensembles of glass states and measured values of S, S RTI , and Z for a quasicontinuous set of packing fractions up to packing fraction φ f in the glassy regime. The second set of simulations measured S shell for many glass samples in order to obtain an equilibrium measure of the vibrational entropy at different packing fractions and to normalize S RTI .
The parameters used in the first simulation are shown in Table I. The main simulations began at φ 0 = 0.3, where it is still possible to efficiently sample configurations by randomly placing spheres in a box using a rejection method. The entropy is normalized relative to the ideal gas at φ = 0 with T = 1 and λ th = 1, see App. A. To get equilibrium samples at higher densities, we ran population annealing Monte Carlo (PA) with Event Chain Monte Carlo (ECMC) as the equilibrating dynamics. Our ECMC simuluations had chain length equal to  Table I. Parameters for equilibrium population annealing runs: N is the number of particles, R is the population size of each run, Sweeps is the total number of ECMC sweeps per replica per run, M is the number of independent simulations, is the average culling fraction, and φ f is the highest packing fraction in the runs. For N = 60, there were two population sizes, 3 × 10 6 and 5 × 10 6 , whose simulations were combined using weighted averaging. For N = 30, the two sets of simulations were analyzed separately. a fixed fraction, 0.618, of the box size. This choice of chain length is significantly longer than that of Ref. [19]. In retrospect, the choice in Ref. [19] would have been preferable. For most of the simulations, the chain schedule is given by, The idea is to do enough chains to equilibrate in the fluid regime and partially equilibrate near the dynamic transition. At φ > 0.59, when the simulation is beyond the dynamic transition, we no longer depend on using ECMC to equilibrate the entire population. Instead, we only perform a few ECMC sweeps in an attempt to equilibrate each population member locally within its own glass state. For N = 100, the number of event chains at each annealing step was halved, but the number of annealing steps was doubled so that the number of chains per change in packing fraction was kept constant. This was done so that the N = 100 culling fraction remained moderate. We also performed 10 simulations with N = 30, R = 10 5 and ×10 the number of chains as in Eq. 34, see Table I for details.
The parameters used in the shell vibrational entropy simulations are presented in Table II. We are able to measure the vibrational entropy of individual glass states with relatively small population sizes, R shell , because the motion in configuration space is confined within a single glass state. Larger population sizes were tried as well and were found to produce numerically identical results which are not presented here. In order to equilibrate after each annealing step, we use the Metropolis algorithm with a fixed step size proposal chosen at each annealing step such that the acceptance ratio is always between 40% and 45%. The constraining shells are chosen such that their initial size is larger than the box containing the particles and are shrunk so that 10% of the population is culled at each annealing step. Error bars for the  shell method are obtained by bootstrapping over all of the shell integrated vibrational entropies.

IV. RESULTS
In this section we present results from a large-scale computational study of the binary hard sphere system using population annealing Monte Carlo to sample equilibrium glass states at high density. We begin by presenting equilibrium measurements of the dimensionless pressure and the entropy as a function of packing fraction in Sec. IV A. Following this, we present our measurements of the configurational entropy and a detailed comparison between the two methods used to estimate the vibrational entropy in Sec. IV B. We then compare the locations of the estimated jamming density, φ J , and Kauzmann transition location, φ K , in Sec. IV C. Finally, we present several metrics to assess the equilibration of the simulations in Sec. IV D.

A. Pressure and total entropy
The dimensionless pressure, Z, is shown in Fig. 4  of state. Our simulation results deviate from the phenomenological equation of state at high packing fractions after the dynamic glass transition. Given that our simulations for sizes N = 30 and N = 60 are believed to be in statistical equilibrium for values of φ past φ d and that the BMCSL equation of state is not correct for high density, as is clear by its divergence at φ = 1, the deviation is at least partially due to the proximity to the true divergence of the pressure at random close packing, see Sec. IV C.
In addition to the pressure, we have also measured the total entropy per particle, S, using ECMC, Eq. 20. For N = 60, these results can be seen in comparison to measurements from previous works in Table. III.  Table I. The red curve (top) corresponds to an N = 30 simulation with ten times the number of Monte Carlo sweeps and exhibits a sharper and slightly later dynamic transition than the other simulations.

B. Vibrational and configurational entropies
We measured the configurational entropy using two different methods to obtain the vibrational entropy: replica thermodynamic integration (RTI) and shell integration. The results for these two methods are shown in Fig.  5, where the continuous curves correspond to the RTI method and the points correspond to the shell method. The constant of integration for RTI is obtained from the shell method at φ = 0.59, that is, we use S For the N = 30 system, the two methods are consistent with each other for the entire range of the simulation while for N = 60, the two methods agree for φ 0.595. The difference between the two methods for N = 60 and φ > 0.595 may be due to an inadequate number of sweeps for the RTI method to estimate S vib .
Although equilibrium quantities such as Z and S are largely independent of the annealing schedule provided sufficiently large population sizes are used, the values of S RTI and φ d , defined roughly as the location of the shoulder of the S − S RTI curve, depend on the annealing schedule. This can be seen explicitly in Fig. 6, where the blue, green, and purple curves correspond to S − S RTI measurements using the parameters of Table I for simulations of N = 30, 60, and 100 particles, respectively. The blue, green and purple curves come from simulations with the same number of ECMC chains per unit change in packing fraction. The outlying red curve uppermost in the plot corresponds to an N = 30 simulation with ×10 the number of ECMC chains as that of the standard N = 30 simulation (see Eq. 34). Here φ d is significantly increased due to the increased ECMC equilibration. For φ > φ d the curve is nearly parallel to the other curves  and, after setting the constant of integration using the shell method, the configurational entropy obtained from this high sweep number simulation is nearly the same as obtained in the main, lower sweep number simulations. Finally, Fig. 6 illustrates an important feature of population annealing. For φ < φ d , we see that S − S RTI is nearly zero. This is because each replica is equilibrated by ECMC and has the same properties as the ensemble of replicas. On the other hand, for φ > φ d ergodicity is broken and each replica is confined to a single glass state so that the ensemble average no longer equals an average obtained from a single replica using ECMC. Nonetheless, the resampling step in PA ensures that ensemble averages represent equilibrium properties well beyond φ d , as discussed in more detail in Sec. IV D. Figure 5 shows that values of S vib obtained from RTI and the shell method give consistent results when averaged over many glass states. However, looking at the distribution of values from the two methods we can see some interesting and, as yet, not fully explained differences. Figure 7 shows the measured probability distribution functions of S shell and S RTI for N = 30 at φ = 0.60. The histogram corresponds to S shell values from 300 samples chosen randomly from the PA simulations and the curve corresponds to S RTI values from 6 × 10 7 samples. The data was normalized so that the averages of S RTI and S shell are equal at φ = 0.59, which is the same as the normalization used for S c . As seen in the plot and Table IV, both distributions exhibit the same general fea- tures, but there is a significant discrepancy between the two in the low entropy tail where S shell has significantly more weight than S RTI . Another interesting characteristic of the distributions is the presence of flat steps in the RTI distribution and a second maximum in the shell histogram at high entropies.
A different perspective on the relationship between S shell and S RTI is seen in their joint probability distribution, shown in Fig. 8. The different sets of colored points correspond to glass states at different packing fractions, the large black points correspond to the average entropies at each packing fraction, and the red line is a line of best fit for the φ = 0.59, 0.60, and 0.61 data. At high packing fractions, S shell and S RTI are strongly correlated and the resulting best fit line has slope 0.99 and intercept -0.02. At these densities, the joint distributions are approximately symmetric about the linear fit, but there are two distinct features of the joint distribution that require explanation. The first is that there are clusters of points at high entropy that have a narrow distribution of S shell values and a wide distribution of S RTI values. These clusters are most pronounced for high packing fraction. The second feature is that the low entropy tail of the distribution is skewed toward higher values of S RTI . This feature is most pronounced for low packing fraction. The high entropy clusters correspond to the rightmost plateau in the histogram shown in Fig. 7. These high entropy clusters require more investigation.
The points in the low entropy skewed part of the joint distribution that are apparent for lower packing fractions in Fig. 8 may represent configurations of particles that are at least partially fluid-like having diffusive rather than caged particle motions. While the RTI method should capture the full entropy of these configurations, the limited number of Metropolis sweeps used in the shell method may fail to fully explore the configuration space of these fluid-like configurations. On the other hand, for high packing fraction, the shell method finds high entropy glass states for which the RTI method finds substantially lower entropies. This may be due to the small number of ECMC chains per annealing step used in the RTI method.
One possible issue with S shell is that it is numerically poorly behaved at high density. When initializing the shells for a glass state, there is a possibility that two spheres will be nearly touching. The resulting shells will be only slightly larger than the sphere sizes and, accordingly, the ideal gas entropies associated with those shells will be large due to the logarithm of the volume present in Eq. 22. The numerical integration of S shell will compensate for the logarithmic term by having many more annealing steps, however, this amounts to subtracting one large number from another and is, therefore, errorprone. This numerical instability is inherent to Frenkel-Ladd techniques and becomes more problematic as the density approaches jamming.
Vibrational entropy, and thus also configurational entropy, is not uniquely defined because glass states are metastable and vibrational entropy slowly increases as the time allowed to explore configuration space increases. Figure 9 is a plot of the shell vibrational entropy as a function of sweeps for glass samples at φ = 0.59 and φ = 0.60, where each data point corresponds to an average over twenty different samples. The plot exhibits a linearlog behavior, which is consistent with the known logarithmic relaxation dynamics of configurational glasses.

C. Transition locations
The jamming density, φ J , can be estimated by making a free volume fit [43] to the dimensionless pressure at densities beyond equilibrium, We make this fit for φ > 0.61 and obtain values of d equal to 2.83 and 2.85 and φ J equal to 0.676 and 0.673 for N = 30 and N = 60, respectively. We also performed this for N = 30 with ten times the number of sweeps and obtained an estimated φ J of 0.677 with d = 2.91. Because our simulations have fallen out of equilibrium over most of the range of the fits, the measured dimensionless pressures are higher than their equilibrium values and the estimated values of φ J act as lower bounds to the random close packed density, φ rcp . Our measured values of φ J = 0.673 and 0.676 are slightly larger than those found in Ref. [11] and within error bars of those found in Ref. [19]. We can estimate the location of the Kauzmann transition by extrapolating S c → 0 with our curves obtained with the RTI method. As shown in the left panel of Fig. 10, the dS c /dφ curves are linear for φ > φ d as long as the simulation remains in equilibrium. We make a quadratic fit to S c for N = 30 with 0.58 < φ < 0.595 and N = 60 with 0.58 < φ < 0.59 and obtain φ K estimates of 0.653 and 0.649, respectively. We also performed this for N = 30 with ten times the number of event chains over the range of 0.585 < φ < 0.595 and obtained an estimated φ K of 0.663. These fits are shown as dotted lines in the right panel of Fig. 10. The fit ranges were chosen to be after the dynamic glass transition but within the equilibrated regime. The resulting extrapolations are generally consistent and support the inequality φ K < φ J ≤ φ rcp . Note that if the RTI data is normalized using a shell simulation with more sweeps, see Fig. 9, then φ K will shift to a yet lower value. Thus the simulations suggest that a thermodynamic glass transition occurs at finite pressure. This is the central physics result of this work.

D. Equilibration
As discussed in Sec. III D, one can estimate the systematic errors of a weighted average of many PA simulations using the quantity ρ * f (see Eq. 29). Following Ref. [19], we set ρ * f < 0.01R tot as the threshold for equilibration of the weighted average of the simulations. Here R tot is the total population of the combined simulations. The solid lines in Fig. 11 show ρ * f as a function of packing fraction for N = 30 (blue, lower) and N = 60 (green, upper). The equilibration threshold is ρ * f = 6 × 10 5 , 9.4 × 10 5 for N = 30, 60 , respectively. These thresholds are shown by red dots in the figure and suggest that the N = 30 simulations are equilibrated over the entire range up to φ = 0.63 while the N = 60 results fall out of equilibrium at φ = 0.60. Figure 11 shows a scatterplot of the joint distribution of the dimensionless pressure, Z, and total entropy, S, at φ = 0.60 for the three system sizes. Each point represents a single simulation included in the weighted average of observables and in the computation of ρ * f . The negative slope of the joint distribution shows that runs with higher values of S have lower values of Z, as expected. In weighted averaging, these high entropy runs are more heavily weighted so that the high entropy tail of the distribution must be well-sampled to accurately estimate Z, ρ * f and S. Furthermore, for insufficient sampling, ρ * f and S will tend to be underestimated and Z overestimated. It is clear that the high entropy tail is poorly sampled for the N = 60 and N = 100 runs but perhaps adequately sampled for N = 30. Thus, it seems likely that the range of equilibration suggested by the red dots in Fig. 11 is overly optimistic due to insufficiently many runs used in estimating ρ * f . A more conservative approach would be to extrapolate ρ * f using its nearly pure exponential behavior before the knee of the curve. An exponential fit in the range before the knee of curve is insensitive to the precise fitting range. Setting the fitted function to the equilibration threshold yields more conservative estimates that equilibration is achieved for φ 0.618 for N = 30 and φ 0.598 for N = 60.
In Fig. 11, we also plot e −N ∆Sc vs φ, where ∆S c (φ) = S c (φ) − S c (φ * ). In this expression we used the fitted value of S c for φ > 0.59, described in the previous section instead of the measured value since we believe the extrapolation of S c is more accurate than the measured values deep in the glassy regime and leads to a more conservative criterion for equilibration. The initial packing fraction φ * = 0.58 is chosen to be close to the dynamic transition, φ * ≈ φ d , where each replica in the population is assumed to become trapped in a distinct glass state. For values of φ such that e −N ∆Sc < 0.01R tot (shown in the plot by the height of the red dot), the combined simulations have sampled more than 100 equilibrium glass states. We see that this occurs over the whole range for which the ρ * f < 0.01R tot criteria is satisfied.

V. DISCUSSION
In this work we have introduced a new method to compute the configurational entropy of structural glasses and presented several new results relating to the glassy regime of binary hard sphere mixtures. Most importantly, we found new estimates of the Kauzmann density, φ K , and bounds to the random close packing density, φ rcp . The jamming density was obtained by fitting the pressure divergence and provides a lower bound on φ rcp that is slightly larger than previous estimates. The Kauzmann density was extrapolated from the configurational entropy measured in the equilibrium fluid regime beyond  the dynamic glass transition, φ > φ d . Our results suggest that the Kauzmann transition occurs prior to random close packing, φ K < φ rcp , so that a thermodynamic glass transition does indeed exist in this system.
To obtain these results we introduced two new computational methods for calculating the vibrational entropy of a configurational glass: the shell method, which is a variant of the Frenkel-Ladd method, and replica thermodynamic integration, which relies on integrating the entropy of individual replicas in population annealing starting from low density. Replica thermodynamic integration produces a continuous curve of vibrational entropy but requires an alternative method to set an additive constant.
Population annealing has been shown to be capable of equilibrating configurational glasses beyond the dynamic glass transition [19]. Here the dynamic transition is associated with the computational method, either Markov chain Monte Carlo or molecular dynamics, that drives equilibration in the algorithm at each packing fraction. The limitation in probing equilibrium properties with population annealing is related to the precipitous decline in configurational entropy as the density is increased. In population annealing, each member of the population is frozen in a glass state at φ d so that upon entering the glass regime there are R glass states, where R is the population size. As the packing fraction increases, the number or distinct glass states in the population decreases exponentially as R e N [Sc(φ)−Sc(φ d )] and collapses into a single glass state when this number reaches one. The exponential dependence on N explains why the method is restricted either to small systems or packing fractions only slightly above φ d . On the other hand, population annealing and closely related techniques such as parallel tempering are the only methods known to us for going beyond the dynamic transition of a Monte Carlo or molecular dynamics scheme acting at a fixed packing fraction. It should be noted that the decrease in glass configurations with density imposes the same limitations on parallel tempering in going beyond the dynamic transition. However, in continuously polydisperse systems for which the swap algorithm is effective, population annealing (or parallel tempering) combined with the swap algorithm would permit direct measurements of equilibrium properties somewhat beyond the dynamic transition for the swap algorithm, which is already deep within the glassy regime.

ACKNOWLEDGMENTS
The work was supported in part by NSF Grant No. DMR-1507506.

Appendix A: Entropy normalization
To find the normalization of the entropy per particle, we start with the thermodynamic definition, where the Boltzmann constant has been set to unity and V is the volume per particle. This equation can be integrated in volume to get where V I is the "ideal gas volume" with the property that V I 1. In this limit, we can make the approximation that the system begins as a binary mixture of two ideal gasses with entropy S(V I ) = log V I λ 3 th /2 where λ th is the thermal deBroglie wavelength. By setting λ th = 1, the resulting entropy at V becomes S(V ) = lim VI→∞ V VI P T dV + log(V I ) + 5 2 + log(2).
(A4) We can rewrite this in terms of packing fraction by using the relationship between φ and V , φ = 4π 3V where r 0 and r 1 are the radii of the small and large particles which, in our system, have values r 0 = 1 and r 1 = 1.4. Changing variables from V to φ gives where the log(φ I ) term cancels with the logarithmically diverging pressure in the integral. By rearranging, we get a form without explicit divergences, Typically we do not integrate from φ = 0 and instead start at a non-zero initial packing fraction. Starting at φ 0 , the total entropy in a population annealing simulation is given by whereZ is the equilibrium pressure obtained during the population annealing simulation. The normalization entropy is given by Eq. A7,