{'fraction_non_alphanumeric': 0.03586309143426804, 'fraction_numerical': 0.013343217197924388, 'mean_word_length': 4.747318936389012, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Treating random sequential addition via the replica method

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][9][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][12][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][17][18][19][20][21][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][17][18][19][20][21][22][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][25][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][28][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][31][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 seminal work of Giorgio Parisi wherein nonphysical complications stemming from the replica trick were resolved through a phenomena called replica symmetry breaking. [35][36][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 include 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][31][32][35][36][37][38][39] In the past decade, the utility of the replica method has also been demonstrated in traditional structural glasses. [40][41][42][43][44][45][46][47][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 and the tracking of the glassy state as it approaches jamming upon compression. [40][41][42][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 an extremely 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-57 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. 58 RSA has some interesting properties, including a terminal (also called saturation) density [49][50][51][52]54 beyond which the process cannot be continued and an unusual logarithmic form 51,54 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. 59,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. 61 Other work has focused on the extension of replica integral equation theory [62][63][64] to RSA. 61,65,66 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][41][42][43][44][45][46][47][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. 65,66 The remainder of the manuscript is organized as follows. In Section 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 Sect. 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 Sect. 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][31][32] Originally developed for spin glasses, 33,34 it provides a recipe for extracting properties of the real quenched disorder 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. [62][63][64] 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. 67 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, [68][69][70] 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 are 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 g κ,n (r) = g κ,n (r|m → 0).
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. [62][63][64] 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. 71,72 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 which 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 broken and still generate a finite contribution. 19,[73][74][75][76] Beyond the fourth coefficient, the complexity grows rapidly; therefore, we truncate the expansion at fourth order. The final fourth order expansion, derived in Appendix C, after taking the where η i ≡ v D σ D ρ i and ρ i are the total volume fraction and number density, respectively, after the ith 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,73 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 Appendices D-E. 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,[77][78][79] 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,[77][78][79] The CS form is virtually exact over the entire fluid regime at all tested dimensions (and polydisperse mixtures) for equilibrium hard spheres. 43,80,81 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 Eqn. 8 to agree with that in Eqn. 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 Eqn. 8. Solving for equality of the quadratic density terms yields and As discussed in Sect. III, this approximate re-summed form (Eqn. 8) has a larger domain of validity than the low density expansion (Eqn. 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, 54 and the predictions do not improve with dimension. For convenience, tabulated values for Q 3 , Q 4 , Q 4 , Z 1 , Z 2 , and A from one to six dimensions are provided in Table I.
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 We will use both the simple expansion (Eqn. 3) and the CS improved form (Eqn. 8) to compare to exact simulation results.

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 Sect. II. To evaluate the accuracy of the above theory, direct simulations of the RSA process for systems in six different spatial dimensions were performed. 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.
A larger number O(100) of separate realizations were used to gather partial radial distribution functions in 3D. 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. 54 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 Sect. 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. 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). 82,83 The corrected (CS based) replica theory that approximately includes some of the missing higherorder 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 to 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 Sect. 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 Fig. 3a-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. 3c, we plot g(σ) from simulation between the first decile and the n th decile (going across the first row of the heat maps), and in Fig. 3d, 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 Fig. 3a-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.

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. 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,54 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][41][42][43][44][45][46][47][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. 65,66 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 Eqn. 7), but the resultant terminal density is far too high compared to 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.

Appendix A: Sequential replica trick
The RSA process for hard spheres can be interpreted as a sequential addition, equilibration, and positional freezing (quenching) protocol. 61 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 ith 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 κ (R 1:κ−1 ) ≡ Rκ exp[−U κ (R κ |R 1:κ−1 )] 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 F n ≡ · · · −lnZ n n−1 n−2 · · · 2 1 .
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 m is the partition function for a complex, equilibrium non-additive mixture of spheres. This is easy to see from the form of Eqn. 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) nonadditivity comes from the various powers of 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 Eqn. A6 with respect to u κ,n (|r 1 − r 2 |) 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 are 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. 71,72 To obtain an expansion in density we leverage the standard virial where B i is the i th virial coefficient. 19,71,73,74 Substituting Eqn. B2 into Eqn. B1 yields The composition dependence of the virial coefficients is apparent from the following decomposition into the species dependent VCs where n C is the number of components. 19,71,72 Eqn. B3 requires the derivative of Eqn. B4 with respect to σ a,b . Taking the derivative and collecting identical terms via the permutation symmetry of the species labels yields Substituting Eqn. B5 into Eqn. B3 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,73,74 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 functionf (r|σ) ≡ ∂f (r|σ)/∂σ which is related to the Dirac delta function, δ(r), viaf (r|σ) = −δ(r − σ). 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 onẽ f (r|σ) (dashed bond) where the integration is over a Cartesian coordinate associated with each node specifying a specific pair of species. 19,73,74 Specifically, the third order graph in (r a,b |σ a,b )f (r a,c |σ a,c )f (r b,c |σ 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 fourth order. Specifically, the fourth order term is more complex, 19,73,74 possessing multiple graphs with varying degrees of connectivity: As some of the graphs possess broken bonds, they can support "cross replica" contributions to the density expansion, as they are non-interacting.
Appendix C: Density expansion of the replica tree mixture contact value

Graphical description of replicated mixture interactions
Applied to the replicated mixture, the summation in Eqn. 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 (Eqns. 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 Γ indicates the position of node κ along the primary backbone relative to any "summed" nodes within the rectangular shaded "plate". The second virial term has no summed nodes so there is only one graph (hence the lack of the plate notation) whereas for the third virial graph κ can come before or after the summed α 3 node. Both graph sets in Eqn. C1 represent a primary backbone in the replica tree which is just one of the m 2 m 3 · · · m n continuously descending paths selected by a specific choice of replicas at levels as they have too many "broken" interactions between any one species. Specifically, α 4 and n have two species they do not interact with in the first and second graph types, respectively.
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.

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 Eqn. B6 according to the various graphically described contributions in Eqns. C1 and C2 yielding g κ,n (m) = 1 + Q 3 n α 3 =κ where i) Q i are defined by Eqns. 4-6; ii) we have changed from number density (ρ) to volume fraction (η) and recognized the replica species densities in Eqn. B6 correspond to incremental densities added in the RSA process (hence the ∆); iii) the factors of two account for permuting α 3 and α 4 ; iv) the factor of m α 4 − 1 comes from the dangling α 4 leaf in B Γ of Eqn. 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 η n ≡ n α=1 ∆η α further simplification yields g κ,n ≡ g κ,n (m → 0) = 1 + Q 3 (η n − η 1 ) The result in Eqn. C5 does not assume anything about how many additions are performed or what is the size of each increment. 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 η 1 = ∆η → 0, and summing the remaining sequence dependent terms in Eqn. 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 Eqn. 3.

Appendix D: Diameter derivatives of the fully-interacting composition dependent virial coefficients
The main text requires derivatives of the form ∂B a,b,α 3 ,...,α i /∂σ 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 σ dB a,b,α 3 ,...,α i dσ = This equivalence is a consequence of the definitions of partial and total derivatives and the functional form of B a,b,α 3 ,...,α i and B HS i -which are identical apart from explicit species labels and corresponding diameter labels. In the same limit, the right hand side of Eqn. D1 yields i(i − 1)/2 identical terms, which combined with Eqn. D2 yields 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 ∂f κ,γ (r κ,γ ) ∂σ κ,γ = −δ(r κ,γ − σ κ,γ ). 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. 74 Finally, we will utilize the definition of the vector norm, or equivalently the law of cosines, to write |r 1 − r 2 | 2 = r 2 1 + r 2 2 − r 1 r 2 cosθ (E7) where here θ is the angle between vectors r 1 and r 2 . Using these simplifications we arrive The integral in Eqn. E9 can be evaluated analytically for σ i,j = σ, which after dropping species labels yields 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 ρ 2 g(r) ≡ n i=1 n j=1 ∆ρ i ∆ρ j g i,j (r) We will find it useful to rewrite Eqn. F1 such that self and cross-terms are separate ρ 2 g(r) ≡ n i=1 ∆ρ 2 i g i,i (r) + 2 n i=1 n j=i+1 ∆ρ i ∆ρ j g i,j (r) (F2) In the continuous addition limit of RSA: ∆ρ i = ∆ρ = ρ/n and n → ∞, the first term in Eqn. 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.