Ekaterina B Zhulina

Oleg V Borisov

Bottlebrush polymer gels: architectural control over swelling and osmotic bulk modulus

Swelling behaviour and bulk moduli of polymer gels comprising of crosslinked bottlebrush subchains enables fine tuning by varying polymerization degrees of the main and side chains of the bottlebrush strands as well as their grafting densities. By using scaling approach we predic power law dependences of structural and elastic properties of swollen bottlebrush gels on the set of relevant architectural parameters and construct phase diagrams consisting of regions corresponding to different power law asymptotics for these dependences. In particular, our theory predict that bulk elastic modulus of the gel exhibits nonmonotonous dependence on the degree of polymerization of side chains of the bottlebrush strands.

Introduction

Architectural control of the performance polymer material properties by varying topology of the macromolecules without changing their chemical nature is a popular paradigm in contemporary material science.

Comb-shape polymers and molecular brushes comprising of long main chain with multiple side chains tethered to it at regular intervals present one of the most popular class of branched macromolecules used for fabricating more complex, including supramolecular architectures. [START_REF] Sheiko | Cylindrical molecular brushes: Synthesis, characterization and properties[END_REF][START_REF] Verduzco | Structure, function, selfassembly of bottlebrush copolymers[END_REF][START_REF] Müllner | Cylindrical polymer brushes -anisotropic building blocks, unimolecular templates and particulate nanocarriers[END_REF][START_REF] Rzayev | Molecular Bottlebrushes: New Opportunities in Nanomaterials Fabrication[END_REF][START_REF] Liang | Supersoft Polymer Networks with Brushlike Strands[END_REF][START_REF] Yuan | In Polymer Science: A Comprehensive Reference[END_REF] Macromolecules with even more complex topology based on the bottlebrush motive, e.g. barbwire macromolecules [START_REF] Liang | Synthesis and characterization of novel barbwire-like graft polymers poly(ethylene oxide)-g-poly(caprolactone)4 by the 'grafting from' strategy[END_REF][START_REF] Uhrig | Synthesis of Combs, Centipedes, and Barbwires: thinspace Poly(isoprene-graft-styrene) Regular Multigraft Copolymers with Trifunctional, Tetrafunctional, and Hexafunctional Branch Points[END_REF][START_REF] Pelras | Polymer Nanowires with Highly Precise Internal Morphology and Topography[END_REF] or dendronized polymers [START_REF] Rosen | Dendron-mediated self-assembly, disassembly, and selforganization of complex systems[END_REF][START_REF] Costanzo | Rheology and Packing of Dendronized Polymers[END_REF] have recently become available as well. Progress in chemical synthesis allowed for a wide variety polymers and copolymers with elaborately controlled branched architectures [START_REF] Zamurovic | Micellization Behavior of Complex Comblike Block Copolymer Architectures[END_REF][START_REF] Wang | Self-assembly of giant bottlebrush block copolymer surfactants from luminescent organic electronic materials[END_REF][START_REF] Li | Synthesis of Hetero-Grafted Amphiphilic Diblock Molecular Brushes and Their Self-Assembly in Aqueous Medium[END_REF][START_REF] Fenyves | Aqueous Self-Assembly of Giant Bottlebrush Block Copolymer Sufactants as Shape-Tunable Building Blocks[END_REF] that can, in turn, serve as precursors for more complex constructs, e.g., swelling networks with branched strands that could mimic elements of natural systems. [START_REF] Vashahi | Injectable hydrogels with tissue-adaptive gelation and mechanical properties[END_REF] Covalent cross-linking of comb-shaped polymers or molecular brushes gives rise to solvent-free elastomeric networks [START_REF] Self | Dynamic Bottlebrush Polymer Networks: Self-Healing in Super-Soft Materials[END_REF][START_REF] Xie | Bates Room temperature 3D printing of super-soft and solvent-free elastomers[END_REF] or polymer gels with novel unusual mechanical properties that can be tuned by proper variation in the architectural codon of the constituent strands (polymerization degrees of the main chains, i.e., cross-linking density, and side chains as well as grafting density). Such materials exhibit very unusual mechanical properties, and allow for independent tuning of elastic moduli in different regimes and in a wide deformation range. [START_REF] Daniel | Solvent-Free, Supersoft and Superelastic Bottlebrush Melts and Networks[END_REF] Remarkably, living nature amply employes bottlebrush structure for adjusting mechanical properties of different tissues [START_REF] Klein | Polymer in living systems:from biological lubrication to tissue engineering and biomedical devices[END_REF][START_REF] Kiani | Structure and function of aggrecan[END_REF][START_REF] Button | A periciliary brush promotes the lung health by separating the mucus layer from airway epithelia[END_REF] Recently we examined the swelling behavior of networks formed by barbwire molecular brushes with multiple short side chains emanating from each branching point in the loosely crosslinked network subchains. [START_REF] Zhulina | Polymer Networks Formed by Molecular Brushes: Scaling Theory[END_REF] The aim of the present work is to extend the scaling analysis for bottlebrush networks to the cases of arbitrary cross-linking density and degrees of polymerization (DPs) of the side chains, and to predict the power law dependences for the degree of gel swelling and its bulk elastic modulus as a function of relevant architectural parameters. chains each comprising n monomer units that are attached by one end to each m-th monomer unit of the main chain. Total number of monomer units in one strand of the network is N = M (1 + n/m). The main chain and side chains are assumed to be chemically identical, with Kuhn segment length on the order of monomer size a taken below as a unit of length. Below we address the case n m which corresponds specifically to the bottlebrush limit. Repulsive (in athermal solvent) interactions between densely grafted side chains lead to stretching of both side chains and spacers in the network strands and swelling of the gel. Two regimes of the bottlebrush gel swelling are distinguished below depending on the set of architectural parameters {M, m, n}.

The bottlebrush gel model

The semi-dilute C * -regime

In this regime the characteristic mesh size R (the average end-to-end distance of the main chain of a bottlebrush network strand) in the swollen gel exceeds by far thickness D of the individual bottlebrush, which is controlled by extension of the side chains, Figure 1a. The swollen gel can be treated as a gel formed by cross-linked linear chain strands with renormalized contour length and thickness (or, equivalently, effective excluded volume parameter). The equilibrium degree of gel swelling, Q = V /V dr R 3 /N , which is the ratio of volumes, V R 3 , and V dr N in the swollen and dry states of the network, respectively, can be found by applying the c * -theorem, [START_REF] De Gennes | Scaling Concepts in Polymer Physics[END_REF] according to which the equilibrium mesh size R of the gel is on the order of the size R eq of a single strand.

The individual bottlebrush strand can be assimilated [START_REF] Zhulina | Polymer Networks Formed by Molecular Brushes: Scaling Theory[END_REF][START_REF] Birshtein | Conformations of comb-like macromolecules[END_REF][START_REF] Feuz | Bending rigidity and induced persistence length of molecular bottle brushes: A selfconsistent-field theory[END_REF][START_REF] Theodorakis | Computer simulation of bottle-brush polymers with flexible backbone: Good solvent versus theta solvent conditions[END_REF][START_REF] Fredrickson | Surfactant-induced lyotropic behavior of flexible polymer solutions[END_REF][START_REF] Subbotin | Spatial Self-Organization of Comb Macromolecules[END_REF][START_REF] Zhulina | Solution and Melts of Barbwire Bottlebrushes: Hierarchical Structure and Scale-Dependent Elasticity[END_REF] to a selfavoiding wormlike chain with thickness D determined by extension of the side chains, and the effective contour length L ∼ = (M/m)h, where h is the average extension (end-to-end distance) of a spacer (segment of the main chain connecting neighboring branching points). Equivalently, the bottlebrush can be presented as a flexible chain composed of L/D impermeable superblobs each of size D. Repulsive (under good solvent conditions) interactions between densely grafted to the backbone side chains impose stretching of both side and main chains and impermeability of the superblobs. Hence, the characteristic size of a swollen network strand is

R R eq ∼ = L ν D 1-ν = ( M m ) ν h ν D 1-ν (1)
where [START_REF] Birshtein | Conformations of comb-like macromolecules[END_REF][START_REF] Zhulina | Solution and Melts of Barbwire Bottlebrushes: Hierarchical Structure and Scale-Dependent Elasticity[END_REF] h

∼ =    m ν (n/m) ν(1-ν)/2 , m ≥ n ν/(2+ν) , C * a m, m ≤ n ν/(2+ν) , C * b (2)
D ∼ = n ν (n/m) ν(1-ν)/2 , m ≥ n ν/(2+ν) n 2ν/(1+ν) m -(1-ν)/(1+ν) , m ≤ n ν/(2+ν) (3)
and ν ≈ 0.6 is the Flory exponent. The first and the second lines in eqs 2 and 3 correspond to partial and full (up to the contour length, m) stretching of spacers of the main chain, respectively. Below we refer, respectively, to subregime C * a (corresponding to m ≥ n ν/(2+ν) , partial spacer stretching) and subregime C * b , (corresponding to m ≤ n ν/(2+ν) , full spacer stretching). Scaling expressions, eqs 2 and 3, were derived in ref, [START_REF] Birshtein | Conformations of comb-like macromolecules[END_REF] based on the representation of a cylindrical brush of side chains with average axial distance h per chain as an array of densely packed concentration blobs of size ξ(ρ) ∼ = (ρh) 1/2 , where ρ is radial distance from the main chain, and volume fraction of monomer units of the side chains

φ(ρ) ∼ = ξ -(3ν-1)/ν (ρ) (4)
The condition of conservation of the number of monomer units in the side chain

h D 0 φ(ρ)ρdρ = n (5)
provides scaling relation for the brush thickness

D ∼ = n 2ν/(1+ν) h -(1-ν)/(1+ν) .
When the spacers of the main chain do not reach limit of their extensibility, i.e., h ≤ m, the equilibrium extension h of the spacers is found by balancing the elastic free energy of the extended spacer 24

F elastic k B T ∼ = h 1/(1-ν) m ν/(1-ν) (6)
with the free energy (in k B T units) evaluated as a number of concentration blobs per side chain

F conc k B T ∼ = h D O ρdρ ξ 3 (ρ) ∼ = (D/h) 1/2 (7)
Using eqs 1, 2 and 3 we can express the equilibrium dimensions of one network strand (the gel mesh size) as

R eq ∼ =    N ν (n/m) ν(1-3ν)/2 ∼ = M ν ( n m ) 3ν(1-ν)/2 , C * a N ν (n ν /m) (1-3ν)/(1+ν) ∼ = M ν n 2ν(1-ν)/(1+ν) m -(1-ν) 2 /(1+ν) , C * b (8)
As seen from eq 8, replacing of linear network strands by bottlebrushes with the same total degree of polymerization N leads to weaker swelling of the gel. The mesh size decreases as a function of n/m or n ν /m in regimes C * a and C * b , respectively. At the same time, attaching of side chains to the linear network strands with constant polymerization degree M leads to stronger swelling of the gel, i.e. larger R eq . Remarkably, as long as spacers of the main chains are only partially stretched (subregime C * a ), dimensions of the swollen gel depend only on the ratio n/m, whereas in the regime of full stretching of spacers (subregime C * b ) this universality breaks down. The osmotic modulus G of the swollen gel can be evaluated as 24

G k B T ∼ = R -3 eq ∼ =    N -3ν (n/m) 3ν(3ν-1)/2 ∼ = M -3ν (n/m) -9ν(1-ν)/2 , C * a N -3ν (n/m) 3(3ν-1)/(1+ν) ∼ = M -3ν n -6ν(1-ν)/(1+ν) m 3(1-ν) 2 /(1+ν) , C * b ( 9
) This result can be obtained also from the condition G osm ∼ = π where π is osmotic pressure which can be evaluated as

π/k B T ∼ = (φ * D ) 3ν/(3ν-1) D -3 (10)
and

φ * D ∼ = LD 2 R 3 eq ( 11
)
is volume fraction of D-superblobs in the swollen gel.

4 Concentrated regime (regime C * * )

The C * regime applies as long as the volume fraction φ * D of D-superblobs in the swollen gel is smaller than unity, i.e., D ≤ L. Alternatively, the boundary of the C * a and C * b regimes can be evaluated as equality of the average concentrations of the monomer units in the gel, c * , with

c * ∼ = N R 3 eq ∼ = (12)    N 1-3ν (n/m) 3ν(3ν-1)/2 = M 1-3ν (n/m) (3ν-2)(3ν-1)/2 , C * a N 1-3ν (n ν /m) 3(3ν-1)/(1+ν) = M 1-3ν n (2ν-1)(3ν-1)/(1+ν) m -(2-ν)(3ν-1)/(1+ν) , C * b ( 13
) and the concentration c * * of monomer units within distance D from the main chains of the bottlebrush strands,

c * * ∼ = n + m D 2 h ∼ = ∼ = n D 2 h ∼ =    n 1-3ν (n/m) ν(3ν-1)/2 , C * a (nm) (1-3ν)/(1+ν) , C * b ( 14
)
with h and D specified by eqs 2 and 3, respectively, and the last equality implying n >> m. Hence, C * a and C * b regimes apply as long as c * ≤ c * * , that is

M ≥    m(n/m) ν , m ≥ n ν/(2+ν) , C * a n 2ν/(1+ν) m -(1-ν)/(1+ν) , m ≤ n ν/(2+ν) , C * b ( 15
)
or

N ≥    n(n/m) ν , m ≥ n ν/(2+ν) , C * a n (3ν+1)/(1+ν) m -2/(1+ν) , m ≤ n ν/(2+ν) , C * b (16) 4.1 Subregime C * * a
For shorter main chain or longer side chains of the network strands the c *theorem is not applicable anymore. In the case when the main chains of the stands remain not fully stretched (subregime C * * a ) we can estimate the mesh size (i.e. the extent of the gel swelling) from the minimization of the free energy (per strand or per mesh volume)

F = F elastic + F vol ( 17
)
with the elastic free energy of the main chain 24

F elastic /k B T ∼ = R 1/(1-ν) M ν/(1-ν) (18)
and the excluded volume contribution

F vol /k B T ∼ = π k B T R 3 (19)
The osmotic pressure in the gel is given by π/k B T ∼ = c 3ν/(3ν-1) , and the concentration of monomer units in the gel is

c ∼ = N/R 3 (20)
The result of minimization of the free energy given by eqs 17, 18, 19 with respect to c gives equilibrium concentration of monomer units in the swollen gel

c ∼ = N 1-3ν (n/m) 3ν(3ν-1)/2 = M 1-3ν (n/m) (3ν-2)(3ν-1)/2 (21)
that coincides in terms of scaling dependences with eq 13 at m ≥ n ν/(2+ν) (regime C * a ). The dimensions R ∼ = (N/c) 1/3 of an individual network strand in this regime are given by the first line in eq 8 which grows as a function of n/m and reaches the limit of extensibility,

R ∼ = M , at n/m ∼ = M 2/3ν .
The gel osmotic modulus in this regime C * * a is given by

G k B T ∼ = π k B T ∼ = c 3ν/(3ν-1) ∼ = N -3ν (n/m) 9ν 2 /2 = M -3ν (n/m) 3ν(3ν-2)/2 (22)
which weakly decreases as a function of n/m at constant M .

c ∼ = N M 3 ∼ = M -2 n m ∼ = N -2 ( n m ) 3 (23)
and the bulk osmotic modulus of the gel in the regime C * * b is given by

G k B T ∼ = π k B T ∼ = c 3ν/(3ν-1) ∼ = N -6ν/(3ν-1) (n/m) 9ν/(3ν-1) = M -6ν/(3ν-1) (n/m) 3ν/(3ν-1) (24)

Diagram of states

In Figures 2a,b,c the diagrams of states of the bottlebrush gel are presented in M, n (Figure 2a), M, m (Figure 2b) and m, n (Figure 2c) coordinates, respectively. Regimes C * a and C * b correspond to loose networks, D R eq , where the main chains of the bottlebrush strands are not fully stretched (C * a ) or are fully stretched (C * b ), respectively (see Figure 1a). Regimes C * * a and C * * b correspond to dense networks where the main chains of the bottlebrush strands are not fully stretched (C * * a ) or fully stretched (C * * b ), respectively, (see Figure 1b). At the boundaries between the regimes smooth crossover of the gel mesh size, R R eq , occurs. As for the bulk osmotic modulus, the smooth crossover occurs at the boundary between regimes C * a and C * b , as given by eq 9, and between regimes C * * a and C * * b , as given by eqs 22 and 24. However, there is no smooth crossover for the osmotic modulus between regimes C *

(c -ν/(3ν-1) ) -3 . Therefore, G| C * * a /G| C * a ∼ = (n/m) 3ν/2 (25)
and

G| C * * b /G| C * b ∼ = (n ν /m) 3/(1+ν) (26)
In Figures 3a and3b we present the dependences of the gel mesh size R and the gel osmotic modulus G on the DP of side chains n for the case m ≥ M 1/3 . Similar dependences for the case m ≤ M 1/3 are presented in Figures 4a and4b, respectively. These dependences correspond to crosssections of the diagram in Figure 2c indicated by dash-dotted arrows.

(a) (a (a

n 2 m    m 3 2 3 m    2 3 n M m  * a C * b C ** a C ** b C 3 m c > 1 m 3 m M Comblike strands 2 n M m = 1 1 2 2 n M m       = 1 1 n M m     = = (b) M n n m * a C * b C ** a C ** b C 2 n    3 2 n    1 Comblike strands Unphysical values of parameters, M/m < 1 m M 1 1 1 m M n      = 2 3 m M n   = 1 2 1 1 m M n        = = (c) n M m 1 3 M * a C * b C ** a C ** b C M 2 3 M    2 M Comblike strands 1 n m = 2 n mM = 2 3 n mM  = 1 1 n M m     = 2 n m    = 1 1 2 2 n M m       =
) n R M 1/3 m M  m m 1 1 n M m     2 3 n mM  * a C ** a C ** b C 3 (1 )/2 n R M m           (b) n B G kT 3 M   1/3 m M  m m 1 1 n M m     2 3 n mM  * a C ** a C ** b C 9(1 )/2 3 M M m            3(3 2)/2 3 M M m           9 (1 )/2 3 n M m             3 (3 2)/2 3 n M m            3 /(3 1) 6 /(3 1) n M m              2 M 
) n R M 1/3 m M  2 n m    1 1 2 2 n M m       m m * a C * b C ** b C 3 (1 )/2 n R M m           2 2
(1 )/( 1)

(1 ) /(1 ) n R M m           (b) n B G kT 3 M   1/3 m M  m m * a C * b C ** b C 3 9(1 ) M m      2 3 6 (1 )/(1 ) 3(1 ) /(1 ) M n m             3 /(3 1) 6 /(3 1) n M m              2 n m    1 1 2 2 n M m       9 (1 )/2 3 n M m             3 M  3/2 ( ) Mm 

Discussion and conclusions

In the present paper we have studied conformational and elastic properties, that is, the swelling behavior and bulk osmotic modulus of the gels formed by crosslinked bottlebrushes. The scaling approach is suitable for polymer networks formed by densely grafted bottlebrush strands immersed in low molecular weight solvent with varying thermodynamic quality. This approach properly accounts for correlation in fluctuations in polymer density arising due to the gel hierarchical connectivity with explicit account of the network topology and goes beyond the mean-field models. Two major gel regimes were established: with non-overlapping (subregimes C * a , C * b ) and overlapping side chains belonging to different network strands within a mesh (subregimes C * * a , C * * b ). The results are summarized in the scaling-type diagrams of states (figures 2) with boundaries indicated in the plots. The power law dependences for the equilibrium gel parameters in different regimes are collected in Table 1.

Table 1. Equilibrium parameters of swollen gels with bottlebrush strands in various regimes of the diagram of states

R eq Q = V /V dr R 3 eq /N G π C * a M ν ( n m ) 3ν(1-ν) 2 M 3ν-1 ( n m ) -(3ν-2)(3ν-1) 2 M -3ν ( n m ) -9ν(1-ν) 2 C * * a M ν ( n m ) 3ν(1-ν) 2 M 3ν-1 ( n m ) 9ν(1-ν) 2 M -3ν ( n m ) 3ν(3ν-2) 2 C * b M ν n 2ν(1-ν) (1+ν) m -(1-ν) 2 (1+ν) M 3ν-1 ( n m ) (6ν-1) M -3ν n -6ν(1-ν) (1+ν) m 3(1-ν) 2 (1+ν) C * * b M M 2 ( n m ) -1 M -6ν/(3ν-1) ( n m ) 3ν (3ν-1)
The theory predicts that the gel mesh size, R, either increases upon increasing DP n of the side chains, or remains constant upon the approach to the limit of the backbone extensibility. A more complex behavior is predicted for the bulk osmotic modulus, G, which is expected to vary nonmonotonously (i.e., exhibit the local minimum) as a function of DP of the side chains and undergo abrupt changes at the boundaries between the regimes.

The power law dependences in Table 1 were derived based on the asymptotic behavior of bottlebrush strands with long side chains and backbones. That is, the side chains in the bottlebrush strands were presumed to be long enough to swell and stretch due to the intermolecular repulsions under good and theta solvent conditions. However in the majority of experimental systems the side chains and backbones of the strands are still relatively short that makes the behavior of these polymers closer to that of linear chains. The DPs of grafts in currently simulated bottlebrush networks [START_REF] Jacobs | Dobrynin Nonlinear Elasticity and Swelling of Comb and Bottlebrush Networks[END_REF] remain also ˜10 -30 monomers, which makes it difficult to compare the predictions of the presented scaling model (dealing with asymptotically long side chains and backbones) with the results of available experiments and simulations. Notably, a similar situation took place upon publication of the first scaling models of polymer brushes formed by linear chains, 3637 for which the theoretical predictions preceded the experimental realizations.

However, in some cases the comparison between our predictions and the experimental results is possible. As it follows from Table 1, in the subregimes C * a , C * b , and C * * a , at fixed strand architecture, n/m = const, swelling ratio V /V dr is predicted to vary with the strand backbone length M as V /V dr ∼ M 3ν-1 . Recent experiments on chemically crosslinked gels [START_REF] Sarapas | Compressing and Swelling To Study the Structure of Extremely Soft Bottlebrush Networks Prepared by ROMP[END_REF] implementing mononorbornene macromonomer for ROMP polymerization and dinorbornene cross-linker (both composed of poly(n-butyl acrylate)) demonstrated the increase in V /V dr upon increasing macromonomer to cross-linker ratio, [M M ]/[XL]. The latter correlates with the density of gel crosslinking [START_REF] Sarapas | Compressing and Swelling To Study the Structure of Extremely Soft Bottlebrush Networks Prepared by ROMP[END_REF] to approximate the average DP of the strand backbone as M ∼ [M M ]/[XL]. Each macromonomer contained approximately 30 n-butyl acrylate monomer units, and according to the scaling estimates, sizes of the peripheral concentration blob in the bottlebrush strand (modelled as cylindrical brush) and of the thermal blob are close. That is, the side chains in bottlebrush strands were (presumably) close to theta-solvent conditions, and we used the value of exponent v = 1/2 in comparing the experimental data to the theoretical estimate, V /V dr ∼ M 3ν-1 ∼ M 1/2 ∼ ([M M ]/[XL]) 1/2 . In Figure 5 the experimental swelling ratio Q = V /V dr in toluene at 37C o (Table 1 in ref [START_REF] Sarapas | Compressing and Swelling To Study the Structure of Extremely Soft Bottlebrush Networks Prepared by ROMP[END_REF] ) is presented as a function of [M M ]/[XL] in the logarithmic coordinates. The slope 1/2 indicated by the solid line clearly fits the experimental data better than the slope 4/5, corresponding to athermal solvent conditions with v = 3/5.

Another set of experimental data on physically cross-linked (self-assembled) networks formed by micro phase segregated triblock copolymer PNIPAM-PEO-PNIPAM with PEO bottlebrush block in toluene allows us to compare the theoretical and experimental swelling ratios, Q = V /V dr , more accurately. In contrast to a small number of strands emanating from each cross-link, and structural heterogeneity of chemically cross-linked gels, micro phase segregated ABA triblock copolymers give rise to networks with well defined domains (cross-links) and numerous bridges (strands) connecting neighboring domains. At low contents of condensing blocks A, the spherical domains surrounded by swollen bottlebrush blocks B can be assimilated to starlike micelles with aggregation number P 1, and the corona thickness H much larger than the core radius. Scaling model of starlike micelle with bottlebrush predicted by eqs (8), V /V dr ∼ R 3 eq /N ∼ M 3ν-1 , is indicated by solid line for ν = 1/2 (theta solvent conditions) and dashed line for ν = 3/5 (athermal sovent conditions). corona [START_REF] Zhulina | Micelles Formed by AB Copolymer with Bottlebrush blocks. Scaling Theory[END_REF] allows us to express the network swelling ratio V /V dr as

 20 0 4 Q 1 i  Q 2 i  Q 3 i  4 20 0 20 0 4 Q 0 i  Q 0 i   Q 0 i   4  20
Q theor = V V dr theor = φ L 4π 3
H 3 P n L + 1 = φ L 4π 3 k (M/2) 9/5 n 33/25 L + 1 (27)
Here, the first term in brackets is ratio of volumes of block B in the swollen corona of micelle and of the condensed block A, the corona thickness H ∼ (M/2)

φ L = n L n L + M (n/m + 1)/2
is volume fraction of A in ABA triblock copolymer, and monomer units A and B are assumed to have equal sizes.

In Figure 6 we plot the experimentally measured swelling ratio Q = V /V dr for series of samples in Table S8 in ref. [START_REF] Vashahi | Injectable hydrogels with tissue-adaptive gelation and mechanical properties[END_REF] versus Q theor calculated according to eq 27 with prefactor k = 1, and values of φ L , n L , and M/2 from Table S8. (Note that in Table S8 DP of the bottlebrush block, M , is denoted as n bb ). As it is seen in Figure 6, the experimental and theoretical values of Q are consistent (follow the straight line with slope 1 in log-log coordinates, and (as expected) the numerical difference between them is on the order of unity.

Overall, the correspondence between the theoretical predictions and the experimental data is reasonable taking into account the number of approximations involved. Our theoretical predictions could be therefore used for rationalizing the experimental data on samples with systematically varied parameters and might promote experimental and computational studies of structural and elastic properties of bottlebrush gels. The developed scaling model provides guidelines for rational control of structural and mechanical properties of polymer gels by tuning architecture of the subchains, and thus might contribute to the design concept of new materials with programmed properties. theoretical swelling ratio V/V dr experimental swelling ratio V/V dr Figure 6: Equilibrium swelling ratio V /V dr of physically crosslinked gel [START_REF] Vashahi | Injectable hydrogels with tissue-adaptive gelation and mechanical properties[END_REF] formed by microsegregated triblock copolymer PNIPAM-PEO-PNIPAM with PEO bottlebrush block in toluene at 37C o (data for samples in Table S8 in ref [START_REF] Vashahi | Injectable hydrogels with tissue-adaptive gelation and mechanical properties[END_REF] with varied block copolymer composition φ L = n L /[M (n + 1)/2 + n L ] with DP of the side chains n = 9, DP n L of linear PNIPAM block, DP M of the PEO bottlebrush backbone) versus theoretical dependence V /V dr = [(4π/3)H 3 /P n L + 1]φ L predicted in ref [START_REF] Zhulina | Micelles Formed by AB Copolymer with Bottlebrush blocks. Scaling Theory[END_REF] for spherical starlike micelles with aggregation number P ∼ n 4/5 L and corona thickness H ∼ M 3/5 P 1/5 (regime ST, all other prefactors are asigned unity). Slope 1 is indicated by solid line.

Figure 1 :

 1 Figure 1: Schematics of a swollen bottlebrush gel in semi-dilute, C * (a) and concentrated, C * * , (b) regimes. The number of monomer units in the main chain is M . The numbers of monomer units in a side chain is n, the number of monomer units in a segment of the main chain separating neightbouring side chains (spacer) is m.

4. 2

 2 Subregime C * * b In the regime C * * b the main chains of the bottlebrush strands are fully stretched, R ∼ = M , the concentration of monomer units in the gel is

  a and C * * a , as well as between regimes C * b and C * * b . Indeed, in regimes C * a and C * b the modulus at the boundary C * a /C * * a and C * b /C * * b is proportional to D -3 , whereas in regimes C * * a and C * * b at the same boundary it is proportional to the inverse volume of the concentration blob ∼

Figure 2 :

 2 Figure 2: Scaling diagrams of state of a swollen bottlebrush gel in M, n (a), M, m (b) and m, n (c) coordinates. Equation for boundaries between regimes are indicated. The region m ≥ M in panel (a) is unphysical. The regions n ≤ m in panels (b) and (c) correspond to comb-shaped network strands with loosely grafted side chains not considered here. The regions shaded in blue correspond to polymer volume fraction larger than unity, that is, are unphysical.

Figure 3 :

 3 Figure 3: Gel meshsize R (a) and the gel modulus (b) as a function of DP of side chains n for the case m ≥ M 1/3

Figure 4 :

 4 Figure 4: Gel meshsize R (a) and the gel modulus (b) as a function of DP of side chains n for the case m ≤ M 1/3

  3/5 P 1/5 with M as DP of the backbone in block B, the aggregation number P ∼ n L with n L being DP of block A (regime ST in ref 34 ). Prefactor k incorporates thickness D B of the bottlebrush, while

	4/5

Acknowledgements

This work was supported by Russian Science Foundation grant 20-13-00270


Frank Smallenburg

Efficient event-driven simulations of hard spheres

Efficient event-driven simulations of hard spheres

I. INTRODUCTION

Hard spheres are perhaps the most commonly studied model system in soft matter physics. In addition to being excellent model systems for certain colloidal suspensions [1][2][3], their simplicity makes them a fundamental model study for many-body systems at any scale. As a result, simulations of hard spheres have been instrumental in shedding light on many aspects of statistical physics, including phase behavior in bulk [4][5][6][7][8][9][10][11] and in confinement [12][13][14][15][16][17][18], crystal nucleation [19][20][21][22][23][24][25][26], and glassy materials [27][28][29][30][31][32][33][34][35][36].

Due to their popularity as a model system, it should come as no surprise that significant efforts have been made to design efficient hard-sphere simulations, typically based on either Monte Carlo (MC) methods or event-driven molecular dynamics (EDMD). Monte Carlo methods have been applied to hard-sphere systems since the 1950s [37,38] and continue to be in active use. For a hard-sphere system in the canonical ensemble, the basic MC algorithm consists of simple single-particle displacement moves, which are accepted as long as no overlaps are created [START_REF] Frenkel | Understanding Molecular Simulations: From Algorithms to Applications[END_REF]. This approach can readily be adapted to a variety of thermodynamic ensembles, or to incorporate biasing potentials or external fields, and scales well to large systems. Moreover, intelligently chosen MC trial moves, such as cluster moves [START_REF] Dress | [END_REF][41][42] or swap moves [43], can lead to significant speedups in sampling speed. In particular, event-chain algorithms, which move chains of particles collectively in a rejection-free move, are highly efficient for hard-sphere systems [44,45]. Depending on the goal of the simulations, one potential downside of MC simulations is the fact that, unless special care is taken [46,47], the dynamics do not typically correspond well to that of a real-world system.

Event-driven molecular dynamics (EDMD) simulations are a second widely used method for simulating hard-sphere systems, with a similarly long history [48]. These simulations solve for the motion of the particles using Newton's equations of motion, effectively model-ing the movement of hard spheres in a vacuum 1 . Following Newton's third law, hard spheres move in simple straight lines until they collide with a neighboring particle. EDMD simulations rely on predicting future collision events between particles, and then resolving each of these collision events in order. In comparison to standard MC simulations, hard-sphere EDMD simulations are extremely efficient, rivaling e.g. event-chain methods [45], with the added benefit of reproducing physically realistic dynamics. This is ideal for the study of glassy dynamics or spontaneous nucleation, where observing the dynamics of interest requires simulation over long time scales.

Here, we present a light-weight and fast hard-sphere EDMD code, which combines an efficient event-queueing system [51], neighbor lists [52], and "lazy" event evaluation [53]. We benchmark the simulation code on a range of monodisperse and bidisperse hard-sphere systems, and find significant performance gains when comparing the results with other commonly used simulation codes [54,[START_REF] Rapaport | The art of molecular dynamics simulation[END_REF].

In the remainder of this article, we discuss some of the more detailed aspects of the efficient implementation of EDMD simulations (Sec. II), show benchmarks of three variations of the presented simulation code (Sec. III), and discuss the results (Sec. IV).

II. IMPLEMENTATION

In this section, we outline the implementation of an efficient event-driven hard-sphere simulation algorithm. After briefly outlining the typical setup of a hard-sphere EDMD, we focus on three implementation choices: the event calendar, the neighbor selection, and the choice of how many events to schedule for a given particle.

For the simulations described in this article, we always consider periodic systems in an orthogonal simulation box.

A. Basic outline

The basic premise of a hard-sphere EDMD is straightforward. Given a configuration of N particles with positions {r N } and velocities {v N }, the algorithm determines which pair of particles will be the first to undergo a collision, and when this collision will take place. Once the event is determined, the simulation time t jumps forward to the collision time, particle positions and velocities are updated to reflect the impact of the collision, and the simulation proceeds to predict the next collision.

Since we are simulating purely hard particles moving in a vacuum, we know that in between collision events the particles experience no forces, and hence move in straight lines. As a result, for hard spheres, predicting the collision time for a pair of particles can be done analytically based on the initial distance vector r ij = r j (t) -r i (t) and the relative velocity v ij = v j (t) -v i (t) between the colliding particles i and j. Specifically, the time until the next collision t col is given by [START_REF] Rapaport | The art of molecular dynamics simulation[END_REF]

t col = -b -b 2 -v 2 ij (r 2 ij -σ 2 ij ) v 2 ij , (1)
where b = r ij • v ij and the minimum approach distance σ ij = R i + R j is the sum of the radii R i , R j of the two spheres. Note that a collision is only expected to occur when b < 0 (i.e. the particles are currently moving towards each other), and the discriminant

b 2 -v 2 ij (r 2 ij -σ 2 ij ) ≥ 0.
When a collision occurs, the two particles involved will instantaneously undergo a change in their velocities. Assuming a perfectly elastic collision, this change conserves the total momentum and kinetic energy of the particles, resulting in velocity changes:

δv i = - 2m j m i + m j br ij σ ij , δv j = 2m i m i + m j br ij σ 2 ij , (2)
where m i is the mass of particle i, rij is a unit vector in the direction of r ij , and the values of b and rij are taken at the moment of collision.

Using these equations for predicting and resolving collisions, a naive implementation of the algorithm could simply proceed as follows:

1. Check all particle pairs for possible collisions, and record the first collision time t first .

2. Set t = t first and update all particle positions.

3. Process the collision by updating the velocities of the two colliding particles 4. Go to step 1.

In practice, this is rather inefficient, and there are a number of standard optimizations that can be implemented to improve performance. First, after a collision between particles i and j, any other collisions previously predicted for particle i and j are no longer valid, and hence the future collisions for these particles need to be recalculated. However, all predicted collisions that only involve particles other than i and j are still valid, since the trajectories of those particles did not change. To make use of this, we keep track of an event calendar : a list of predicted future events. Whenever a collision occurs, we then remove any future events that involve the colliding particles from the event calendar, but leave the others in place, avoiding the need to re-predict them. We will discuss the implementation of the event calendar in Sec. II B.

A second useful optimization lies in the fact the positions of any other particles are irrelevant for handling a collision between two particles. Hence, there is no need to update particle positions for any collision they are not involved in. Instead, for each particle i we can simply keep track of the last time t i that the particle was updated, and its position and velocity at that time, and update its position only when needed. This avoids the vast majority of particle updates in the algorithm described above.

Finally, we can save additional time by realizing that particles can only collide with particles in their immediate environment. Hence, predicting collisions between particles that are separated by distances many times larger than the particle size is rather wasteful: most likely these events will be invalidated by intervening collisions long before they would take place. Efficiently dealing with particles with a finite interaction range is a common problem in computer simulations, and the standard approach for addressing this is to only consider a limited selection of nearby neighbors when calculating interactions or predicting collisions [START_REF] Frenkel | Understanding Molecular Simulations: From Algorithms to Applications[END_REF]. In EDMD simulations, the most commonly used method is the use of cell lists [54,[START_REF] Rapaport | The art of molecular dynamics simulation[END_REF], but neighbor list (also called Verlet lists) have also been implemented [52]. In both cases, we introduce a new event type, a neighborhood event, which triggers an update of the neighbor list or cell list position for a given particle. These methods will be discussed in more detail in Section II C.

As a result, the basic outline of a modern EDMD simulation is as follows.

1. Initialize the cell list or neighbor lists.

For each particle, predict the point in time when

either it moves to a different cell or its neighbor list needs an update, and add this event to the event calendar.

3. For all neighboring particle pairs, predict future collisions and add them to the event calendar.

4. Locate the first collision or neighborhood event in the event calendar, and move the simulation time t forward to the associated event time.

5. If the current event is a neighborhood event, do the necessary bookkeeping on the particle environment, and predict new neighborhood and collision events for the updated particle.

6. If the current event is a collision, update the two colliding particles to time t, and process the collisions. Then, predict new neighborhood and collision events for both particles.

7. Go to step 4.

In the following subsections, we look at some of the implementation details in more detail.

B. Event calendar

The implementation of the event calendar in EDMD simulations has been discussed extensively in literature [51,53,[START_REF] Rapaport | The art of molecular dynamics simulation[END_REF][START_REF] Marın | [END_REF][57]. We are storing a large number of events (scaling with the system size N ), each containing (at least) information about the event type, the time at which the event will occur, and the particles involved. The data structure we choose for this has to allow for 1. rapidly finding the event that is first chronologically (i.e. with the lowest event time), 2. efficient addition of new events, and 3. efficient deletion of invalidated events.

The first requirement strongly suggests that the events in the event calendar should be kept sorted based on their event time. Using an unsorted event calendar, finding the first event would require a search through all stored events (O(N ) operations), which would be expensive to perform for every collision. A sorted data structure avoids this issue by ensuring that the next event in line is always easy to find. In particular, event calendars in the form of a binary search tree are commonly used [START_REF] Rapaport | The art of molecular dynamics simulation[END_REF], but other variations of binary tree structures have been investigated as well [57]. Regardless of the exact implementation, in binary tree structures the total cost of operations for handling a collision (i.e. the combination of finding the first event, deleting the old events for the particles involved, and inserting the newly predicted collisions) scales as O(log N ) in system size. Since the number of collisions per simulation time unit in an EDMD trivially scales as O(N ), this results in a total computational cost per time unit of the simulation of O(N log N ).

A more efficient event calendar strategy was introduced in Ref. 51, which proposed combining a binary tree structure for events in the near future with an array of linear lists for the remaining events. In particular, this strategy divides the simulation time into intervals of length ∆t, and chooses to only store events within the current time interval in a sorted binary tree structure as described above. The remaining events are stored inside a set of "event buckets", with each bucket corresponding to an unsorted list of all events that occur within the same time interval t ∈ [k∆t, (k + 1)∆t) (with k an integer). At the end of the current time interval (i.e. when no more events remain in the binary tree), the events from the next event bucket are added to the binary tree, and the simulation continues. Since the events in each bucket are implemented as unsorted doubly linked lists, addition or removal of an event from the buckets can be done at a fixed cost of O(1). Moreover, the cost of updating the initial binary tree scales only as O(log n), with n the typical number of events scheduled in a time interval ∆t. By choosing ∆t ∝ 1/N , this number n can be chosen to be independent of N , resulting in a O(1) cost for event scheduling and an overall EDMD simulation performance of O(N ) operations per simulation time unit. This yields a significant performance improvement over implementations that rely on a pure binary tree, especially for larger systems (i.e. N > ∼ 10 3 ) [51]. Note that in this approach, care needs to be taken that a sufficiently large number of event buckets are available. The number of buckets should typically be on the order of t max /∆t, where t max is the maximum expected difference between the current time and a predicted event time. In practice, in an EDMD simulation this maximum time is not bounded, as two particles that happen to move slowly or nearly parallel can experience a collision arbitrarily far in the future. However, as described in Ref. 51, this issue can be circumvented by the use of an "overflow" bucket which contains the (rare) events that are too far into the future. In all our simulation codes presented here, we follow the algorithm described in Ref. 51 for the event calendar, with one minor modification: we use a binary search tree rather than a complete binary tree [57] as our initial event tree. Since the scheduling of events into the (small) binary tree takes up only a small fraction of the overall simulation time, we do not expect this choice to make a significant difference.

C. Neighbor selection

The majority of the computational cost of an EDMD simulation is spent on the prediction and scheduling of collisions. Hence, reducing the number of collision partners to check for a given particle can drastically affect simulation times. Broadly, there are two main strategies for determining which particles to check.

Most commonly used are cell lists: the simulation box is divided up into a grid of cells that are sufficiently large to ensure that collisions can only occur between particles in adjacent cells (i.e. particles that are at most one cell away in along any of the three spatial dimensions). We can then limit ourselves to only predicting collisions with particles that are either in the same cell as the particle being checked, or in any of the 26 (in 3D) adjacent cells. Since the number of particles per cell is limited, finding the future collisions for a given particles then only requires checking a limited number of neighbors, rather than all N particles in the system. This efficiency comes at the (small) cost of having to keep track of the cell each particle is in. To this end, an extra event type is added that occurs when a particle crosses from one cell to the next. Like collisions, these cell-crossing events are predicted and added to the event calendar each time a particle is involved in an event. When a cell crossing occurs, particles in the newly adjacent cells (i.e. 9 cells in 3 dimensions) are checked for collisions.

Another common approach for neighbor selection is via the implementation of neighbor lists [START_REF] Frenkel | Understanding Molecular Simulations: From Algorithms to Applications[END_REF]. An efficient implementation of this approach for event-driven simulations was shown in Ref. 52. In this strategy, each particle maintains a list of particles that are considered to be "close enough" for possible interactions. In particular, we store for each particle i the center of its current neighborhood r neigh,i , and the set of other particles j such that

|r neigh,i -r neigh,j | < (1 + α)σ ij , (3)
with α > 1 a parameter that controls the size of the neighborhood of a particle. We can then conclude that particle i cannot collide with any particle k not in its neighborhood, unless either i has moved by at least αR i or k has moved by at least αR k from their respective neighborhood centers. Hence, we can safely limit ourselves to predicting collisions only between particles that are considered to be neighbors, as long as for each particle, we update its neighbor list when it reaches a distance of α times its radius from the center of its environment.

In addition to particle collisions, we then also include an event for each particle that takes place when

|r i -r neigh,i | = αR i . (4)
The prediction of this event is analogous to a normal sphere-sphere collision, with an event time

t neigh = -b n,i + b 2 n,i -v 2 i (|r i -r neigh,i | 2 -α 2 R 2 i ) v 2 i , (5)
with b n,i = r i • r neigh,i . The neighbor list update affects not only the list of neighbors of particle i, but will also cause i to be added to (or removed from) the neighbor lists of nearby particles. Note that a cell list is still used to find possible neighbors of a given particle, in order to limit the number of distances that need to be checked. The cell size in this case is chosen to be large enough to ensure that only particles in neighboring cells can be in each others neighbor list. Updates to the cell list are done as part of the neighbor list update.

Analogous to Ref. [52], we implement the neighbor list update of particle i as follows:

1. Update particle i to the current simulation time.

2. If particle i has left the simulation box, apply periodic boundaries.

3. Set r neigh,i = r i 4. Update the cell list with the new particle location.

5. For all particles j that were previously neighbors of i, remove i from the neighbor list of j and vice versa.

6. Check all 27 neighboring cells for particles j that satisfy: |r neigh,j -r neigh,i | < (1 +α)σ ij . For each of these, add i to the neighbor list of j and vice versa.

7. Delete any future events for particle i from the event calendar.

8. Predict and schedule new collision and neighbor list events for particle i.

We implement the neighbor list as a fixed-length array for each particle, allowing a maximum number n max neigh of neighbors that is constant during the simulation and should be chosen large enough to accommodate the maximum number of neighbors that a particle can have at the chosen shell size α. For the hard-sphere systems studied here, we found that a shell size α = 0.5 gave near optimal results (see SI). For this choice, n max neigh = 24 was sufficient in all cases 2 .

It should be emphasized that the cell list and neighbor list approaches will perform optimally under different circumstances. In particular, neighbor lists have excellent performance in systems where particles are likely to remain in the same environment for a long time, as is the case in solid phases and glassy fluids. Under these circumstances, the neighbor lists require few updates, and the number of calculated collisions is typically lower than in a cell list approach. For example, in a monodisperse hard-sphere crystal at a number density ρσ 3 ≈ 1, the neighbor list for each particle is likely to only contain its 12 nearest neighbors, while checking 27 cells of a size of at least σ 3 each requires checking on the order of 27 particles. In contrast, in low-density fluids, the high particle mobility leads to much more frequent neighbor list updates. As neighbor lists updates are more expensive than cell list updates, this favors the cell list method. Since long simulation times are more likely to be a concern in systems with relatively slow dynamics (e.g. glassy systems), we focus here on the neighbor list approach, but include benchmarks for one code (labeled CellList) using cell lists.

D. Lazy event scheduling

A key point of consideration in the implementation of an EDMD simulation code is the question of how many events to schedule per particle. In the above, we have implied that all predicted collisions (and the neighbor list update) for each particle are added to the event tree, but this is not necessarily the most efficient choice. Lubachevsky [53] noted that when predicting many events per particle, the later events are likely to be pre-empted by the earlier ones, and hence the computational cost of adding them to the event calendar might well be a waste. As a result, many studies (e.g. Refs. 52, 53, and 58) choose to schedule only a single event per particle. If the event scheduled for particle i is invalidated by a collision between two other particles, the event is typically converted into an update event, which simply rechecks for future collisions and neighborhood events without changing the particle trajectory. Although this means that sometimes the same collision is predicted more than once, the reduction in computation cost due to scheduling fewer events can result in significant speedups [58]. An alternative strategy, used in e.g. Refs. [54,[START_REF] Marın | [END_REF] does maintain all predicted events per particle, but still schedules only one event per particle in the event tree (limiting its size), while separately maintaining lists of sorted events per particle.

Here, we explore both the option of scheduling multiple events per particle and scheduling only a single event per particle. In particular, we consider two EDMD codes, both using neighbor lists and the event scheduling approach from Ref. 51, with the first (labeled Multi-Event) scheduling multiple events per particle, and the second (labeled SingleEvent) scheduling only one event per particle.

Specifically, in the MultiEvent code, we schedule all collision events that take place before the next neighbor list update event for both colliding particles. Since our implementation of the neighbor list update of particle i completely reevaluates the future events of i, there is no benefit to scheduling any events beyond this time.

In contrast, in the SingleEvent code, we only schedule the first (collision or neighbor list update) event for each particle. In order to check for invalidated events, each particle keeps track of the total number of collisions n col it has experienced. In particular, when updating the event for particle i (i.e. after a collision or neighbor list update), we do the following:

1. Predict the new neighbor list update for particle i, and set t min equal to the predicted time.

2. For all neighbors of i, predict a possible future collision. If a predicted collision time is smaller than t min , set t min equal to the collision time, and remember the collision partner.

3. Schedule the new (neighbor list or collision) event for particle i at t min .

4. If the predicted event was a collision (with particle j), store n col,j as additional information with the event.

Then, when handling the collision that was predicted for particle i with collision partner j, we:

1. Remove the event from the event calendar.

2. Check that the value of n col,j stored with the event predicted for particle i is equal to the current value of n col,j .

3. If not, then the collision is no longer valid. Predict the next event for particle i, and move on to the next event (ignore the steps below).

4. Otherwise, resolve the collision between the two particles, and increment n col,i and n col,j .

5. Remove the scheduled event for j from the event calendar, and predict new events for both particles.

Our approach is similar to the one followed in Ref. 52, with the distinction that we make no attempt to ensure that collision predictions are symmetric between collision partners. For example, if the first collision found for particle i is with particle j, then this does not affect the predicted collision for particle j, regardless of the predicted collision time. Since the i -j collision will be handled when the event associated with particle i is handled, there is no need to also schedule the same event in reverse. Moreover, the original prediction for j may still take place if the i -j collision is invalidated, so it is beneficial to leave the original event prediction for j in place.

In addition to saving significant amounts of scheduling effort and memory usage, this scheduling approach allows for a few small optimizations in the collision prediction. First, the strict time limit on predicted events can sometimes avoid the full calculation of the collision time in Eq. 1, by noting that t col is always greater than (r 2 ij -σ 2 ij )/2b. Hence, if this lower bound is already greater than the current t min , there is no need to calculate the exact collision time. Additionally, the total number of events that are scheduled is now equal to the number of particles, plus a small number of special events that are used for e.g. performing measurements, writing simulation output, or applying a thermostat. This eliminates the need for separate data structures for events and particles, allowing us to combine both into a single structure. This may help in optimizing memory access. Finally, the typically small times for predicted events means that the number of event buckets used in the event calendar can be quite limited, reducing memory usage.

III. BENCHMARKS

We test three versions of our simulation code (Cell-List, MultiEvent, and SingleEvent) on both monodisperse and bidisperse hard-sphere fluids, and compare their performance to other available EDMD codes [54, 55]. Some additional practical details of our simulation codes are described in the Supplemental Information (SI). Initial disordered configurations were created using a variation of the SingleEvent code where all particles grow in size during the simulation. This results in a rapid compression from a dilute random configuration to the desired packing fraction. After reaching the desired packing fraction, the system was then equilibrated using the SingleEvent code. In order to cut down on the computational cost of equilibration (which can be slow in the dense binary system), all initial configurations for systems larger than N = 10 5 particles were created by duplicating a smaller configuration either 2 × 2 × 2 times or 3 × 3 × 3 times.

For all benchmarks, the simulation is run for at least 10 7 collisions, and we measure the average number of collisions processed per second, starting the measurement after the initialization of the system. All three simulation codes are implemented in C. Compilation and optimization details are given in the SI. All benchmarks were performed on a workstation containing two 3.3GHz Intel Xeon Gold 6234 CPUs, with 8 CPU cores each, using Ubuntu 20.04 as the operating system . For each benchmark, both CPUs were filled by running 16 identical jobs in parallel. Hyperthreading was disabled.

A. Monodisperse hard spheres

First, we simulate a monodisperse hard-sphere fluid at packing fraction η = 0.49, for a range of system sizes N ranging from 500 to 1.35 × 10 6 particles. This packing fraction is slightly below the freezing packing fraction [START_REF] Frenkel | Understanding Molecular Simulations: From Algorithms to Applications[END_REF]. We compare all three codes to both the code by Rapaport, provided with Ref. [START_REF] Rapaport | The art of molecular dynamics simulation[END_REF] 3 , and the simulation package Dynamo developed by Bannerman et al. [54], version 1.7.6.

In Fig. 1, we show the number of collisions per second processed by each simulation. We observe a significant variation in simulation speeds, which is largely consistent as a function of system size. It is interesting to note that even for this relatively mobile hard-sphere fluid, the introduction of a neighbor list (i.e. switching from the CellList to MultiEvent code) results in a speedup of around 50% accross all system sizes. Combining this with the choice to store only one event per particle (i.e. the SingleEvent code) provides another similar speedup. In total, the SingleEvent code outperforms Dynamo by a factor 4-10, depending on system size.

All benchmarked simulations show a noticeable decrease in efficiency around a system size of ≈ 10 4 , which we attribute to cache misses. The simulation methods examined here typically use a few hundred bytes of memory per particle, and the total cache size of the used processors is approximately 2.5MB per core. This means that for small systems (up to a few thousand particles), all data relevant to the simulation fits inside the processor's cache. For larger systems significant data exchange between CPU and the main memory is needed, which drastically slows down performance. As some of the CPU cache is shared between cores, this performance drop can be shifted to slightly larger system sizes by running only a single job at a time (see SI). Note that in EDMD simulations, particle data is typically accessed in an essentially random order: two consecutive collisions can take place anywhere in the system, and particles that are close together in space may not be stored close together in memory. As a result, optimizing EDMD simulations with respect to cache usage is inherently difficult.

As Fig. 1 shows, the Rapaport code shows a stronger decrease in performance for large systems than the other methods. This is due to its use of a O(log N ) binary search tree as the event calendar, in contrast to the O(1) scheduling methods used by the other codes.

B. Bidisperse hard spheres

As our second test system, we take a binary hardsphere system with a size ratio q = 0.85 and a composition x L = N L /N = 0.3 (where N L is the number of large spheres in the system). This system has been explored in several studies of glassy dynamics [34][35][36]. Note that since the Rapaport code [START_REF] Rapaport | The art of molecular dynamics simulation[END_REF] implements only monodisperse systems, we did not include it in the comparison here. For this binary mixture, we first benchmark the simulation performance as a function of system size at a fixed packing fraction η = 0.58. The results are shown in Fig. 2. We observe similar results as for the monodisperse system, with again significant performance gaps between the different simulation codes. In this system, dynamics are glassy: particles are trapped by their neighbors for significant amounts of time before escaping their cage. As a result, neighbor updates are relatively rare, and hence the effect of incorporating a neighbor list is more significant than in Fig. 1: it results in a speed-up of approximately a factor two. Additionally, the typical number of events predicted per particle in the MultiEvent code is higher in this system than in the lower-density monodisperse system, resulting in an increased performance boost from switching to the SingleEvent code.

In order to investigate the effect of packing fraction on performance, we also benchmark each simulation code at a fixed system size of N = 104 particles for a range of different packing fractions. As shown in Fig. 3, we find a significant variation in the efficiency of the different codes as the packing fraction changes. As already hinted at in Sec. II C, the CellList code performs excellently at low packing fractions, but slows down as the packing fraction increases. This is understandable: as the packing fraction increases, more and more particles have to be checked for collisions. As most of these events will be pre-empted by an earlier collision, many of these calculations are wasted, resulting in an overall slowdown with increasing packing fraction. The step-wise performance decreases that can be seen for both the CellList code and Dynamo in Fig. 3 are due to the discrete reduction the number of cells that fit in the simulation box. A step occurs at a packing fraction where the number of cells is reduced by 1, which results in a discontinuous increase in the average number of particles per cell, and hence an increase in the computation cost of checking for collisions. Similar steps can be seen even for the MultiEvent and SingleEvent codes, but only at lower packing fractions (η < ∼ 0.5) where neighbor list updates (which rely on the cell list) account for a significant part of the simulation time.

In contrast to the CellList code, in terms of collisions per second, the neighbor-list based codes become faster with increasing packing fraction. As the packing fraction increases, particle mobility decreases and hence fewer neighbor list updates are required 4 . As a result, both neighbor-list driven codes outperform the CellList code at packing fractions η > ∼ 0.4. Moreover, at higher packing fractions the number of unnecessary events predicted scheduled per particle in the MultiEvent code is larger, and hence the SingleEvent code becomes the fastest option by a significant margin. Note that most hard-sphere simulation studies where long simulations are required are likely to focus on the high-packing fraction fluid regime, where e.g. glassy behavior or spontaneous crystal nucleation can be observed.

IV. DISCUSSION AND CONCLUSIONS

The simulation codes presented here (and in particular the SingleEvent code) clearly perform well in comparison to other publicly available EDMD implementations. It should be noted that part of this efficiency is due to the fact that the code was written with only hardsphere systems in mind. In contrast, Dynamo [54] offers the possibility to simulate different particle shapes, different boundary conditions, stepped pair interactions, as well as a variety of other features. However, for the common use case of studying hard-sphere systems in a periodic rectangular box, as is often necessary for e.g. glass or nucleation studies, the presented difference in simulation speed (typically a factor 5 to 10) may be a crucial factor, saving significantly on both CPU time and energy consumption.

It should be noted that if the goal is simply to equilibrate a hard-sphere system or to sample equilibrium states, it is not a priori clear that EDMD simulations are the ideal choice. Recent work by Klement et al. demonstrated that hard-sphere systems can be highly efficiently equilibrated using Newtonian Event Chains (NEC) [45], an approach that blends aspects of Monte Carlo simulations and event-driven dynamics. Simulating a monodisperse system of hard spheres at η = 0.49, they observed that their NEC simulations were significantly faster than EDMD simulations in terms of particle displacements per hour, which directly translated into a similar speedup in terms of the simulation time required to simulate e.g. particle diffusion or crystal melting. In Table I, we compare the number of displacements per second achieved by the currently available implementation of NEC [59] in HOOMD [60] to the number of collisions per second of the SingleEvent code for several system sizes. We see that a sufficiently efficient EDMD code can indeed slightly outperform this implementation of NEC 5 . However, it should be noted that since the complexity of an NEC displacement is lower than that of an EDMD collision, further optimization of the NEC implementation may be possible. Moreover, as also noted in Ref. 45, the local nature of the NEC algorithm allows for more straightforward parallelization [61]. On the other hand, the EDMD approach has the significant advantage of producing physically realistic dynamics.

Although here we focus purely on single-processor simulation codes, simulations of large systems can often be parallelized to decrease simulation time, although this typically carries an overhead cost in total CPU time.

For Monte Carlo methods, highly parallelized approaches have been shown to be capable of extremely fast equilibration of 2D hard-disk systems [62,63]. A similar approach may be effective for 3D hard spheres as well.

In contrast, EDMD simulations are notoriously difficult to parallelize [64], due to the fact that two events that are spatially far apart could influence each other even if they occur at nearly the same time, as long as they are connected by a chain of (nearly) touching particles. Nonetheless, several effective methods for parallelization of EDMD simulations have been demonstrated (see e.g. [58,64]), and can in principle be extended to the EDMD codes presented here.

As a final note, it is interesting that in all of the benchmarks, we found that the SingleEvent code performed at least as efficiently as the MultiEvent code. This observation is likely to be general for hard spheres and any extensions of the code to other spherically symmetric systems. In these systems, the cost of predicting a collision is low, and hence the cost of having to sometimes recalculate the same prediction is acceptable. However, it is worth noting that in EDMD simulations of anisotropic particles, such as polyhedra, ellipsoids, or patchy particles [52,[65][66][67][68][69], the computational cost of predicting a collision event is typically much higher, since collisions for these particles cannot be predicted analytically. This may favor an approach where multiple events are stored per particle, to avoid repeated predictions of the same event. On the other hand, the restriction to a single event per particle also helps to limit the time window to search for possible collisions, hence providing an efficient "early out" for numerical collision-prediction routines. Hence, the benefit of single-event scheduling approaches may depend on the system under consideration.

In summary, we have presented a set of light-weight EDMD simulation codes for hard-sphere systems that appear to perform very efficiently in comparison to the the simulations currently used in literature on e.g. the glassy dynamics of hard spheres. Simulating these systems more efficiently allows for a significant saving of CPU time and energy.

The simulation codes are freely available at: https://github.com/FSmallenburg/EDMD

FIG. 1 .

 1 FIG. 1. Benchmark results on monodisperse hard-sphere fluids for different simulation codes and a range of system sizes N . The label Rapaport indicates the code provided with Ref. [55]. The label Dynamo indicates the simulation package Dynamo [54].

FIG. 2 .

 2 FIG.2. Benchmark results on bidisperse hard-sphere fluids with size ratio q = 0.85, composition xL = 0.3, and packing fraction η = 0.58, as a function of the system size N .

FIG. 3 .

 3 FIG.3. Benchmark results on bidisperse hard-sphere fluids with size ratio q = 0.85, composition xL = 0.3, and system size N = 10 4 , as a function of the packing fraction η.

5

  It should be noted that the number of NEC displacements per second in Table I is significantly lower than those reported in Ref.45, especially for larger systems. We attribute this to our protocol of filling up the CPU with an identical number of jobs when performing benchmarks. When running only a single job per CPU, we observe performance in line with Ref.45, with EDMD still outperforming NEC in terms of displacements per second.

	N NEC displacements/s EDMD collisions/s
	1 000	1.0 × 10 6	1.6 × 10 6
	10 000	8.8 × 10 5	1.2 × 10 6
	100 000	4.1 × 10 5	4.3 × 10 5
	1 350 000	2.2 × 10 5	3.7 × 10 5

TABLE I .

 I Comparison between Newtonian Event Chain (NEC) simulations[45, 59], and the SingleEvent EDMD code, for a selection of the systems shown in Fig.2. Following Ref.45, we compare the number NEC particle displacements per second to EDMD collisions per second.

It should be noted that variations on EDMD that simulate Brownian[49] or Langevin[50] dynamics have been developed as well.

Note that in principle, it is possible to implement a variation on this approach where the shell size of a specific particle is temporarily reduced if too many neighbors are encountered[52], but this was found to be unnecessary for the systems considered here, where a good upper bound on the maximum number of nearest neighbors can be estimated.

Note that small modifications to the initialization and output of the code were made to allow the code to read in an initial configuration and report the elapsed time during the simulation.

Note that this effect does not hold for the cell list updates in the CellList code: for a homogeneous system, the average (bidirectional) flux of particles through any cell wall is only a function of the particle density and the typical speed (i.e. temperature) of the particles. While some particles may be trapped for a long time in a single cell, this is compensated by other particles trapped vibrating around the edge of a cell and hence undergoing many cell crossings.

V. ACKNOWLEDGEMENTS I would like to thank Michael Engel, Marco Klement, and Joshua Anderson for helpful discussions and aid in the comparison to the Newtonian Event Chain algorithm, and Laura Filion for many useful discussions.