While many physical processes are non-equilibrium in nature, the theory and modeling of such phenomena lag behind theoretical treatments of equilibrium systems. The diversity of powerful theoretical tools available to describe equilibrium systems has inspired strategies that map non-equilibrium systems onto equivalent equilibrium analogs so that interrogation with standard statistical mechanical approaches is possible. In this work, we revisit the mapping from the non-equilibrium random sequential addition process onto an equilibrium multi-component mixture via the replica method, allowing for theoretical predictions of non-equilibrium structural quantities. We validate the above approach by comparing the theoretical predictions to numerical simulations of random sequential addition.

## I. INTRODUCTION

Theoretical treatment of non-equilibrium problems represents an important and formidable challenge in the modeling of physical phenomena. Driven systems^{1–7} and active matter^{8–10} are examples of non-equilibrium processes of interest in the field of material science. Such systems also display a diverse array of complex phase transitions.^{11–13} Furthermore, the biological processes relevant to life are inherently non-equilibrium.^{14,15} Non-equilibrium processes are more complex than their equilibrium analogs in that one must consider an ensemble of dynamical trajectories (instead of an ensemble of states) and the history of the trajectory is relevant.^{14,15}

Despite the significance and sheer abundance of non-equilibrium systems, there is no comprehensive theoretical framework for their modeling. By contrast, for systems in equilibrium, a broad array of statistical mechanical tools has been developed. Examples of such tools include theories of the mean-field, re-normalization group, and liquid-state closure varieties.^{16–22} In addition to the aforementioned theoretical tools, the relationships between thermodynamic variables and how such quantities relate to phase transitions are also well established.^{16–23}

In this work, we wish to leverage the vast body of work on equilibrium statistical mechanics to better understand and describe non-equilibrium systems. One path forward in this regard is to formulate a thermodynamic framework for non-equilibrium problems; several works have formulated definitions of entropy for non-equilibrium systems, for instance.^{24–26} Maximum caliber is a generalization of this idea, where the distribution of dynamic trajectories (instead of the distribution of states in equilibrium) is inferred from the maximum entropy principle.^{27–29}

One intriguing alternative possibility for theoretically describing non-equilibrium phenomena is to discover an approximate mapping from the non-equilibrium process to an equivalent equilibrium system. For a subset of non-equilibrium problems characterized by the presence of quenched disorder (i.e., degrees of freedom not in thermal equilibrium but rather frozen in place), the replica method (also known as the replica trick) provides a path forward.^{30–32} While the replica method gained recognition with its first applications to spin glasses,^{33,34} a more complete appreciation of its power followed from the seminal work of Giorgio Parisi wherein nonphysical complications stemming from the replica trick were resolved through a phenomena called replica symmetry breaking.^{35–37} The establishment of replica symmetry breaking and the unique hierarchical structure for breaking the symmetry was a central aspect of the 2021 Nobel Prize in Physics. Despite the counterintuitive mathematics of the replica trick, which includes creating *m* copies (replicas) of the thermal degrees of freedom and then sending *m* → 0, it has enabled the solution of complex spin glass problems, in some cases yielding provably exact results.^{30–32,35–39}

In the past decade, the utility of the replica method has also been demonstrated in traditional structural glasses.^{40–48} Unlike spin glasses, structural glasses do not have any imbued quenched disorder. Nevertheless, for hard spheres, the replica method enables the identification of glassy basins from the equilibrium fluid equation of state (EOS) and the tracking of the glassy state as it approaches jamming upon compression.^{40–43} This is a remarkable demonstration of the replica method’s ability to handle what is, nominally, considered a non-equilibrium phenomenon using purely equilibrium statistical mechanics. The theory also yields the complexity, the analog of configurational entropy, which is a count of the number of glassy states, as well as a provocative prediction of an ideal glass, which is the densest amorphous glass packing and is akin to a disordered crystal. While the existence of the ideal glass is still debated, the replica method provides a comprehensive and microscopic predictive theory of structural glasses and jamming.

In this work, we leverage the replica method to approximately map the non-equilibrium random sequential addition (RSA) process^{49–58} onto an equilibrium problem. RSA is conceptually simple: one particle is added to a box in a random position and frozen in place. A second particle is added at a random position, subject to the constraint that it does not overlap with the first particle. This procedure of adding particles randomly, so long as they do not generate particle overlaps, is repeated iteratively until no more particles can be placed in the box. Since the entire history of the process influences the end result, the procedure must, practically, be repeated until statistics converge.^{59} RSA has some interesting properties, including a terminal (also called saturation) density^{49–56} beyond which the process cannot be continued and an unusual logarithmic form^{52,55,60} of the contact peak of the radial distribution function (RDF). Furthermore, unlike equilibrium hard spheres, correlations between spheres differ depending on the time point at which they were added. RSA is also a canonical example of a sequential exclusion physical process. Processes in this general family have been used to model real-world phenomena such as traffic flow and cell transport.^{61,62} Specifically, in 1D, RSA is known as Rényi’s car parking problem as it is a model for the uncooperative parking of cars in a line.^{54,60} The terminal density has also been directly calculated in 1D as well as analytical expressions for the radial distribution function and structure factor.^{60}

Previous work has recognized the utility of the replica method for RSA. In one case, a free energy for RSA was derived in the grand canonical ensemble and fit to a rational function approximation in two-dimensions in an attempt to extract the terminal density.^{63} Other work has focused on the extension of replica integral equation theory^{64–66} to RSA.^{63,67–69} We significantly extend this body of work in several important ways. First, we apply the replica method to develop an expansion for structural correlations instead of the free energy. This formulation allows for examination of hard sphere contact correlations as a function of the order in which they were added to the system, allowing for predictions on a per particle level. We also show results from one- to six-dimensions, and we provide a clear description of (and justification for) which graphical terms are included in the theory, providing a road map for further theoretical developments. Finally, the results of this work are timely with respect to recent replica theory developments in the structural glass community;^{40–48} specifically, this body of work may shed light on the apparent lack of a terminal density in replica theories (including this one) for RSA.^{67–69}

The remainder of the article is organized as follows: In Sec. II, we describe the mapping between the RSA process and an equivalent equilibrium system via the replica method. We defer the bulk of the mathematical details to the Appendix. In Sec. III, we provide computational details for the RSA simulations and compare the results of the theory to RSA numerical simulations, where we show that the agreement between theory and simulation is very good. Finally, we conclude and provide an outlook in Sec. IV.

## II. THEORY

In order to map the RSA process onto an equilibrium system, we employ the replica method—a powerful mathematical tool that allows for the thermodynamic evaluation of systems possessing quenched (frozen) and thermal (ergodic) degrees of freedom.^{30–32} Originally developed for spin glasses,^{33,34} it provides a recipe for extracting properties of the real quenched disordered system from a fictive isomorphic system whereby the quenched degrees of freedom are treated on the same footing as the thermal analogs. As a relevant example, consider the multi-step process of equilibrating hard spheres at some finite number density (*ρ*_{1}), freezing these spheres in place, and then adding and equilibrating a second “batch” of hard spheres with density *ρ*_{2} in the presence of the frozen spheres. The thermodynamic properties of this hybrid frozen/ergodic system are complicated and are not equivalent to a system of equilibrium hard spheres at density *ρ*_{1} + *ρ*_{2}. This is where the replica method enters. The isomorphic equilibrium system can be thought of as a single copy of the frozen spheres in the presence of *m* copies (replicas) of the mobile spheres.^{64–66} Within a single copy, the particles of the mobile system mutually interact; however, particles in different replicas are non-interacting. The entire system (even the originally frozen spheres) is fully thermalized. The replica method is then used to compute the relevant thermodynamic quantities at integer values of *m* and then *m* is analytically continued to zero to recover the original quenched disorder system.

The RSA process is related to (but more complicated than) the simple example above, possessing an infinite hierarchical form of quenched disorder. Each addition is a quenched disorder problem, where the particles already placed in the box are frozen and the particle that is being added is thermalized. Therefore, the equivalent equilibrium system in the thermodynamic limit is an infinite-component mixture with a tiered structure that can be imagined as follows. The first particle added is represented as a single *A* particle that interacts with *m*_{B} copies of a single *B* particle. The copies of *B* are mutually pairwise non-interacting, but they all interact with the single *A* particle. Similarly, each *B* particle gets its own *m*_{C} copies of a single *C* particle. None of the *m*_{C} × *m*_{B} copies of the *C* particles interact directly with each other, and they only directly interact with “their” *B* copy. All *C* replicas interact with the *A* particle. This structure is repeated infinitely.^{70} This interaction hierarchy is depicted graphically in Fig. 1 for *m*_{B}, *m*_{C} = 2, where the only particles that directly interact are connected by a contiguous pathway of downward-facing arrows.

Because the above system is fully thermalized and amounts to a multi-level Widom–Rowlinson mixture,^{71–73} we use liquid-state theory to compute the quantities of interest [here, the contact value of the radial distribution function *g*(*σ*), where *σ* is the hard core diameter] and the replica method is invoked to compute the values for the RSA process. As derived in Appendix A, the relationship between the RDF of the real system and the replicated system between spheres added at addition *κ* out of a total of *n* additions is

where *g*_{κ,n}(*r*|** m**) is the partial RDF between particles at level

*κ*and

*n*in the replicated tree structure that is connected by a continuously descending path,

**≡ {**

*m**m*

_{2},

*m*

_{3}, …,

*m*

_{n}} is the set of the number of replicated copies at each level, and

*g*

_{κ,n}(

*r*) is the real partial RDF between particles. This further simplifies to

Correlations between species not connected by a continuously descending path (sometimes called blocking correlations) also have a physical connection to the real RSA process, though it is more obscure. Such replica “blocking” correlations can provide the real correlations to the following example.^{64–66} In the real RSA process, we can look across separate realizations of particle additions, where we add particles identically up to some density and after which we follow different addition sequences. Particles added after the randomization step will be correlated across realizations, but only by virtue of their shared history. We do not pursue blocking correlations in this study and instead reserve their treatment for future work.

In formulating our theoretical approach, we pursue a virial expansion^{19} of the replicated mixture for *g*_{κ,n}(*r*|** m**) at contact according to standard liquid-state theory in terms of 2-, 3-, 4-body interactions.

^{74,75}For convenience, the standard liquid-state theory virial expansion is discussed in Appendix B within the context of this work. Non-zero contributions to the virial coefficients can be enumerated using graphs, as discussed in Appendix C. It can be shown that only pairs (triplets) of particles in the above equilibrium system that directly interact with each other contribute to the second (third) virial coefficient. The fourth-order virial coefficient is more complicated to compute because some of the interactions in the quartet of particles can be vanishing and still generate a finite contribution.

^{19,76–79}Beyond the fourth coefficient, the complexity grows rapidly; therefore, we truncate the expansion at the fourth order. The final fourth order expansion, derived in Appendix C, after taking the

**→ 0 limit is**

*m*where *η*_{i} ≡ *v*_{D}*σ*^{D}*ρ*_{i} and *ρ*_{i} are the total volume fraction and number density, respectively, after the *i*th RSA addition, *v*_{D} and *s*_{D} are the volume and surface area of a unit sphere in *D* dimensions, respectively, *σ* is the hard sphere diameter, and the coefficients are

where *B*_{1,2,…,n} are the standard species dependent virial coefficient from liquid state theory^{19,76} and *σ*_{a,b} are the diameters between species *a* and *b* in the virial coefficients. Numerical evaluation of the diameter derivatives with respect to the third and fourth virial coefficients are discussed in Appendixes D and E.

The terms *Q*_{3} and *Q*_{4} in Eq. (3) are just the standard terms expected in the equilibrium hard sphere fluid expansion at contact whereas the $Q\u03034$ term captures the non-equilibrium effects of RSA. An interesting corollary of this is that, at the lowest order non-equilibrium correction when *η*_{κ} = *η*_{n} the pair correlations at contact are equivalent to an equilibrium fluid at *η*_{n}. This is almost certainly a limitation of the fourth order expansion, but our analysis in Sec. III of the contact value between different additions does validate this result at low to moderate density and suggests an approximate interpretation.

We further leverage liquid-state theory to attempt to correct the truncated expansion above. For hard spheres, the analogous virial expansion underpredicts the entropy loss as a function of density (i.e., the available space is over-predicted more dramatically with increasing density).^{19} Part of the issue is that a virial series expansion is not rapidly convergent.^{19,22,80–82} The Carnahan–Starling equation of state (CS-EOS) circumvents this difficulty by approximately re-summing the terms in the virial expansion as a geometric series that can be analytically evaluated, resulting in a nearly exact expression for the contact value (and all other thermodynamic properties) of the equilibrium hard-sphere liquid phase,

where *η* = *v*_{D}*σ*^{D}*ρ* in the volume fraction, and *A* has a simple analytical form for all dimensions.^{22,43,80–82} The CS form is virtually exact over the entire fluid regime at all tested dimensions (and polydisperse mixtures) for equilibrium hard spheres.^{43,83,84} We leverage the CS relation as an approximate way to also “re-sum” higher order effects for RSA via the following ansatz:

where *Z*_{1} and *Z*_{2} are yet to be determined coefficients. We choose to set the unknown coefficients by forcing the series expansion of Eq. (8),

to agree with that in Eq. (3) for each term in density. By design, the lowest order density term from the CS relation recovers the third virial coefficient, hence our neglect of a scalar to multiply the *η*_{n} term in Eq. (8). Solving for equality of the quadratic density terms yields

and

As discussed in Sec. III, this approximate re-summed form [Eq. (8)] has a larger domain of validity than the low density expansion [Eq. (3)] alone. It also has a terminal density, by construction; however, it is far larger than the observed values. For example, our theory in 2D has a terminal density at *η* = 0.714, which is much larger than the known value of *η* = 0.547,^{55} and the predictions do not improve with dimension. For convenience, tabulated values for *Q*_{3}, *Q*_{4}, $Q\u03034$, *Z*_{1}, *Z*_{2}, and *A* from one to six dimensions are provided in Table I.

D
. | Q_{3}
. | Q_{4}
. | $Q\u03034$ . | Z_{1}
. | Z_{2}
. | A
. |
---|---|---|---|---|---|---|

1 | 1 | 1 | −1/4 | 1/4 | 1/4 | 0 |

2 | 1.5640 | 2.1289 | −0.878 08 | 0.561 43 | 0.562 01 | 0.435 99 |

3 | 2.5000 | 4.5912 | −2.885 7 | 1.154 3 | 1.190 8 | 1/2 |

4 | 4.0507 | 9.7181 | −9.220 6 | 2.276 3 | 2.156 6 | −0.050 720 |

5 | 6.6250 | 19.449 | −29.033 | 4.382 4 | 3.827 5 | −1.625 0 |

6 | 10.910 | 34.164 | −90.631 | 8.307 1 | 6.813 3 | −4.910 1 |

D
. | Q_{3}
. | Q_{4}
. | $Q\u03034$ . | Z_{1}
. | Z_{2}
. | A
. |
---|---|---|---|---|---|---|

1 | 1 | 1 | −1/4 | 1/4 | 1/4 | 0 |

2 | 1.5640 | 2.1289 | −0.878 08 | 0.561 43 | 0.562 01 | 0.435 99 |

3 | 2.5000 | 4.5912 | −2.885 7 | 1.154 3 | 1.190 8 | 1/2 |

4 | 4.0507 | 9.7181 | −9.220 6 | 2.276 3 | 2.156 6 | −0.050 720 |

5 | 6.6250 | 19.449 | −29.033 | 4.382 4 | 3.827 5 | −1.625 0 |

6 | 10.910 | 34.164 | −90.631 | 8.307 1 | 6.813 3 | −4.910 1 |

From the partial radial distribution functions at contact for RSA, it is easy to compute the total radial distribution function at contact. As elaborated upon in Appendix F, the calculation is a straightforward double integral over the continuous sequential additions,

## III. RESULTS AND DISCUSSION

In this section, we validate predictions of the replica theory of this study by comparison to direct numerical simulations of the RSA process. In particular, we compare the theory and simulation contact values for both the total RDF and the partial RDFs. The partial RDFs are grouped on the basis of the order in which they are added to the simulation box, which is equivalent to the alphabetic labels for the equivalent equilibrium system described in Sec. II. To evaluate the accuracy of the above theory, direct simulations of the RSA process for systems in six different spatial dimensions were performed. The computational cost grows rapidly with increasing dimensionality, necessitating the use of cell lists to speed up the simulations. At each density, statistics for the contact value of the total RDF were found to be well converged after roughly $O(10)$ separate realizations of a 10 000 particle simulation. The contact value was calculated using a linear extrapolation of the first three binned points after contact (*r*/*σ* = 1.0125, 1.0375, 1.0625) back to literal contact (*r*/*σ* = 1) using a bin spacing of 0.025*σ*. A larger number $O(100)$ of separate realizations were used to gather partial radial distribution functions in 3D with bins of width 0.025*σ*. Our simulations allowed up to 1 000 000 insertion attempts before we stopped simulating. As dimensionality increases, approaching the terminal density becomes more difficult; as such, we do not get as close to the terminal density in higher dimensions (though this does not inhibit the validation of the theoretical framework). Virtually exact terminal densities, up to 8D, are known from a study using a more sophisticated algorithm aimed at probing the terminal density directly.^{55}

Figure 2 compares contact values of the total RDF computed in various ways for one- to six-dimensional hard-sphere RSA processes. The x-axis in each subplot of Fig. 2 extends only up to the terminal volume fraction for the RSA process as determined by prior simulations for that dimensionality.^{55} The RSA simulation results are plotted as black open circles. For comparison, the dotted-dashed teal lines show the CS-EOS contact values for equilibrium hard spheres. At lower to intermediate packing fractions, the contact value for RSA simulations is lower than for equilibrium hard spheres due to the lack of two-body correlations in the random insertion process. However, as the density increases, the RSA process runs out of free volume more quickly because there is no correlated motion or rearrangement, which results in the contact value swiftly increasing, eventually crossing over the equilibrium hard-sphere values.

The theory derived in Sec. II follows in spirit from the liquid-state theory treatment of equilibrium hard spheres, where the approximations induced by truncating the virial expansion break down at higher packing fractions when higher-order correlations become more influential. This limitation also manifests in the uncorrected replica theory, which is plotted as short orange dashed lines in Fig. 2. As expected by analogy to equilibrium hard spheres, the agreement between theory and RSA simulation is good at low densities, but then breaks down as the packing fraction increases, with the theory underpredicting the contact value; the exact result for 1D shown in Fig. 2 also re-iterates this observation and validates the simulation results. Predictions do seem to improve with increasing dimensionality as one would suspect from the increasing ideality of equilibrium hard spheres with increasing dimension (i.e., at infinite *D* only the second virial correction is required for equilibrium hard spheres).^{85,86} The corrected (CS based) replica theory that approximately includes some of the missing higher-order terms, plotted as solid orange lines, is in better agreement with the simulation results. There is still some minor discrepancy at very high densities, possibly due to the missing RSA corrections at fifth- and higher order; discussion of other possible interpretations and future avenues for research along these lines is deferred in Sec. IV. Interestingly, it seems possible from Fig. 2 that the CS corrected (and uncorrected) replica theory may improve with increasing dimension, though further work is required to fully assess this. The CS corrected theory is probably more rapidly convergent than the uncorrected virial expansion though.

Because the uncorrected theory significantly differs from the simulated results at higher packing fractions and the corrective methodology that brings the results into alignment is somewhat *ad hoc*, we provide additional support that the theory is meaningfully capturing the physics of the RSA process by temporally decomposing the particles on the basis of the order in which they are added to the system. (In the equilibrium theory, the addition order corresponds to the “level” label described in Sec. II.) At *η* = 0.15 and *D* = 3, where the uncorrected theory, CS-corrected theory, and simulations are all in excellent agreement, we compare *g*(*σ*) of the temporal self- and cross-terms in Figs. 3(a) and 3(b) for the simulation and CS corrected theory, respectively. The particles are grouped into deciles: the first 10% of the particles added to the system, the second 10% of the particles added, and so on. In Fig. 3(c), we plot *g*(*r*) from simulation between the first decile and the *n*th decile (going across the first row of the heat maps), and in Fig. 3(d), we plot the *g*(*r*) between the 10th decile and the *κ*^{th} decile (going down the last column of the heat maps). Along with the RDFs, we show the contact value predicted by theory for each RDF as a horizontal dotted line. Across all panels of Fig. 3, we see near quantitative agreement between the simulated and theoretical contact values. The excellent agreement between theory and simulation in Fig. 3 provides strong evidence that the theoretical agreement with simulation is not fortuitous as it also captures the relatively fine-grain metric of temporally specific partial RDFs.

Our convention in Figs. 3(a) and 3(b) is that *κ* ≤ *n*, though the plot is symmetric diagonal. When *κ* < *n*, the *κ* particles were frozen when the *n* particles were added. For both the theory and the simulations, as *n* increases for any value of *κ*, the contact value also increases noticeably; that is, particles that are added later in the RSA process have stronger correlations (in a two-body sense) with frozen particles. This trend is easily understood. As the simulation box fills up, it is increasingly likely that subsequent particles will be placed in close proximity to a frozen particle. As *κ* increases (particularly for larger values of *n*), the contact value decreases, though the magnitude of the effect is much weaker. The origin of this effect is less obvious but can be imagined as follows. As the background density increases, there are a decreasing number of void spaces large enough to accommodate two particles. Therefore, at larger *κ* values, particles that are added in close succession are actually less likely to be in close proximity to each other. Note that while we can rationalize the trends in Fig. 3 by leveraging physical intuition about the non-equilibrium RSA process, the same quantitative trends are present in the theoretical predictions as well. This phenomenon also offers a physical motivation (though not a full explanation) for the observation that when *κ* = *n* the fourth order theory [Eq. (3)] predicts that the contact value is identical to an equilibrium hard sphere fluid at *η*_{n}. At low to moderate densities, the structural correlations between particles added in close succession are perturbative from equilibrium, and the overall non-equilibrium effects are stored between additions with larger temporal differences.

## IV. CONCLUSIONS AND OUTLOOK

In this work, we developed a theory to describe the non-equilibrium RSA process by mapping RSA onto an isomorphic equilibrium system via the replica method. We validate the theory by comparing to direct simulations of RSA, showing good agreement between the RDFs at contact.

This work suggests several directions for future inquiry. The first is to reduce the degree of physics lost in the current theory by the truncation in the virial expansion. For instance, some other derivation, such as a Ree-Hoover expansion, could potentially have terms with complexity that grows less rapidly with the order of the expansion.^{78,79} The second is to probe the infinite-dimension limit, where it may be possible to derive an exact expression via a full re-summation of ring diagrams^{85,86} yielding a new high *D* packing law based on RSA processes. A comparison of this scaling to the known result for the ideal glass and related jamming transition would be very interesting. Other potential extensions would be to modify the theory to account for additional complexities such as a time-dependent rate of addition in the random sequential process or particle size polydispersity. Ultimately, by building up a comprehensive theory for RSA, it might be possible to develop a comprehensive theoretical framework for all types of sequential exclusion processes.

There are also various questions that this theoretical framework invokes by way of analogy to the large body of work on the replica method as applied to structural glasses and jamming. First, the development of an expansion about the contact value would be informative to see if the unusual logarithmic form of the contact peak is recovered.^{51,55} Replica theory for structural glasses has shown remarkable success in predicting the near contact behavior—obtaining nearly quantitative predictions in supercooled soft-sphere systems.^{44,45} Second, it is known that one-step replica-symmetry breaking corresponds to the onset of configuration space fragmenting into separate basins (glassy states) in equilibrium fluids.^{40–48} By analogy, it seems reasonable that replica symmetry breaking may be required to capture the more rapid growth of the contact value in RSA as the critical density is approached. The fragmentation of configuration space in RSA (if found) would likely be due to the previously quenched particles creating localized islands of configuration space for any new thermalized addition. Such a finding would support a more fundamental link between RSA insertion saturation and regular fluid jamming. Interestingly, the possible need for replica symmetry breaking is supported by liquid state replica symmetric integral equation theory studies of RSA wherein theory was found to vastly underestimate the RDF at densities near the saturation point and seemingly avoid any singularity entirely.^{67–69} The same avoidance is found in a non-replica derived integral equation approach.^{53} We note that our CS corrected theory has a singularity [by way of the denominator in Eq. (7)], but the resultant terminal density is far too high compared with the true values. Furthermore, the singularity in the theory is by construction and not emergent.

The general strategy employed in this paper is potentially applicable to certain other non-equilibrium processes as well, though they are more complex and emergent in nature. Diffusion limited aggregation^{87} or colloidal gelation^{88} may be approximated by repeated thermalization and quenching processes. Random organization (RO), a non-equilibrium model for colloidal shearing, is another such process that also has an element of quenched disorder.^{89–92} In RO, particles are randomly placed in a box. Particles that do not overlap are fixed, and particles with overlaps are active. The active particles move randomly: they could become fixed if they move such that they do not have overlaps, or they could remain active if the move does not relieve their original overlap or if their move generates a new particle overlap. (Similarly, inactive particles can become active if an active particle moves such that it overlaps with it.) The model has interesting phase behavior as a function of particle density and move size where, depending on these variables, simulations either relax to a quiescent state where all particles are inactive or to a steady state where active particles persist indefinitely. It would be interesting to see if the replica method could be adapted to the RO model, where particles can alternate between frozen and mobile (i.e., emergent quenched disorder). This current work in combination with future efforts could further the development of a rigorous theoretical treatment of non-equilibrium statistical mechanics.

## ACKNOWLEDGMENTS

We acknowledge support from the Welch Foundation (Grant No. F-1696) and the Texas Advanced Computing Center (TACC) at the University of Texas at Austin. B.A.L. acknowledges support from the Los Alamos National Laboratory through the Darleane Christian Hoffman Distinguished Postdoctoral Fellowship.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Ryan B. Jadrich**: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Writing – original draft (supporting); Writing – review & editing (lead). **Beth A. Lindquist**: Formal analysis (supporting); Funding acquisition (supporting); Investigation (supporting); Writing – original draft (lead); Writing – review & editing (supporting). **Thomas M. Truskett**: Conceptualization (supporting); Funding acquisition (lead); Supervision (lead); Writing – original draft (supporting); Writing – review & editing (supporting).

## DATA AVAILABILITY

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

### APPENDIX A: SEQUENTIAL REPLICA TRICK

The RSA process for hard spheres can be interpreted as a sequential addition, equilibration, and positional freezing (quenching) protocol.^{63} At each step (indexed by *κ*), new particles interact with one another, and with the previous particles, via hard-core interaction potentials. For convenience of notation in this section, all energies (potentials, free energies, etc.) will be implicitly per unit of thermal energy, *k*_{B}*T*, where *k*_{B} is Boltzmann’s constant and *T* is the temperature. Also for convenience, we assign book-keeping indices to the hard-core potential between particles added at steps *κ* and *γ* as *u*_{κ,γ}(*r*|*σ*) where *r* is the center-to-center distance between the two particles and *σ* is the hard-core diameter. Thus, at step *κ*, the energy for the added particles is broken into a self-term and a contribution from the new particles interacting with all the previously added, and now frozen, particles,

where *r*_{i,κ} is the position of the *i*th particle from the *κ*th addition, and *R*_{λ} and *R*_{1:λ} are shorthand for the set of positions for addition *λ* and 1 through *λ*, respectively. Thus, the equilibrium configurational probability distribution is

where $Z\kappa (R1:\kappa \u22121)\u2261\u2211R\kappa exp[\u2212U\kappa (R\kappa |R1:\kappa \u22121)]$ is the equilibrium configurational partition function. Furthermore, we will denote an average over *P*_{κ} as ⟨⋯⟩_{κ}. Only the configurational contributions to the free energy and partition function need to be considered in this section as we seek only structural correlations.

To model the thermodynamics of a macroscopic RSA system that is self-averaging (i.e., thermodynamics does not depend on the realization of quenched disorder), we require the quench-disorder averaged configurational Helmholtz free energy for the whole RSA process,

This is a formidable quantity to compute theoretically as it is not amenable to treatment via the standard tools of equilibrium statistical mechanics. To map this problem onto the domain of equilibrium statistical mechanics, we leverage the replica trick. First, we define a replicated partition function,

where ** m** = [

*m*

_{2}, …,

*m*

_{n}] are variables that can assume any real value. Defining the corresponding replicated free energy as

it can be shown that the real free energy can be obtained from the replicated free energy via

For general ** m**, this does not simplify the calculation. However, for the special case of all positive integer

**,**

*m**Z*

_{n}(

**) is the partition function for a complex, equilibrium non-additive mixture of spheres. This is easy to see from the form of Eq. (A4): (1) every average is multiplied by the corresponding partition function of equal power, effectively converting the average to a simple summation (integration) over the particle coordinates and (2) non-additivity comes from the various powers of**

*m**m*that effectively create

*m*non-interacting clones of the newly added particles at each level in the addition sequence. However, all of the clones interact identically with previously added particles.

Ultimately, the complex mixture can be described by a branched tree encoding the hierarchical relationship among species. At level *κ* in the tree, there are *m*_{2}*m*_{3}⋯*m*_{κ} nodes that represent non-interacting copies (replicas) of the set of particles added at stage *κ* in the RSA process. Any one replica at level *κ* has a parent node (replica) at level *κ* − 1 that is common to its *m*_{κ} − 1 siblings. Parent replicas interact with all of their descendants via a hard core repulsion. More specifically, any two replicas interact via a hard core if they are related via a continuously descending (or ascending) path in the tree; otherwise, they are non-interacting.

As we seek to predict the contact value of the radial distribution function, we require a relationship between the RDFs of the final added particles (level *n*) in RSA and that of earlier analogs at some arbitrary level *κ* ≤ *n*. We obtain the relationship relating the real RDF to the replicated RDF by taking the functional derivative of Eq. (A6) with respect to $u\kappa ,n(r1\u2212r2)$, which yields

The replicated RDF is the radial distribution function between any pair of replicas at level *κ* and the final level *n* (replica symmetry has been assumed) that is connected by a continuously descending path (there are *m*_{2}*m*_{3}⋯*m*_{n} of them). After application of the derivatives and limits, one finds

The “trick” is to derive an expression for the mixture in the case of all integer ** m** and assume that this can be continued to all real values of

**.**

*m*### APPENDIX B: GENERAL DENSITY EXPANSION FOR THE MIXTURE CONTACT VALUE

Working in the canonical (NVT) ensemble, it is straightforward to show that the contact value between species *a* and *b* [*g*_{a,b}(*σ*_{a,b})] is given by

where *f* is the excess Helmholtz free energy per particle and thermal energy, *δ*_{a,b} is the Kronecker delta, *x*_{a} is the particle (mole) fraction of species *a*, *s*_{D} is the surface area of a unit D-sphere, *σ*_{a,b} = *σ*_{b,a} is the cross-diameter between species *a* and *b*, and *ρ* is the total number density.^{74,75} To obtain an expansion in density, we leverage the standard virial expansion

where *B*_{i} is the *i*th virial coefficient.^{19,74,76,77} Substituting Eq. (B2) into Eq. (B1) yields

The composition dependence of the virial coefficients is apparent from the following decomposition into the species dependent virial coefficients,

where *n*_{C} is the number of components.^{19,74,75} Equation (B3) requires the derivative of Eq. (B4) with respect to *σ*_{a,b}. Taking the derivative and collecting identical terms via the permutation symmetry of the species labels yields

The virial coefficient derivatives are related to the standard Mayer-f function *f*(*r*) of equilibrium statistical mechanics and can be expressed in a convenient graphical form.^{19,76,77} For hard spheres, *f*(*r*) depends only on the core diameter (*σ*) and is trivially related to Heaviside step function, *H*(*r*), via *f*(*r*|*σ*) = −*H*(*σ* − *r*). For succinctness, we define the additional function $f\u0303(r|\sigma )\u2261\u2202f(r|\sigma )/\u2202\sigma $, which is related to the Dirac delta function, *δ*(*r*), via $f\u0303(r|\sigma )=\u2212\delta (r\u2212\sigma )$. Graphical expressions can be defined using these two functions. The second and third order terms are fairly simple,

and

where each graph represents an integrated product of *f*(*r*|*σ*) functions (solid bond) and one $f\u0303(r|\sigma )$ (dashed bond) where the integration is over a Cartesian coordinate associated with each node specifying a specific pair of species.^{19,76,77} Specifically, the third order graph in Eq. (B8) is formally $\u221d\u222b\u222b\u222bdradrbdrcf\u0303(ra,b|\sigma a,b)f(ra,c|\sigma a,c)f(rb,c|\sigma b,c)$. Importantly, if any bond (f-function) vanishes (i.e., corresponds to non-interacting species pair) then the whole graph vanishes. This property will be key to identifying the finite species contributions from the replica tree and is particularly relevant starting at the fourth order. Specifically, the fourth order term is more complex,^{19,76,77} possessing multiple graphs with varying degrees of connectivity,

where

and

As some of the graphs possess nodes that are not directly linked, they can support cross-replica contributions to the density expansion, since they will not automatically vanish if these nodes correspond to non-interacting species.

### APPENDIX C: DENSITY EXPANSION OF THE REPLICA TREE MIXTURE CONTACT VALUE

#### 1. Graphical description of replicated mixture interactions

Applied to the replicated mixture, the summation in Eq. (B5) extends over all of the species described by the RSA replica tree. There is an infinite number of combinations to consider; however, by adopting replica symmetry (the assumption that any group of replicas with the same hierarchical relationship in the replica tree posses the same statistical correlations), the summation can be reduced to a sum over a finite number of realizable hierarchical relationships among *i* species, each with a weighting that counts the number of equivalent possibilities. The various relationships can be summarized by an abbreviated graphical notation.

The second and third order virial coefficients [Eqs. (B7) and (B8), respectively] are composed of a single fully connected graph and can, thus, not support any non-interacting species pairs. As such, the only hierarchical relationship that is compatible is all species in a single descending path (and, thus, fully interacting). The continuously descending second and third relationships are expressed as

where circled nodes are fixed and non-circled nodes are summed over from level 1 to level *n*. Both graph sets in Eq. (C1) represent a primary backbone in the replica tree, which are just one of the *m*_{2}*m*_{3}⋯*m*_{n} continuously descending paths selected by a specific choice of replicas at levels *κ* and *n*. The summed nodes, while graphically shown as having a position, can take on any position along the backbone from level 1 to *n* (this includes before *κ*). All replicas along the primary path fully interact with one another via just hard-sphere interactions. Things become a bit more complicated at the fourth virial level with the allowed graphs,

where the first is just the primary path graph relevant at all virial levels and the second is a new branched graph with a single dangling species that resides one step off of the main path. In $B4$, the node *α*_{3} can sweep from level 1 to *n* − 1 and node *α*_{4} must be at least one level lower than either *α*_{3} or *κ*, depending on which is lower (note: *α*_{3} and *α*_{4} can be permuted, and must be for a full accounting of possible relationships). Replicas on a branch do not interact with those on the primary path that comes after the branch point. Examples of graphs that do not contribute at the fourth virial level are

as they have too many vanishing pairs of interactions. Specifically, *α*_{4} and *n* have two species they do not interact with in the first and second graph types, respectively. The fourth virial coefficient can support at most one disconnect for any of the species. Each higher order virial coefficient can support one more break, allowing for more complicated relationships. Finally, graphs with branches more than one node deep are irrelevant as they vanish in the ** m** → 0 limit of RSA, as discussed below.

#### 2. Density expansion

Using the results of the previous section, we can now calculate the contact value between particles added at different times during the RSA process. Replica symmetry is assumed at every step of replication in the tree, which is equivalent to assuming that any group of replicas with the same hierarchical relationship in the replica tree posses the same statistical correlations. To highlight the contributions from the specific contributions, we “re-sum” terms in Eq. (B6) according to the various graphically described contributions in Eqs. (C1) and (C2), yielding

where (i) *Q*_{i} are defined by Eqs. (4)–(6); (ii) we have changed from number density (*ρ*) to volume fraction (*η*) and recognized the replica species densities in Eq. (B6) correspond to incremental densities added in the RSA process (hence the Δ); (iii) the factor of two in front of $Q\u03034$ accounts for permuting *α*_{3} and *α*_{4}; (iv) the factor of $m\alpha 4\u22121$ comes from the dangling *α*_{4} leaf in $B4$ of Eq. (C2) that is one removed from the primary descending path. Cases where *α*_{4} is two or more deep vanish in the limit ** m** → 0 as multiplicative factors of

*m*get accrued that are not offset by any finite value. Using the definition of the total volume fraction $\eta n\u2261\u2211\alpha =1n\Delta \eta \alpha $ further simplification yields

The result in Eq. (C5) does not assume anything about how many additions are performed or what is the size of each increment. The first three terms are expected as they correspond to the equilibrium hard sphere fluid. The fourth term is unique to the non-equilibrium sequential addition nature of the problem and thus depends on the density at additions *n* and *κ* and retains a term that cannot simply be summed up to a total density (i.e., it depends on how additions are performed). In this work, we limit our study to the constant rate addition of infinitesimal amounts characteristic of what is typically referred to as random sequential addition (though recognizing it is a subset of a family of processes). Setting Δ*η*_{α} = Δ*η*, using Δ*η* → 0, and summing (integrating) the remaining sequence dependent terms in Eq. (C5) yields

providing an exact low density expansion for the structural correlations between the sets of particles added in the RSA process. A slightly regrouped form of this equation is shown as Eq. (3).

### APPENDIX D: DIAMETER DERIVATIVES OF THE FULLY-INTERACTING COMPOSITION DEPENDENT VIRIAL COEFFICIENTS

The main text requires derivatives of the form $\u2202Ba,b,\alpha 3,\u2026,\alpha i/\u2202\sigma a,b$ evaluated at the point of all equivalent diameters. Using {⋯} to denote a set, the point of equivalent diameters formally means *σ*_{κ,γ} = *σ* for *κ*, *γ* ∈ *C*({*a*, *b*, *α*_{3}, …, *α*_{i}}) where *C*({⋯}) generates the pair combinations of the entries in an arbitrary set, and *σ* is the single desired diameter. For brevity, we will sometimes use the shorthand notation {*σ*_{κ,γ}} = *σ* to indicate the aforementioned conditions.

The first step in the derivation is to take the total diameter derivative of the mixture virial coefficient with respect to *σ*,

Under the condition that *∂σ*_{κ,γ}/*∂σ* = 1 and equivalent diameters, {*σ*_{κ,γ}} = *σ*, the quantity in Eq. (D1) is trivially related to the total diameter derivative of the analogous hard sphere virial coefficient for monodisperse hard spheres of diameter *σ*$(BiHS)$ as

This equivalence is a consequence of the definitions of partial and total derivatives and the functional form of $Ba,b,\alpha 3,\u2026,\alpha i$ and $BiHS$—which are identical apart from explicit species labels and corresponding diameter labels. In the same limit, the right hand side of Eq. (D1) yields *i*(*i* − 1)/2 identical terms, which combined with Eq. (D2) yields

Furthermore, as $BiHS$ is of the form $BiHS\u221d\sigma D(i\u22121)$, we have

### APPENDIX E: DIAMETER DERIVATIVE OF THE SINGLY-NON-INTERACTING COMPOSITION DEPENDENT FOURTH VIRIAL COEFFICIENT

As discussed in the main text, the first term that captures contributions involving unrelated species in the replica tree is that coming from the fourth virial coefficient. Specifically, we must evaluate *∂B*_{a,b,c,d}/*∂σ*_{a,b} at the point *σ*_{a,b} = *σ*_{a,c} = *σ*_{a,d} = *σ*_{b,c} = *σ*_{c,d} = *σ* and *σ*_{b,d} = 0. The interaction conditions of the replicated fluid make many of the Mayer-f graph contributions in Eq. (B9) vanish, requiring the evaluation of only

where each solid line indicates a Mayer-f function between the species, a dashed line indicates a derivative of the Mayer-f function with respect to the particle diameter. Referring to the left and right graph as $G4$ and $G5$, respectively, we can explicitly write them in a unified form as

where

*D* = 1 is a special case with the closed form results $G4=4\sigma 2$ and $G5=\u22123\sigma 2$. For the case *D* ≥ 2, we will show that Eq. (E2) can be reduced down from a 3 × *D* dimensional integral to that of a double integral, which is easily evaluated by quadrature.

To make progress, we will utilize a few simplifications. First, the derivative of the Mayer-f function for hard spheres is simply a Dirac delta function,

Second, we will convert D-dimensional Cartesian integrals to analogous D-dimensional spherical coordinate-based integrals

and

where *θ* is an arbitrary angle, *s*_{D} is the surface area of a unit *D*-dimensional sphere, and *H*(*r*) and *H*(*r*, *θ*) are arbitrary radial and polar functions, respectively.^{77} Finally, we will utilize the definition of the vector norm, or equivalently the law of cosines, to write

where here *θ* is the angle between vectors *r*_{1} and *r*_{2}. Using these simplifications, we arrive at

where

The integral in Eq. (E9) can be evaluated analytically for *σ*_{i,j} = *σ*, which after dropping species labels yields

where $F12(a,b;c;z)$ is the ordinary Hypergeometric function and Γ(*x*) is the Gamma function. One final simplification can be achieved by leveraging the properties of convolutions and Fourier transforms to arrive at the two dimensional integral,

where

and *J*_{v}(*x*) is the Bessel function of the first kind of order *v*.

### APPENDIX F: RANDOM SEQUENTIAL ADDITION STRUCTURAL HISTORY INTEGRAL

The random sequential addition (RSA) process of the main text can be viewed as a sequence of *n* steps whereby hard spheres are added to a volume *V* until reaching the total density *ρ* via increments of Δ*ρ*_{i}, where 1 ≤ *i* ≤ *n*, and are forevermore frozen in place for any subsequent additions. Defining *g*_{i,j}(*r*) as the radial distribution function between particles added during addition *i* and *j*, respectively, it is trivial to compute the total radial distribution function via the density weighted average,

We will find it useful to rewrite Eq. (F1) such that self and cross-terms are separate,

In the continuous addition limit of RSA: Δ*ρ*_{i} = Δ*ρ* = *ρ*/*n* and *n* → *∞*, the first term in Eq. (F2) vanishes, yielding

where *g*(*r*|*ρ*_{1}, *ρ*_{2}) is the radial distribution function between particles added when the density reaches *ρ*_{1} and *ρ*_{2} > *ρ*_{1} during the RSA process.

## REFERENCES

A single realization of an infinite sized system would also suffice assuming self-averaging.

In actuality, the theory is formulated as an aliquot of A particles in an infinite volume, and so on. Since the statistics are convergent, we can convert the single-particle representation into an incremental density increase.