Igor Shvab

Laurent Brochard

Hegoi Manzano

Enrico Masoero

Precipitation mechanisms of mesoporous nanoparticle aggregates: off-lattice, coarse-grained, kinetic simulations

Controlling the evolution of morphology and properties during the precipitation of mesoporous materials is a key technological challenge. This entails the scientific challenge of modelling and simulating a complex mesoscale kinetics. We present an original off-lattice Kinetic Monte Carlo approach to simulate various precipitation mechanisms at the mesoscale of nanoparticle aggregates. The simulations are based on novel coarse-grained rate expressions of nanopartice precipitation/dissolution, accounting for both solution chemistry and mechanical interactions. The precipitation of ordered and amorphous domains is simulated, showing how particle-particle and particle-substrate interactions determine various mechanisms: layer-by-layer precipitation, islands formation, Cahn and Avrami nucleation and growth, and gel-like precipitation. The simulations clarify how the total precipitation rate depends on the triggered mechanism, and therefore on solution chemistry and on mechanical interactions. This brings together chemical kinetics and nanoparticle simulations for a more fundamental understanding of mesostructure development, towards a computer-aided design of mesoporous materials.

I. INTRODUCTION

Mesoporous materials, with pore size between 2 and 50 nm, are very important for our economy. 1,2 Their production often involves nanoparticle precipitation from supersaturated solutions, during which some key structural features develop at the mesoscale between the nanometre and the micrometre. [START_REF]Nanostructured Soft Matter: Experiment, Theory, Simulation and Perspectives Zvelindovsky[END_REF] Being able to tailor this mesostructure would help improving the strength, functionality, and durability of important engineering materials, such as mesoporous metal oxides, zeolites, and cement paste. [START_REF] Van Vliet | [END_REF]5 The current approaches to mesostructural optimization are mostly empirical [6][7][8] . What is missing is a theoretical understanding of the relationship between mechanisms of nanoparticle precipitation/aggregation and experimental synthesis variables, such as solution chemistry and temperature. 9,[START_REF]Mesoscale Chemistry: A Workshop Report[END_REF] Several challenges limit the current models and simulations of nanoparticle precipitation and aggregation, hindering a more fundamental understanding of mesostructure formation. [START_REF]Mesoscale Chemistry: A Workshop Report[END_REF] Here we will address two of them:

• The timescale of nanoparticle precipitation is usually too large for molecular dynamics simulations [START_REF] Yip | [END_REF] . Kinetic Monte Carlo overcomes this limitation, but it requires a priori knowledge of the rates of all possible events [START_REF] Voter | Introduction to the Kinetic Monte Carlo Method[END_REF] , and typically there is an infinity of possible events (e.g. infinite possible positions for the nucleation of a new nanoparticle).

• The rates of nanoparticle precipitation or dissolution depend on free energy changes from both chemical reactions and mechanical interactions. Typically, the models of chemical kinetics, e.g. Classical Nucleation and Growth, neglect or strongly simplify the effect of mechanical interactions between particles and the role of entropy-induced fluctuations. [START_REF]Mesoscale Chemistry: A Workshop Report[END_REF][START_REF] Avrami | [END_REF][14][START_REF] Lasaga | Kinetic Theory in the Earth Sciences[END_REF] Simulations based on the Phase Field Method [START_REF] Boettinger | [END_REF] have started to address the coupling between mechanical stress and chemistry, but they become extremely onerous for systems with thousands of particles with local interactions. [START_REF] Boettinger | [END_REF]17 Vice versa, particle-based simulations consider aggregation and entropy-induced fluctuations in detail, but they usually simplify the chemical evolution of the system using effective rates of precipitation that are not rigorously linked to the underlying chemical reactions. 18,19 To address these two challenges, we develop an original, off-lattice Kinetic Monte Carlo (KMC) simulation approach for nanoparticle precipitation and aggregation. The rates of nanoparticle formation and dissolution are coarse-grained from the molecular scale and account for both chemical reactions and mechanical interactions. The simulations are first benchmarked against simple scenarios with analytical solution. The complexity is then increased to consider the precipitation of idealised, ordered mesoporous materials. Finally we simulate the precipitation of amorphous mesoporous domain, referring to cement hydrates to illustrate how the rates can be parametrised for real materials. The simulations predict various mechanisms of mesostructure formation and show how chemical drive and mechanical interactions can affect the overall precipitation rate.

II. METHODOLOGY

We propose two original methodological developments: (i) a Kinetic Monte Carlo approach that reduces the infinite number of possible positions for nanoparticle formation to a finite off-lattice sampling; (ii) new coarse- (b) Same configuration with trial particles (yellow) that can move locally within their Vc and interact with existing particles but not with each other. (c) Molecular process behind nanoparticle formation: n molecular units of solid form radially via chemical reactions that change the chemo-mechanical equilibrium free energy of the system ∆Geq and have an associated activation energy made of two contributions: ∆Ga,0, which depends on the chemical reaction and on the solution chemistry, and a fraction χ of ∆Geq (more details in the text).

grained expressions for the chemo-mechanical rates of nanoparticle insertion and deletion.

The formation of a nanoparticle is considered as the product of multiple chemical reactions in an implicit background medium, e.g. aqueous solution of salts (Figure 1). For simplicity, we will move within the boundaries of the four assumptions listed below, but the KMC sampling and coarse-grained rates that we will develop will be a starting point to treat cases also beyond these assumptions.

• The chemical composition of the medium is uniform in space but can change with time. Hereafter we will always consider the medium as a liquid solution from which solid nanoparticles can precipitate.

• The nanoparticles to be inserted or deleted are spherical, isotropic in terms of their internal molecular structure, monodisperse in size, and have all the same chemical composition.

• The particle size, shape, and chemical composition do not change during the simulations. Therefore only particle insertion and deletion will be considered; no partial growth, dissolution, or leaching.

• The position of the existing particles is changed only by energy-minimisation between subsequent KMC steps of particle insertion/deletion. This means that we consider particle displacements as instantaneous compared to the slow kinetics of solidification and dissolution.

• The mechanism behind nanoparticle formation is classical growth starting from preexisting nuclei. We assume that the background medium contains a certain concentration C cn of already-formed nuclei of solid consisting of one molecular unit (red in Figure 1c) and that these grow via sequential attachment of n other molecular units of solid in all radial directions. This corresponds to assuming that the critical nucleus of solid is as small as a single molecular unit, and that the nucleation time is much smaller than the time for the nucleus to grow until reaching the size of the considered nanoparticle. An explicit account of nucleation and other possible non-classical mechanisms is left to future work 20,21 (a discussion on including classical nucleation can already be found in ref. [START_REF] Masoero | International RILEM Conference on Materials[END_REF] ). Despite this assumption, the rate expressions that we will obtain here will display general features of coarsegrained rates that should emerge irrespectively of the molecular mechanisms considered. These features will be highlighted and discussed.

A. Off-lattice Kinetic Monte Carlo KMC stochastically explores state transitions, selecting events (here particle insertion or deletion) in proportion to their rates. The rates must be determined beforehand depending on the current configuration of particles and implicit medium. We will refer to the classical rejection-free KMC algorithm [START_REF] Serebrinsky | [END_REF] , whereby the time increment associated to any event is the inverse of the cumulative rate of all the events that are possible at that time. In this work we will consider events whose rates explicitly depend on time: this requires a slight modification of the original algorithm, as per ref. 24 Particle deletions are easy to consider. For example, the N orange particles in Figure 1a are all possible candidates for deletion, and the list of their associated rates is built by trying to remove each of them separately from the same reference configuration. In the next subsection II B we will see that the deletion rates depend on the change of interaction energy ∆U that the deletion itself would cause in the system, hence the deletion rate depends on the position of the particle to be deleted.

Particle insertions are more difficult to deal with because insertion can occur at infinite possible positions. Considering only a finite number of positions on a lattice would require extremely fine lattices, because the interactions in nanoparticle systems are typically very short ranged 25 (precisely, one would have to consider a lattice size of the order of or less than the interaction length between particles, i.e. Angstroms). A small particle displacement can thus result in substantial changes of ∆U caused by particle insertion, and therefore of the rate. We propose instead an off-lattice approach: M trial particles are first placed on a cubic lattice and, before calculating the rates, we let each trial particle move within its lattice cell volume V c to minimize its interaction energy with the other existing particles (no interaction between trial particles; see Figure 1b). The position of a trial particle after local minimization is considered as representative of all other positions within the same V c . In appendix A we will discuss how this approach ensures convergence of mechanisms and rates with the fineness of the lattice. It is worth noting that "off-lattice" in the KMC community refers to the possibility of particles to access any position in space, which is typically achieved by some sort of relaxation [START_REF] Boateng | Multiscale Model[END_REF] . This does not remove the need of creating a finite catalog of possible events, thus only a finite set of representative positions must be considered in the end. One may argue that "KMC with adaptive lattice" would be a better name for this approach, but here we will retain the conventional terminology ("off-lattice").

Each KMC step performs either a particle insertion or a deletion. Before moving to the next step all the existing particles, that thus far were kept fixed in their positions, are allowed to rearrange and minimize their interaction energy. Time is advanced only by the KMC steps. All the simulations in this work have been performed with a specifically written code that uses LAMMPS as a library. 27

B. Coarse-grained rates

The mechanisms of nanoparticle formation/dissolution in Figure 1c involve a series of chemical reaction, each causing the addition or removal of one molecular unit of solid in radial direction. We call rate r i + the rate of the i th reaction of molecular unit addition, r i -the one-unit dissolution. Transition State Theory (TST) provides expressions for the rates (per unit area) of molecular units forming on a solid surface or dissolving from it [START_REF] Lasaga | Kinetic Theory in the Earth Sciences[END_REF] :

r i ± = k B T h 1 γ * exp - ∆G a,i ± k B T . ( 1
)
k B is the Boltzmann constant, T the temperature, γ * the activity coefficient of the activated complex at the transition state (in units of area), and h the Planck constant. ∆G a,i + and ∆G a,i -are the activation free energies for the i th addition and dissolution respectively.

TST implies assumptions about the underlying kinetics, e.g. that only one reaction pathway exists for each transition, that the activated state is unique, and that reactants are at equilibrium with activated complexes. In complex liquid-solid reactions these conditions are often unmet and more advanced kinetic theories and rate expressions are needed 28 . Here we use TST because we will not refer to specific material systems for which an alternative theory is known to be preferable, and because the simplicity of TST will allow us to analytically obtain coarse-grained rates which nevertheless display some important features that will be listed and discussed later in this section. Once developed under the assumption of TST, our approach to derive coarse-grained rates will provide a methodology for future works that may be targeted to specific material systems and include a richer descriptions of their molecular-scale kinetics. ∆G a,i ± are determined by two contributions (Figure 1c): ∆G 0 a,i which depends only on the chemical reaction under consideration and on the concentration of reactants in the implicit solution, and a fraction χ (between 0 and 1) of the change of free energy ∆G eq of the system in (metastable) equilibrium just before and just after the i th reaction. ∆G 0 a,i for dissolution coincides with the activation energy of dissolution in standard state, ∆G * , whereas for solidification one must add -k B T ln(β) to account for the increased energy level of the starting configuration due to the supersaturation β of reactants in the implicit medium [START_REF] Lasaga | Kinetic Theory in the Earth Sciences[END_REF] . ∆G eq instead is due to the change of interaction energy U and particle surface Ω caused by the i th reaction (times the solidsolution interfacial energy γ i , which in general depends on the radial particle size i 29,30 ). All this leads to:

∆G a,i + = ∆G * -k B T ln(β) + χ γ i ∆Ω i i-1 + ∆U i i-1 , (2)
∆G a,i -= ∆G * + (1 -χ) γ i ∆Ω i-1 i + ∆U i-1 i , (3)
where ∆Ω i i-1 should be read as the change of particle surface area caused by the reaction that takes the particle from size "subscript" (i-1) to size "superscript" (i).

The atypical part of ∆G a + and ∆G a -in eqs 2 and 3 is their dependence on the change of equilibrium free energy ∆G eq caused by the reaction, viz. the terms following the χ parameter. A value of χ intermediate between 0 and 1 means that ∆G eq occurs progressively during the chemical reaction, and not entirely before (χ = 1) or after (χ = 0) the transition state. This is analogous to the Prandtl-Eyring model of non-Newtonian flow, which considers that the strain energy released by the activation of a deformation mechanism contributes to reducing the activation energy of the mechanism itself 31,[START_REF] Gutzow | The vitreous state[END_REF] . To understand why this aspect is important for our approach, we will use the rest of this paragraph to compare eqs 2 and 3 to classical crystal growth theory. In classical crystal growth, χ = 0 and therefore γ i ∆Ω and ∆U affect only the dissolution rate [START_REF] Lasaga | Kinetic Theory in the Earth Sciences[END_REF] . This implies that a one-molecular solidification causing a large increase of ∆U , e.g. forming a new molecular unit in the same place as an existing one, would have the same rate as a solidification reaction causing a negative, hence favourable, ∆U . This approximation is manageable in algorithms with molecular-scale resolution, in which the unit that caused the unfavourable positive ∆U would have a very high dissolution rate in the next step. The approximation is more problematic at the coarse-grained level, where one particle formation represents n molecular reactions of solidification. In this case, accepting the formation of a particle with large positive ∆U means disregarding that dissolution would immediately block the progression of n unfavourable molecular reactions. The issue is avoided in classical growth theory because net growth rates are adopted: r i,net = r i + -r i -. Net rates however are not well suitable to nanoscale simulations because they limit the effect of entropy-driven fluctuations. To understand why, notice that eqs 1, 2 and 3 with χ = 0 would lead to r i,net ∼ β -exp (γ i ∆Ω i-1 i

+ ∆U i-1 i
)/(k B T ) . Now consider for example the formation of non-interacting (∆U = 0) two-molecular metastable clusters with γ i = 0, from a solution containing single molecules. If the solution is such that β = 1 (equilibrium with respect to the growth of a large crystal), then r i,net = 0 and no cluster could form. This is incorrect because some clusters should still form driven by the entropy of mixing, which the net rates average out. This approximation is acceptable for large crystals but problematic at the nanoscale. Another issue is that net rates can be negative: this is not a problem when net rates are used in the integration of rate equations at a continuum level, e.g. in Phase Field approaches, but in a KMC algorithm all the rates must be positive. Hereafter we show that eqs 2 and 3 with χ intermediate between 0 and 1 penalise events increasing the free energy of the system while leading to coarse grained rates of particle insertion and deletion that are always positive and leave room for the entropy of mixing. Overall this is more appropriate for nanoparticle simulations than classical growth theory.

To turn the reaction rates per unit solid per unit area in eq 1 into coarse-grained rate expressions (per unit time only), we assume that particle insertion entails n straight formations of one-molecular units, whereas particle deletion entails n straight dissolution reactions. This gives:

R in = C cn V c (1 -φ loc )a 2 n i=1 r -1 i + -1 , (4)
R del = a 2 n i=1 r -1 i - -1 . (5)
R in and R del are the coarse-grained rate of particle insertion and deletion, C cn is the concentration of single molecular units of solid that we assumed to be always present in the solution, φ loc is the volume fraction of V c occupied by already precipitated solid particles (not the trial ones), and a is the linear size of the molecular unit of solid (Figure 1). Overall, C cn V c (1 -φ loc ) is the number of single-molecule nuclei in the portion of V c that is not already occupied by solid particles. Other sequences assuming both solidification and dissolution reactions behind particle insertion (or deletion) may be possible but would require more than n events, taking more time and thus having lower probability in the KMC algorithm.

Combining eqs 1 to 5 we obtain:

R in = kr * 0 βa 2 n i=1 exp χ γ i ∆Ω i i-1 + ∆U i i-1 k B T -1 , (6)
R del = r * 0 a 2 1 i=n exp (1 -χ) γ i ∆Ω i-1 i + ∆U i-1 i k B T -1 , (7)
where we defined

k = C cn V c (1 -φ loc ) and r * 0 = k B T h 1 γ * exp -∆G * k B T
. These coarse grained expressions reflect the assumptions that we listed in the introductory part of this section II. However, they address both the challenges that we presented in the introduction of this manuscript and display general features to be expected from any coarse-grained rate of nanoparticle formation and dissolution. Specifically:

• k accounts for the fact that particle insertion involves not only growth, but also nucleation (here simply captured by the constant C cn , but elaborations are possible [START_REF] Masoero | International RILEM Conference on Materials[END_REF] ). Furthermore, V c in k ensures that increasing the number of trial nuclei M to sample possible positions of particle insertion (Figure 1.b) does not lead to a systematic increase in the total rate of particle insertions

M j=1 R in j (be- cause V c = V box /M
, where V box is the volume of the simulation box).

• r * 0 sets the timescale, linking it to the chemical reactions taking place in the system via γ * and ∆G * . The resulting chemistry-controlled timescale can be many orders of magnitude larger than that of phonons in molecular dynamics simulations.

• β accounts for the chemistry of the solution and how it drives precipitation. This is important because β is what the experimentalists typically can control, e.g. the concentration of ions in a solution.

• n captures the fact that forming or deleting a large particle (large n) take longer and thus have lower rate compare to forming or deleting a small one.

• ∆Ω i i-1 and ∆U i i-1 account for changes in equilibrium free energy caused by the i th chemical reaction, or more specifically in the χ-related portion of equilibrium free energy that affects the activation energy of the reaction. ∆U is important because it captures how mechanical interactions impact mesoscale precipitation in terms of both rates and mechanisms (thanks to ∆U , R in depends on the position of particle insertion). In particular, ∆U disfavours the occurrence of any particle insertion/deletion event that causes an increase of local stress, and favours events causing stress relaxation.

Later in this manuscript we will use interaction potentials with interaction strength scaling as the surface area of the interacting particles, viz. U ∼ Ω. In this case the coarse-grained rates simplify to:

R in = kr * 0 βa 2 n i=1 exp χ γ i + ∆U in Ω k B T 8a 2 i (2i -1) 2 -1 , (8)
R del = r * 0 a 2 1 i=n exp (1 -χ) -γ i + ∆U del Ω k B T 8a 2 i (2i + 1) 2 -1 . (9)
∆U in and ∆U del are the total changes in interaction energy that would be caused by a whole particle insertion or deletion. Ω is the whole surface area of the particle, and the i-dependent fractions are the absolute value of the changes of surface area caused by the i th molecular reaction ∆Ω i i-1 and ∆Ω i-1 i for a sphere (NB: ∆Ω i i-1 is the difference in Ω between a particle of size i -1 and the same particle of size i, divided by the number Ω i-1 /a 2 of molecular reactions that occur on the particle's surface).

C. Implications for equilibrium

Equilibrium is reached when the total rate of particle insertion equals on average the total rate of particle deletion:

R in tot /R del tot = M R in /(N R del ) = 1.
The angle brackets indicate averages over all possible insertions and deletions. Considering eqs 6 and 7, this gives:

β kM N 1 i=n exp (1 -χ) γi ∆Ω i-1 i + ∆U i-1 i k B T n i=1 exp χ γi ∆Ω i i-1 + ∆U i i-1 k B T := β app = 1 . (10)
For a more convenient notation, we call µ all what comes after β, thus equilibrium entails

β app = βµ = 1.
In classical crystal growth, equilibrium of a large crystal occurs at β = 1. eq 10 shows that other terms are also important for the equilibrium of nanoparticle systems. Consider for example a scenario in which the formation of a new particle causes on average an increase of interfacial or interaction energy (positive, thus unfavourable, ∆Ω i i-1 and ∆U i i-1 ). This implies µ < 1 in eq 10 and the equilibrium condition

β app = 1 will require β = µ -1 > 1.
The µ term depends on the mesoscale morphology of the nanoparticle aggregate, thus the equilibrium of a mesoporous material is related to the mesoscale morphology too. When one experimentally sets the equilibrium of a mesoporous material forming e.g. by precipitation, the condition that is actually created in the experiment is β app = 1, not β = 1. Later in this manuscript we will refer to an experimentally measured supersaturation for the precipitation of mesoporous cement hydrates: it will be important to remember that the supersaturation from that experiment is β app and must be converted to β before using it in our rates. This conversion requires knowing µ, because β = β app µ -1 . Obtaining µ directly from eq 10 can be tedious. A simpler approach, which we will use, is to start a preliminary simulation with a sufficiently large β, which stimulates particle insertion, and when the number N of inserted particles is sufficiently large to have a representative morphology, reduce β progressively until N stays constant on average. This β corresponds to the condition β app = 1, thus equals µ -1 .

Another important implication of eq 10 concerns entropy. Consider again the example that we discussed in section II B in relation to classical growth theory: a system of two-molecular clusters with γ i = 0 and U = 0, forming from a solution of single molecules with concentration C cn . In this case, µ in eq 10 is simply kM/N and, when β = 1, the equilibrium condition becomes

N = kM = C cn V box (1 -φ loc )
. This corresponds to maximising the entropy of mixing of an ideal binary mixture, i.e. the number of two-molecular clusters N equals the number of single molecules in the volume of solution that is not occupied by two-molecular clusters. This means that, differently from the previously discussed case of classical growth theory using net rates, our approach accounts for the entropy of mixing as an additional driving force for particle formation and dissolution.

D. Mesh effects of the trial particle lattice

We pointed out that the total insertion rate does not scale with the number of trial particle lattice cells M (Figure 1b), because V c in k, in eq 6, scales as M -1 . Other errors however can occur due to the finite size of the trial particle lattice, in particular because only one position for each lattice cell is taken as representative of the interaction energy anywhere in that cell (section II A). This approximation is good as long as one position is overwhelmingly favourable in terms of energy compared to all other possible positions within the cell, and as long as the trial particle can find that position. In appendix A we discuss in more detail the possible issues related to mesh effects, concluding that they should be evaluated case by case depending on the specific process to be simulated and on the quantities of interest. In appendix A we show how we evaluated mesh effects for the cases considered in this manuscript, concluding that a trial lattice with linear size equal to one particle diameter is sufficient to describe correctly the overall mesoscale mechanism of precipitation, which is our main focus here, even if the absolute value of the rates may be overestimated by several orders of magnitude.

E. Condensation of parameters

The coarse-grained rates in eqs 6 and 7 contain many parameters. This is a necessary consequence of having based our coarse graining on a mechanistic description of molecular scale processes. Nevertheless, these parameters are coupled to the extent that approximate insertion and deletion rates averaged over all possible deletions and insertions can be expressed using three independent parameters only. To do this we take χ = 0.5 and approximate ∆Ω i i-1 and ∆U i i-1 as Ω/( 4 3 πn 3 ) and ∆U /( 4 3 πn 3 ), i.e. all [START_REF] Van Vliet | [END_REF] 3 πn 3 chemical reactions involved in the formation or dissolution of a particle contribute with same change in particle surface and interaction energy. This leads to:

R in ≈ r * 0aγn • α βγ • α U , (11)
R del avg ≈ r * 0aγn • α -1 U . ( 12
)
r * 0aγn = r * 0 a 2 n -1 exp γΩ 2k B T 4 3 πn 3
is a kinetic constant that sets the timescale but, being the same for insertion and deletion, it does not affect the mechanism of precipitation nor the evolution of the rate with time.

α βγ = β exp -γΩ k B T 4 3 πn 3
depends only on the chemistry of the system, whereas α U = exp ∆U 2k B T [START_REF] Van Vliet | [END_REF] 3 πn 3 accounts for the relative magnitude of free energy changes caused by mechanical interactions and thermal fluctuations. ∆U is the same but opposite in sign for average particle insertion and deletion. Overall, this means that our model is governed by three main elements: the kinetic timescale, the chemical driving force, and the change in free energy caused by mechanical interactions. In section III we look into the effect of the mechanical interactions, α U , while section IV explores also the chemical drive α βγ .

III. RESULTS: MECHANISMS OF ORDERED PRECIPITATION

We consider a 3D system of particles that precipitate and aggregate on an insoluble substrate (orange and grey particles in Figure 2). To favour the formation of ordered precipitates, we build the substrate as an ordered fcc array of spherical particles and we allow for the precipitation of spherical particles with only one possible diameter, same as for the substrate particles. In eqs 11 and 12 we set r * 0aγn = α βγ = 1 and focus only on: (i) the strength of inter-particle cohesion compared to the intensity of thermal fluctuations α U , and (ii) the strength of the interactions between precipitate and substrate.

In the coarse-grained rate expressions we use χ = 0.5 and pairwise Lennard-Jones interactions:

U ij (r) = 4 σ r 2α - σ r α , (13)
with α = 6 and σ = 1. The interaction strength is initially set to 70k B T , which is in the range of covalent bonds. In a maximum-density fcc packing, monodisperse spherical particles have 12 neighbours. On a growing surface we can expect that half of these contribute to the change of interaction free energy, ∆U = 6ε. This gives α U ≈ exp(50), having taken n = 1 here. The particlesubstrate interactions are initially identical to those between precipitating particles.

The mechanisms emerging from these conditions is the layer-by-layer precipitation in Figure 2a. This is an enthalpy-driven mechanism governed by the strong α U and particle-substrate interactions. Forming a new particle on top of a fully formed layer takes much longer than adding a particle to a horizontally growing layer, because of the reduced number of neighbours, approximately 6 on an growing layer, and 4 on a new layer, giving a difference of ca. ∼ exp (17) in terms of relative rate. Consistently, the temporal evolution of the process in Figure 2 displays horizontal plateaus which are the waiting time to nucleate a new layer, and sharp steps up corresponding to the rapid horizontal growth of a layer. This is known in the literature as Frank-van der Merve regime. [START_REF]Physics of Crystal Growth[END_REF] Predictably, the precipitating layers preserve the fcc order.

Decreasing α U to ∼50, corresponding to ε ≈ 6k B T (in the range of van der Waals interactions), leads to the irregular layering mechanism inFigure 2b. The weaker interactions of a particle forming on a new layer, which favours the mechanism of horizontal layer growth, are now statistically in competition with the smaller number of particles that can form on the perimeter of the growing layer, compared to those that can form on its surface and start a new layer. Therefore the timescales of layer growth and nucleation are similar, the precipitating domain gets rougher, and the temporal evolution of the process is smoother compared to the previously discussed stepwise trend.

A further decrease of α U to ∼10, still in the van der Waals range, leads to an even rougher layer accompanied also by the homogeneous precipitation of particle clusters in the bulk solution (Figure 2c). Even if α U is still significantly larger than 1, homogeneous precipitation occurs because the number of possible locations for particle precipitation in the bulk is much larger than the number of locations on the layer surface: this is an entropic effect. All characteristic timescales are of similar magnitude, thus the reaction proceeds linearly in time and the stepwise trend is lost.

We now move the focus onto the precipitate-substrate interactions. We go back to the original α U ≈ exp(210) and reduce the substrate-particles interactions by 50%. This triggers the mechanism of island growth in Fig- ure 2d, known as the Volmer-Weber regime. [START_REF]Physics of Crystal Growth[END_REF] The rate of forming an isolated particle on the substrate is much smaller than the rate of forming a particle next to another precipitated one. Because of this, new particles aggregate preferentially with as many existing ones as possible and not as much on the substrate, which favours the 3D growth of the aggregated domain rather than horizontal growth of layers. Furthermore, precipitating the first FIG. 2. Effect of substrate-particle and particle-particle interactions on the precipitation of ordered domains. The size of the cubic simulation boxes explored are between 20 3 and 50 [START_REF]Nanostructured Soft Matter: Experiment, Theory, Simulation and Perspectives Zvelindovsky[END_REF] particle diameters. Snapshots created with OVITO [START_REF] Stukowski | [END_REF] .

particle takes a long (induction) time, and then the process accelerates as more and more particles are formed. The structure of the island appears to have mainly fcc symmetry with a small fraction of hcp regions.

The same conditions that led to the island mechanism applied to an initial configuration with amorphous substrate structure trigger the interesting mechanism in Fig- ure 2e. At first, the aggregated domain grows disorderly as a hemisphere. At some point however the local arrangement of some precipitated nanoparticles happens to display some ordering, which triggers the much faster growth-by-aggregation of an ordered domain on top of the amorphous one. This delayed ordered domain displays a comparable incidence of hcp and fcc regions. The disorder of the substrate has two further implications: (i) it generates a sufficiently rough energy landscape to allow for multiple particles to form on the substrate before the large domain eventually takes over, and (ii) it removes separation of timescales between first precipitation and island growth, hence the induction time discussed for the ordered island mechanism disappears. The appearance of the ordered domain corresponds then to a sharp acceleration of the precipitation.

In all the mechanisms discussed in this section (and in the next one too) we considered a substrate made of monodisperse particles with same size as the new precipitating particles. Releasing this assumption, e.g. to consider coarsely patterned or very finely smoothed substrates, may lead to a substrate that is incompatible with the immediate ordering of new precipitating particles. This is conceptually similar to the case of amorphous substrate discussed above, thus if the precipitating particles are monodisperse we would still expect delayed ordering with details that depend on the morphology of the substrate. Our simulations may thus help predict the important effect of substrates in key technological applications such as thin film deposition.

Overall, our simulations captured some precipitation mechanisms for ordered mesoporous domains, predicting the corresponding temporal evolution of the process. We have shown that interaction strength and substrate structure are key in determining the mechanism, and that entropic effects can be important too.

IV. RESULTS: MECHANISMS OF AMORPHOUS PRECIPITATION

Here we start with an amorphous insoluble substrate, 20 nm thick and made of nanoparticles with diameter of 10 nm (not participating to the KMC insertion-deletion). We allow for precipitation of monodisperse particles, but then we slightly (±5%) and randomly alter the size of a newly inserted particle to avoid ordering. The parameters for the coarse-grained rates in eqs 6 and 7 are matched to the case of cement hydrates (more specifically calcium-silicate-hydrate, C-S-H, which is the main hydrate in ordinary cements), precipitating from aqueous ionic solution and forming an amorphous mesoporous phase. We simulate 4 precipitation mechanisms that appear in the literature on cement and we discuss the associated kinetics predicted by our simulations. We also discuss the effect of an evolving chemical environment, viz. a supersaturation β that changes with time.

The details of the parametrization for cement hydrates are in the appendix: here we present only the most relevant features:

• The diameter of the precipitating particles is D = 10 nm. 35 The linear size of a C-S-H molecular unit is a = 0.65 nm, 36,37 thus forming a particle requires n = 8 units in radial direction.

• We use γ = 87.6 mJ/m 2 as interfacial energy between C-S-H and solution 38 .

• For the interactions between particles we use the spherical Lennard-Jones-like potential in eq 13 with σ = D and α = 14. 39,40 We relate ε to the solidsolution interface that is removed when particles are in contact: ε = γΩ/6 where 6 is half the number of first neighbours in an fcc crystal and Ω is the particle surface area. Since ε ∼ Ω we can use eqs 8 and 9. The resulting energy scale is ε = 660 kCal/mol which is close to the value indicated by recent molecular simulations. 41 • The supersaturation β in a cement solution depends on coupled processes of cement dissolution, ion transport, and hydrates precipitation. Here we focus only on precipitation, thus the temporal evolution of β(t) must be provided as an input. To clarify the role of the chemical drive, we consider two extreme cases. Firstly, β = 100 and constant in time, which represents a solution at equilibrium with very soluble cement grains. Such a large value of β ensures that the nanoparticles will precipitate indefinitely as long as there is space available. Secondly, we start with a very concentrated ionic solution which gets rapidly depleted. This type of supersaturation evolution for a cement solution has been simulated in ref. 38 (high opacity dataset; see Figure 3a). The reference actually provides what we called β app in section II C. With preliminary simulations as described section II C we obtained µ ≈ 1/3, which we used to convert β app (t) into β(t) for our coarse grained rate expressions.

• The kinetic constant r * 0 is unknown for cement hydrates precipitation, therefore we treat it as a free parameter. Notice that in our simulations, in which only one type of particle can form or dissolve, r * 0 only act as a linear scaling factor for the overall rate. Thus we can set r * 0 to obtain a timescale of ∼24 hours that is meaningful for the precipitation of cement hydrates, without affecting the shape of the rate curve as a function of time.

For all the simulations we used lattices of trial particles with spacings between 0.5D and D. To speed up the calculations, we did not calculate the local solid fraction of precipitated solid φ loc in eq 6; instead we assumed a constant value of φ loc = 0.5 everywhere, which is reasonable for the dense assemblies of almost monodisperse nanoparticles that we will obtain. We used the same values of ε, γ, and β in all simulations except one case in which we consider weaker interactions from atomic force microscopy experiments. 42 As for the cases of ordered precipitation, also here we use χ = 0.5 in the coarse grained rate expressions.

A. BNG (or Cahn) mechanism

The Cahn mechanism describes the heterogeneous precipitation of a phase on a substrate starting from defects that provide favourable sites for initial nucleation. 14 The new phase grows hemispherically in 3D, similar to the island mechanisms in Figure 2d, with a constant rate of radial growth (extensions to anisotropic growth and non-constant growth rates exist 38 ). 3D growth causes acceleration of the total reaction rate, i.e. the total volume of precipitated solid per unit time. This continues until adjacent precipitated domains growing from different sites start to impinge, which leads to a peak in the total rate. If there is enough bulk solution space left, the domains will eventually form a thick layer and grow perpendicularly to the substrate with constant total rate.

Today the Cahn mechanism is probably the most commonly used to describe the precipitation of cement hydrates on the surface of slowly dissolving cement grains.

In the cement literature it is better known as the Boundary Nucleation and Growth (BNG) mechanism. 14 In a recent paper the rate of radial growth of the precipitating phase has been expressed as a function of the supersaturation of the solution with respect to C-S-H precipitation. 38 We will use of the supersaturation data and BNG results in that paper to inform our simulations and discuss the results.

The initial condition for our simulation is the substrate in Figure 3c0, which is periodic in plane and represents a cement grain. Cahn's model assumes the presence of favourable nucleation sites on the substrate. To mode them we consider an interaction strength between substrate and cement hydrates particles that is 75% smaller than the ε between hydrate particles, except for some substrate particles within a 30 nm region (the favourable site) which have same ε as for hydrate-hydrate interactions. The density of favourable sites on the substrate is ∼10 µm -1 , from ref 38 . We also run additional tests without the favourable sites and, as for the already-shown island mechanisms in Figures 2d,e, we still obtained the formation of 3D hemispherical domains growing and impinging. Our simulations therefore could be used to clarify the relationship between density of nucleation sites in the Cahn model, substrate morphology, and substrateprecipitate interactions. This however goes beyond the scope of this manuscript and is left to future work.

The simulation results are shown in Figure 3: r * 0 was set to place the rate peak at ∼5 hours, leading to r * 0 = 43.33 ± 0.125 s -1 µm -2 for variable β and r * 0 = 46.66 s -1 µm -2 for constant β. During the first 4 hours, the 3D growth of precipitated domains in Figure 3c1 allows an increasing number of trial particles to find favourable ∆U positions on the surface of the domains. Trial particles with favourable ∆U contribute with high R in to the total rate, and this explains the initial acceleration, which is caused by the mechanical interactions and therefore appears in both cases of variable and constant β. The rate peak occurs because the precipitated domains impinge laterally on each other, as shown in Figure 3c2. The impingement reduces the surface area of domains, changing the growth mechanism from 3D hemispherical to layered, as in Figure 3c3. Layered growth implies that the part of the rate related to mechanical interactions is constant on average, thus the temporal evolution of the rate is controlled by β (at least in our far-from-equilibrium conditions, when particle deletions are infrequent). This explains why the post-peak regime of our simulated rate is constant for the case with constant β and decreases for the case with decreasing β(t). Overall, our results show that initial acceleration and post-rate-peak regime of a BNG mechanism are respectively determined by the mechanical interactions between nanoparticles and by the temporal evolution of the solution chemistry.

B. Layered precipitation mechanism

In the cement literature, the layered growth mechanism is not considered as a likely one for hydrates precipitation. It is sometimes used for its simplicity in models that consider mass and volume balance equations. [START_REF] Bishnoi | [END_REF][45][46] Here we consider this mechanism mainly to show how its associated kinetics, predicted by our simulations, compares with the previous one from BNG. We use the same parameters as in the BNG simulations, but we remove the initial favourable sites and set the substrate-hydrates interaction strength equal to the hydrate-hydrate one ev-erywhere. This is analogous to the case of layered growth in Figure 2a, but now the amorphous substrate and the small size polydispersity of the particles prevent ordering, as shown in Figure 4a. This leads to the same situation as in the final stages of BNG (cf. Figure 3c3), therefore the rate reflects the temporal evolution of β: constant when β is constant, and decreasing when β is taken from Figure 3a (r * 0 = 4.785 s -1 µm -2 and r * 0 = 2.9 s -1 µm -2 respectively for constant and variable beta). Therefore layered precipitation does not leave any room for initial acceleration, which disagrees with the experiments. 47

C. Gel-like mechanism

The gel-like mechanism has been proposed because during the first hours and days of cement hydrates precipitation the interactions may be governed by rather weak non-covalent forces. 18 Atomic force microscopy (AFM) experiments 42 have quantified these forces showing that the resulting energy scale is approximately one order of magnitude weaker compared to what we are using here. Such weak interactions may lead to progressively densifying gel-like structures, as discussed in the literature 18,48,49 . To test the gel-like mechanism we use same initial conditions as in section IV B but with inter- actions taken directly from the above-mentioned AFM experiments 42 (data for a 10.3 mM concentration of calcium hydroxide in solution, which is the largest one reported in ref. 42 ). We fitted the AFM interactions using eq 13, obtaining α = 16 and ε = 70 kCal/mol, much smaller than the 660 kCal/mol that we used so far. Figure 4b shows that these weak interactions induce homogeneous bulk precipitation and aggregation.

The consequences in terms of rate are interesting. In Figure 4b we see that the variable β leads to a steadily decreasing rate similar to the case of layered precipitation (here r * 0 = 10 s -1 µm -2 ). The rate is thus controlled by β(t) because the interactions are so weak that even a small variation of β will overcome them. If β is constant the aggregation of small clusters in the bulk space induce instead some acceleration (still with r * 0 = 10 s -1 µm -2 ). This acceleration culminates in a rate peak when the bulk space starts to get filled. The resulting shape of the rate curve is similar to that from BNG, thus consistent with the experiments. There are however two problems: firstly, it is unlikely that β stays constant during the precipitation of cement hydrates (although this scenario cannot be excluded a priori near the cement surface). Secondly, the position and magnitude of the rate peak in Figure 4b depends approximately linearly on the size of the simulation box perpendicular to the substrate, as shown in the inset. This happens because the weak interactions imply that the rate of precipitating near existing particles is similar to the rate of precipitating far away in the bulk. This causes a direct dependence of the total rate on the volume of the simulation box. The box size is a parameter that can be changed to represents non-flocculated cement solutions with different water-to-cement mass ratios (w/c). Experiments show that even orders of magnitude differences in w/c of stirred suspensions (to minimise flocculation) do not affect much the rate, 50 , which argues against the gel-like mechanism. This problem however would be solved if the hydrates could precipitate only within a limited distance from the cement surface: a "reaction zone" maybe determined by ion diffusion and mean collision free path. 49,50 Overall, the gel-like mechanism displays an interesting physics, to be explored further for the precipitation of cement hydrates and relevant also for other materials. For the case of cement, it presents problems that may be addressed in more detail by combining our simulations with models of ion transport and cement dissolution. 38

D. Avrami mechanism

The Avrami, or JMAK, mechanism is analogous to BNG but the initial nucleation takes place in the bulk solution and not on a surface. [START_REF] Avrami | [END_REF] The resulting kinetics is similar to the BNG kinetics: initial acceleration during the spherical growth of precipitated domains and rate peak caused by the impingement of adjacent domains. Unlike BNG, however, the Avrami mechanism entails exponential post-peak deceleration due to the filling of the bulk space, which is too fast compared to experimental results on cement. 46 Furthermore, as we already discussed for the gel-like mechanism, bulk space filling entails a strong relationship between precipitation rate and w/c. However, recent results from nanoparticle simulations suggest that the Avrami mechanism may still be important at the interface between different cement grains when the available bulk space is controlled by the flocculation of cement grains and not by the w/c. 19 In our simulation, we obtain the Avrami mechanism by taking the same inputs as for the BNG mechanism discussed above and removing the substrate. The rates in Figure 4c display acceleration for both constant and variable β (r * 0 = 52.635 s -1 µm -2 and r * 0 = 49.764 s -1 µm -2 respectively). This is expected because acceleration is governed by the interactions and not by β. However, as typical for the Avrami mechanism, the deceleration is very sharp unlike the gradual slowdown in BNG. The deceleration is particularly abrupt in our results because our simulation box is small and therefore only one aggregated domain of precipitated particles can form in it, impinging very regularly with its periodic images. A bigger bulk space would allow for a statistics of domains at various distances and impinging at different times, making the deceleration a bit more gradual although still exponential and therefore too fast.

V. CONCLUSION

We have presented Kinetic Monte Carlo simulations for various mechanisms of nanoparticle precipitation leading to mesoporous aggregates. Our coarse graining, based on molecular mechanisms of growth and dissolution, provides effective rates of nanoparticle insertion and deletion that involve many parameters. However, all these parameters have physical meaning, thus are measurable experimentally and/or calculable by molecular or ab initio simulations. We have shown that our model can be reduced to only three independent average parameters, which fully determine precipitation mechanism and associated rate. For simplicity, we considered only monodisperse particles with constant size during a simulation. Thus we did not explore the richness of morphologies that could stem from the molecular growth mechanism assumed when deriving the rates. Nevertheless, assuming a molecular mechanism of precipitation rather than postulating a rate expression is essential to the molecular-to-mesoscale coarse-graining approach that we advocate. Choosing growth rather than other possible ones (e.g. Classical Nucleation Theory) has the advantage of setting the ground for follow up simulations allowing for changes in particle sizes and shapes, towards complex morphologies such as foils and needles.

Our kinetic approach is a step forward compared to existing equilibrium-based simulation (e.g. Grand Canonical Monte Carlo 18,19,53 ) in which the temporal out-of-equilibrium evolution of the system is neglected or described qualitatively via heuristic timescales. Three features are particularly noteworthy: (i) our KMC algorithm treats precipitation off-lattice, ensuring a consistent sampling of the most probable positions for particle insertion; (ii) our coarse-grained rates combine free energy contributions from both chemical reactions and mechanical interactions, allowing our simulations to address chemo-mechanical processes; (iii) our proposal to alter the activation energies of precipitation and dissolution based on changes of surface and interaction energies preserves entropic effects that can be important for the precipitation of nanoparticle aggregates.

When applied to the precipitation of ordered domains, our simulations captured the effect of decreasing the strength of mechanical interactions with respect to thermal fluctuations. An associated change of precipitation mechanism was predicted, from perfect layering, to rough layering, and to homogeneous cluster nucleation. The simulations also showed that weak layer-precipitate interactions trigger island formation with a long induction period, and that an amorphous substrate initially generates an amorphous aggregate of precipitated particles that is eventually overtaken by crystallization, similar to how amorphous clusters serve as precursors to nanocrystal formation at the molecular scale. [START_REF]Physics of Crystal Growth[END_REF] Finally we considered the precipitation of amorphous mesoporus cement hydrates, for which we parametrized our model using data from molecular simulations and experiments, plus reasonable assumptions when data were not available yet. A substrate with favourable nucleation sites leads to a BNG mechanism, which today is regarded as a model for the precipitation of cement hydrates. Removing the favourable sites and considering strong substrate-precipitate interactions triggered a disordered layering mechanism, with unrealistic rate evolution. Starting the simulations without substrate triggered an Avrami mechanism with two shortfalls: a max-imum rate that corresponds to the filling of bulk space, which disagrees with experiments on cement pastes at high water-to-cement mass ratio, and a too abrupt deceleration. Finally, weaker interactions induce gel-like precipitation in the bulk solution, with rate evolution that under certain condition could be similar to those given by the BNG model, but that also entail a strong dependence on the water-to-cement ratio. In all the cases that we considered, acceleration was determined by the mechanical interactions while the chemical drive became more important in the later stage of deceleration.

Overall this work is a first step in developing bottomup, chemo-mechanical, coarse-grained simulations in mesoscale chemistry. Many interesting developments are possible, e.g. coarse-grained expressions for non-classical mechanisms of nucleation and growth, partial particle growth and dissolution to capture autogenous stresses, multiple options for particle shapes, sizes, and chemical compositions, and diffusive processes in the implicit solution to address ripening. Our coarse-graining approach might also be adapted to mesoporous materials that do not form by nanoparticle precipitation and aggregation, e.g. zeolites: this may be achieved by adding variables beyond size and shape, which describe the evolving molecular structure inside a particle (e.g. the number and orientation of molecular layers inside a "particle", or more appropriately a discrete unit representing a finite volume of a crystalline material). Even without these features, our simulations already capture some interesting aspects of the complex kinetic interplay between chemical drive and self-organizing mechanical interactions during nanoparticle precipitation. Continuing efforts in this direction may eventually enable mesoscale simulations that explicitly consider real experimental design variables, such as solution chemistry and temperature, using them to guide the experimental development of mesoporous materials with tailored morphologies and properties.

In section II A we explained that before computing R in we let each trial particle minimise its interaction energy with the existing particles by moving locally within its lattice cell, and that the ∆U at the end of this local minimisation is taken as constant everywhere in the cell. If the trial particle finds the global energy minimum in the cell, this approximation leads to an overestimated rate. If instead the particle finds a local minimum in the cell, the rate may be overestimated or underestimated. In terms of mesoscale mechanism of precipitation, this approach gives the correct result if:

• Each trial particle finds a global energy minimum within its cell, and there is only one global minimum in each cell. This is always true if the lattice cell size is smaller than the characteristic width of asperities on the interaction energy landscape.

• Other local minima in the cell are much less probable than the global minimum. This is true for high-energy interactions, because the probability of a state is proportional to its Boltzmann factor.

FIG. 5. Mesh effect on a 1D system: array of nanoparticles forming from an implicit medium.

• The timescale of mechanical vibrations is much smaller than that of chemical reactions. This justifies the assumption of mechanical equilibrium and thus the use of energy minimisation.

These conditions imply that our approach averages out entropy-driven fluctuations at length-scales below the lattice cell size, so one must be careful when such fluctuations can determine the assembly mechanism, e.g. in the nucleation of defects during the growth of an ordered domain of weakly interacting particles. Furthermore, if the interaction energy landscape is smooth (low-energy and medium-ranged interactions), homogeneous and heterogeneous particle precipitation may compete and the error on the rate may bias the mesoscale mechanism. However, in this case the small energy scale would imply small errors in the rate, mitigating the issue.

If the lattice cell size tends to zero (and M to infinity) our method converges to the correct rate and mechanism. This is shown in the 1D example in Figure 5: starting from an initial particle an array is formed via strong Lennard-Jones interactions that favour particle formation near the rightmost existing particle. This example has analytical solution, whose value for a given set of parameter is indicated in Figure 5; for brevity we omit the details of the analytical solution, which results from integrating eq 6 over r with ∆U (r) of Lennard-Jones type. Figure 5 shows that our simulations converge to the analytical solution as the lattice spacing tends to zero. In this example, the conditions are such that each trial lattice cell contains only one energy minimum and therefore the rate is systematically overestimated (and the insertion time underestimated) as the lattice spacing is increased. This trend breaks down at lattice spacing above 2.5σ, which is the interaction threshold that we used in the example. At this point, the lattice spacing is so large that the trial particles may not "see" the rightmost existing particle. Hence the mesoscale mechanism changes from sequential chain growth by particle addition to homogeneous particle insertion anywhere in the box.

In 3D systems, the additional dimensions compared to the 1D case can lead to more local minima in the interaction energy landscape. Therefore in Figure 6 we show the results of our tests of trial particle lattice size effect for one of the cases that we considered in the manuscript. Specifically, we consider monodisperse particles precipitating on an ordered fcc substrate. The interactions are the strong and narrow Lennard-Jones interactions that we used in section IV. We use a strong supersaturation of the solution, β = 1000, in order to avoid particle deletion and focus only on precipitation. Under these conditions, the expected precipitation mechanism is layer-by-layer, as in Figure 2a. Figure 6 shows that the correct mechanism is indeed captured as long as the size of the trial lattice cells is equal to or slightly bigger than the particle size (10 nm). Similar to the 1D example,the rate is increasingly overestimated (thus t 300 is increasingly underestimated) as the lattice cell size is increased, until the cell becomes so large that the correct mechanism is FIG. 6. Mesh effect on a 3D system similar to those that we used in the manuscript. The particle diameter is 10 nm. t300 is the time to form 300 particles, viz. 3 layers since the substrate surface is made of 10 × 10 particles. The trial particle lattice is cubic, with linear size corresponding to the lattice spacing.

missed and the rate drops.

Overall, Figures 5 and6 show that the correct mesoscale precipitation mechanism can be captured already with a quite coarse lattice of trial particles, even if the rate may be overestimated by several orders of magnitude. This observation should be checked case by case, because there are scenarios which may require fine lattices also to capture the correct mechanism. For example, complex interaction potentials may add roughness to the energy landscape and thus require fine lattices. Also, if the focus of the simulations is on local structural details and not on the overall precipitation mechanisms as in this manuscript, it is possible that the local details may change with the lattice size even if the overall precipitation mechanisms is unaffected. For example, the amorphous precipitation mechanisms in section IV may persist even if locally the structure may change, e.g. in the radial distribution function). Furthermore, the lattice cell size can bias the overall precipitation mechanism when multiple competing events can occur, e.g. dissolution of particles of type A and insertion of particles of type B. Therefore our overall recommendation is that one should always run preliminary tests like those in Figures 5 and6 but specific to the system under consideration and looking at the quantities of interest, judging what lattice cell size is suitable and whether correction factors should be applied to the insertion rates in order to mitigate biases. Cement hydrates are produced by chemical reaction between cement powder and water. [START_REF] Taylor | Cement Chemistry[END_REF] Ordinary cements are largely made of calcium silicates that dissolve in water releasing calcium and silicon ions. The dissolution supersaturates the solution with respect to the precipitation of cement hydrates, mostly calcium-silicate-hydrate (C-S-H) and calcium hydroxide Ca(OH) 2 in ordinary pastes. Experiments on the nucleation of cement hydrates showed that the assumption of critical nuclei being as small as one molecular unit is actually reasonable for this material 47 . The hydrates form a mesoporous phase that is often represented as an assembly of nano-units with characteristic size of 5 nm. 35 Here we assume a particle diameter D=10 nm to use fewer particles at the 500 nm scale.

We set the concentration of nucleation sites in solution C cn equal to the concentration of all ions in solution. In a Ca(OH) 2 saturated cement solution the concentration of Ca ions is ∼20 mmol/L and that of Si ions is micromolar. [START_REF] Chen | [END_REF] Consistently we take C cn = 0.012 nm -3 . For the molecular unit size we consider a C-S-H molecular unit, thus a = 0.65 nm. 36 Experiments on C-S-H precipitation suggest that the diameter of critical C-S-H nuclei is sub-nanometric, thus we consider critical nuclei with diameter a and we set their nucleation time t N to zero as mentioned in Section II B. Particle insertion thus implies growing a nucleus from radius a/2 = 0.325 nm to 5 nm by adding n = 8 molecular unit of cement hydrates with linear size ∼a. We work at room temperature T = 298 K and for the scaling of the solid-solution interfacial energy with the particle size we refer to Tolman's theory, whereby

γ i = γ 1 + 8δ T Di-1+Di -1
. γ is the interfacial energy for infinitely large particles, D i is the particle diameter after i growth reactions in radial direction, and δ T is the Tolman's length, a system-specific parameter whose value is usually between 0.1 and 0.35 nm. 29,30 Here we take δ T = a/2 = 0.325 nm and γ = 87.6 mJ/m 2 as in ref. 38 For the interactions between particles we use the spherical Lennard-Jones-like potential in eq 13. This type of potential, sometimes with an additional shoulder, is commonly used in nanoparticle models of cement hydrates. 18,39,40,53 Following literature on cement hydrates 39,40 we take σ = D and α = 14 (other works used α = 12 obtaining similar results. 18,53 ) The energy scale ε must necessarily depend on the size of the interacting particles. It has been shown that a ε scaling as the particle volume can describe well the mechanical stiffness of cement hydrates. 39,40 Here however the focus is not on mechanical deformations but rather on energy changes due to placing particles near each other. For this scenario it is more appropriate to relate ε to the solidsolution interface that is removed when particles are in contact. Thus we set ε = γΩ/6 where 6 is half the num-ber of first neighbours in an fcc crystal. The ε ∼ Ω scaling allows us to use the coarse-grained rate expressions in eqs 8 and 9, and leads to an energy scale = 660 kCal/mol which is close to the value indicated by recent molecular simulations. 41 We will use the same values of , α, and γ in all simulations presented hereafter except for one case in which we consider weaker interactions from atomic force microscopy experiments. 42 The supersaturation β is a key parameter in our simulations. In a real cement solution it depends on the coupled processes of cement dissolution, ion transport, and hydrates precipitation. Here we focus only on precipitation, thus the temporal evolution of β(t) cannot be predicted and rather has to be provided as a separate input. To clarify the role of the chemical drive, we will consider two extreme cases for β. In the first case we consider a β = 100 constant in time; such a large value ensures that the nanoparticles will continue to precipitate indefinitely as long as there is space available. This condition represents a scenario in which the solution is in equilibrium with very soluble cement grains, with the solution being immediately replenished of ions that are consumed to form hydrates. In the second case we consider a solution that starts with ion concentrations well above the equilibrium for precipitation and that gets rapidly depleted, almost reaching equilibrium with the precipitates in few hours. Ref. 38 describes a similar scenario for cement hydrates (there called "high opacity"), providing both the precipitation rate for a BNG mechanism and the associated supersaturation of the whole mesoporous phase treated as a homogeneous continuum, which we call β app and reproduce in Figure 3a. In Section II B we discussed how the ∆U and ∆Ω terms involved in the precipitation of mesoporous materials generate a mismatch between β and β app . Therefore we must convert β app (t) in Figure 3a to an equivalent β(t) for the molecular reaction to use in eqs 8 and 9. To do this, we assume a linear relationship β = κβ app and we exploit the fact that when β app = 1 our simulations should describe a mesoporous system in which the average number of nanoparticles does not change. Therefore we take a domain of precipitated nanoparticles representing the cement hydrates and look for the value of β that leads to such mesoscale equilibrium: we call this β eq . By definition, β eq corresponds to β app = 1, thus we can set κ = β eq and use this κ to consistently convert β app in Figure 3a to the β(t) that our model requires.

FIG. 1 .

 1 FIG. 1. (a) Simulation box containing implicit solution (not shown), insoluble substrate (grey) and several particles (orange) that interact via a distance-dependent potential ∆U (r). (b) Same configuration with trial particles (yellow) that can move locally within their Vc and interact with existing particles but not with each other. (c) Molecular process behind nanoparticle formation: n molecular units of solid form radially via chemical reactions that change the chemo-mechanical equilibrium free energy of the system ∆Geq and have an associated activation energy made of two contributions: ∆Ga,0, which depends on the chemical reaction and on the solution chemistry, and a fraction χ of ∆Geq (more details in the text).

FIG. 3 .

 3 FIG. 3. (a) Temporal evolution of the supersaturation of a cement solution with respect to the precipitation of mesoporous cement hydrates. 38 (b) Comparison of simulated precipitation rates for time-dependent and constant β 43 with BNG result from ref 38 . (c) Our simulated BNG mechanism of precipitation: dark blue areas within the precipitated domains indicate local compressive virial stress; the rest is under tension.

FIG. 4 .

 4 FIG. 4. Mechanisms of amorphous precipitation parametrized for mesoporous cement hydrates: (a) layered, (b) gel-like, and (c) Avrami. The inset in (b) shows the box size effect on the gel-like mechanism with constant β. The BNG result from ref. 38 are shown to favour comparisons.

The authors thank Matthieu Vandamme and Sidney Yip for valuable discussion.


Holger Bru ¨ggemann

Elzbieta Brzuszkiewicz

Diana Chapeton-Montes

Lucile Plourde

Denis Speck

Michel R Popoff
email: mpopoff@pasteur.fr

Genomics of Clostridium tetani

Keywords: Clostridium tetani, Tetanus neurotoxin, Genomics, Plasmid

Genomic information about Clostridium tetani, the causative agent of the tetanus disease, is scarce. The genome of strain E88, a strain used in vaccine production, was sequenced about 10 years ago. One additional genome (strain 12124569) has recently been released. Here we report three new genomes of C. tetani and describe major differences among all five C. tetani genomes. They all harbor tetanus-toxin-encoding plasmids that contain highly conserved genes for TeNT (tetanus toxin), TetR (transcriptional regulator of TeNT) and ColT (collagenase), but substantially differ in other plasmid regions. The chromosomes share a large core genome that contains about 85% of all genes of a given chromosome. The non-core chromosome comprises mainly prophage-like genomic regions and genes encoding environmental interaction and defense functions (e.g. surface proteins, restriction-modification systems, toxin-antitoxin systems, CRISPR/Cas systems) and other fitness functions (e.g. transport systems, metabolic activities). This new genome information will help to assess the level of genome plasticity of the species C. tetani and provide the basis for detailed comparative studies.

Introduction

Clostridium tetani is a Gram-positive, spore-forming and anaerobic rod-shaped bacterium which produces a potent neurotoxin, the tetanus toxin (TeNT). C. tetani forms spores which appear as translucent terminal enlargements having the general appearance of drumsticks. Spores are resistant to desiccation, heat and exposure to oxygen, and thus C. tetani can survive in environments such as soil for long periods of time [START_REF] Popoff | Botulism and tetanus. The prokaryotes: human microbiology[END_REF]. More rarely, C. tetani can be found in the intestinal flora of animals and humans. Usually, C. tetani enters an organism through an open wound. Spores germinate and C. tetani grow in anaerobic conditions in necrotic tissues.

TeNT is synthesized at the end of exponential growth and is released in the surrounding tissues. TeNT enters motorneurons and is retrogradely transported to the central nervous system. Then, TeNT enters inhibitory interneurons and blocks the release of neurotransmitters (glycine, GABA) upon proteolytic cleavage of the SNARE protein VAMP2 [START_REF] Bohnert | Tetanus toxin is transported in a novel neuronal compartment characterized by a specialized pH regulation[END_REF]3]. TeNT is responsible for all the symptoms of tetanus, characterized by spastic paralysis (muscle rigidity and spasms) [START_REF] Hassel | Tetanus: pathophysiology, treatment and the possibility of using botulinum toxin against tetanus-induced rigidity and spasms[END_REF]. In contrast to botulinum neurotoxins (BoNT), which show wide genetic diversity (see companion articles in this issue), only one TeNT type has been reported until now. It is presumed that C. tetani strains exhibit limited genetic variability. A previous study conducted via DNA/DNA hybridization showed that 10 neurotoxigenic and 3 non-toxigenic C. tetani strains were homogeneous [START_REF] Nakamura | Taxonomy of Clostridium tetani and related species[END_REF]. At least two genetic profiles have been evidenced in five C. tetani strains, as tested by pulse field gel electrophoresis [START_REF] Plourde-Owobi | Molecular characterization of Clostridium tetani strains by pulsed-field gel electrophoresis and colony PCR[END_REF]. It was first demonstrated that the tent gene is located on a large size plasmid and that the loss of this plasmid yields non-toxigenic strains [START_REF] Finn | The structural gene for tetanus neurotoxin is on a plasmid[END_REF]. The first complete genome sequence of a toxigenic C. tetani strain (E88) was determined in 2003 [START_REF] Bru ¨ggemann | The genome sequence of Clostridium tetani, the causative agent of tetanus disease[END_REF]. The genome contained a chromosome of approximately 2.8 Mb and the tetanus-toxin-encoding plasmid of 74 kb. One additional C. tetani genome was determined in 2014 [START_REF] Levy | Clostridium tetani osteitis without tetanus[END_REF]. Here we report the genome sequences of three additional C. tetani strains. We performed comparative genome studies in order to define the core and flexible genome content of C. tetani.

Materials and methods

Bacterial strains

The three C. tetani strains used for this study were CN655, 184.08 and strain A. Strain CN655 is a collection strain from Prevot's collection (Institut Pasteur) that originated from the Wellcome Institute (1947). Strain 184.08 was isolated from a suppurative arthritis case subsequent to a traumatism in a 75year-old man in the south of France in 2008. Strain A is a variant of the Harvard strain.

C. tetani strains were grown in TGY broth (pH 7.5) containing trypticase (trypticase-glucose yeast peptone BBL, BD Biosciences; 30 g/L), yeast extract (Bacto Yeast Extract, BD Biosciences; 20 g/L), glucose (5 g/L) and cysteine, HCl (0.5 g/ L) under anaerobic conditions.

DNA extraction and genome sequencing

Genomic DNA from all strains of C. tetani was extracted and purified as previously described [START_REF] Popoff | Taxonomic position of lecithinasenegative strains of Clostridium sordellii[END_REF]. Draft sequences of the genomes were obtained by using Illumina/Solexa GAIIx machines at the Beijing Genomics Institute (BGI) (Shenzhen, China). The assembly of sequence reads was done using SOAPdenovo (version 1.05).

Sequence analysis

Automatic annotation of all contigs was done using RAST (http://rast.nmpdr.org/). To identify homologs in different strains, comparisons of all coding DNA sequences (CDS) were made using a protein sequence-based bidirectional BLAST approach (blastP version 2.2.18). Sequence homologies were only mentioned in this study for proteins with an amino acid identity of >25% and an overlap of the query and subject sequence of >75%. Genome representations were created by a DNA plotter (Sanger Institute).

Sequence accession

GenBank accession numbers for the draft genome sequences are JSWC00000000 for strain CN655, JSWD00000000 for strain 184.08 and JWIX00000000 for strain A.

Results and discussion

Genomes of neurotoxigenic clostridia have been sequenced. Whereas many genome sequences of C. botulinum have been determined in the last decade (over 10 completed genomes and more than 30 draft genomes), only two C. tetani genomes are available thus far. The first genome sequence of C. tetani was determined in 2003 [START_REF] Bru ¨ggemann | The genome sequence of Clostridium tetani, the causative agent of tetanus disease[END_REF]. Recently, the genome of strain 12124569 was sequenced, an isolate that caused a rare case of chronic osteitis [START_REF] Levy | Clostridium tetani osteitis without tetanus[END_REF]. We decided to obtain more information on the genomic diversity of C. tetani by sequencing three additional genomes, those of strains CN655, 184.08 and A. Strain 184.08 was isolated from a suppurative arthritis of a 75-year-old man in the south of France in 2008. CN655 is a collection strain from Prevot's collection (Institut Pasteur) that originated from the Wellcome Institute (1947). Strain A is a variant of a Harvard strain. All these strains produce the tetanus toxin as tested in a mouse bioassay (data not shown).

Tetanus-toxin-encoding plasmids

All five strains possess a tetanus-toxin-encoding plasmid. G þ C content, size and number of coding sequences (CDS) vary slightly, with the exception of strain 12124569 that carries a considerably smaller plasmid (Table 1). Sequence comparison revealed that plasmids of strains E88, CN665 and A (hereafter named pE88, pCN665 and pA) are highly similar to each other (Fig. 1); in contrast, plasmids of strains 184.08 and 12124569 (named p184.08 and p12124569, respectively) strongly differ (Fig. 2). Bidirectional blast (BIBlast) comparison revealed the extent of sequence divergence (Supplemental Table 1). Several genes of pE88, pCN655 and pA are not present on p184.08 and p12124569, including a twocomponent system and an ABC transport system. In contrast, p184.08 contains additional genes, including a type III restriction-modification system, a phage integrase and a putative bacitracin resistance transport system (Supplemental table 1B).

The sequences of TeNT and collagenase ColT are highly conserved: TeNT of pE88, pCN655 and pA are 100% identical at the amino acid level and are 98% identical (1294/1315) to TeNT from strains 184.08 and 12124569. Similarly, ColT of pE88, pCN655 and pA are 100% identical and they are 97% (963/991) identical to ColT of p184.08. Other highly conserved proteins encoded by all five plasmids are the two putative RNA polymerase sigma factors (CTP10 and CTP11 in pE88; 94% and 99% protein identity, respectively), an ABC transport system (CTP30/CTP31) and two putative aminopeptidases (CTP20 and CTP56 in pE88; 97% and 98% protein identity, respectively).

Chromosome of C. tetani

The chromosomes of the five C. tetani strains have a similar size (ca. 2.8 Mb) and G þ C content (28.6e28.8%) (Table 2). They are highly conserved: the core genome comprises 85% and 83% of all CDS in strains E88 and 184.08, respectively, when applying a cutoff of 30% protein identity (Fig. 3). Even with a cutoff of 50% protein identity, 84% and 81.5% of all CDS of strains E88 and 184.08, respectively, belong to the core genome (Supplemental Table 2).

The majority of strain-specific functions are encoded in prophage regions. For example, the chromosome of strain E88 contains three strain-specific prophage regions (Supplemental Table 2). Other strain-specific functions include transposases (in particular, in the genome of E88), cell-surface/adhesion/Slayer proteins, restriction-modification-systems, ABC transport systems and efflux pumps, two-component systems, resistance determinants and defense/repair functions as well as specific metabolic functions. For instance, strains 184.08 and 12124569 contain peptide methionine sulfoxide reductase MsrB. Strain E88 specifically contains a plasmid-like locus comprising an adenine-specific methyltransferase and a toxinantitoxin system. Regarding metabolic functions, for example, strain 184.08 possesses the genes for an allophanate hydrolase that converts allophanate to ammonium and carbon dioxide, an essential system for using urea as a nitrogen source.

Another difference between the strains concerns the clustered regularly interspaced short palindromic repeats (CRISPR)/Cas loci and the diversity of CRISPR spacers. The CRISPR/Cas system constitutes an adaptive immune system that protects the bacterium from bacteriophage infections and from mobile genetic elements. Strain E88 contains two loci: one locus (CTC1142-CTC1148) is composed of the cas genes cas6, csx8 (cas8a1), devR (cas7), cas5, cas3, cas4 and cas1 (type I-A-like system) with a repeat consensus 5 0 -GTAT-TAGTAGCACCATATTGGAATGTAAAT-3 0 and the other one (CTC1463-CTC1469) is composed of cas6, csh1 (cas8b), csh2 (cas7), cas5, cas3, cas4, cas1 and cas2 (type IeB system) with a repeat consensus 5 0 -ATTTAAATA-CAACTCTTGTTATTGTTCAAC-3 0 . Most of the spacers (66 in strain E88) show no similarity to sequences deposited in databases. Some spacers match strain-specific prophage regions in other C. tetani and/or C. botulinum genomes. Interestingly, a number of spacers also match regions within the same genome: for example, spacers in E88 match E88 genomic regions encoding stage IV sporulation protein (CTC1140), sporulation sigma-E-factor-processing peptidase (CTC1125), vancomycin b-type resistance protein VanW (CTC0360) and S-layer protein/N-acetylmuramoyl-L-alanine amidase (CTC0462). This indicates that such regions are located on mobile genetic elements that have been acquired by horizontal gene transfer events. The type IA-like CRISPR/Cas system in strain 12124569 is associated with a different repeat (5 0 -AATTTACATTCCAATATGGTGCTACTAATAC-3 0 ) and contains 22 spacers. Some of them match prophage-like genomic regions in C. tetani and C. botulinum. The genome of strain 184.08 does not contain any CRISPR/Cas region. Interestingly, this genome also contains more prophage-like genomic regions that account for its larger genome size

Fig. 1 .

 1 Fig. 1. A BiBlast-based comparison of five tetanus-toxin-encoding plasmids of C. tetani with pE88 as reference plasmid. The blue outer rings depict the CDS of pE88 (plus and minus strand); the four inner rings depict bidirectional Blast results regarding pCN655, pA, p12124569 and p184.08, respectively. Red: high similarity; orange: medium similarity; yellow: low similarity; gray: absence. The highly conserved tetanus toxin and collagenase genes are depicted in dark green and cyan, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) two inner rings depict the G þ C content and the G þ C skew.

Fig. 2 .

 2 Fig.2. A BiBlast-based comparison of five tetanus-toxin-encoding plasmids of C. tetani with p184.08 as reference plasmid. The green outer rings depict the CDS of p184.08 (plus and minus strand); the four inner rings depict BiBlast results regarding p12124569, pE88, pCN655 and pA, respectively. Red: high similarity; orange: medium similarity; yellow: low similarity; gray: absence. The highly conserved tetanus toxin and collagenase genes are depicted in dark green and cyan, respectively. The two inner rings depict the G þ C content and the G þ C skew. Please note that the map of p184.08 was reconstructed from 9 contigs aligned to pE88; thus, sequence gaps in p184.08 exist. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1

 1 Features of tetanus-toxin-encoding plasmids.

	Plasmid name	Size	G þ C (%)	CDS b
	pE88	74.0 kb	24.5	98
	p184.08	ca. 73.0 kb a	24.9	93
	pA	ca. 71.7 kb a	24.5	90
	pCN655	ca. 69.6 kb a	24.3	90
	p12124569	58.4 kb	24.0	77

a Multiple contigs. b Based on RAST annotation (http://rast.nmpdr.org/).

Table 2

 2 Features of 5 chromosomes of C. tetani strains.

		Size	Contigs	G þ C (%)	CDS a
	E88	2799251	1	28.7	2736
	CN655	ca. 2780622	111	28.6	2734
	A	ca. 2754351	90	28.6	2704
	12,124,569	2807472	1	28.8	2809
	184.08	ca. 2841472	135	28.6	2826

a Based on RAST annotation (http://rast.nmpdr.org/).

H. Bru¨ggemann et al. / Research in Microbiology 166 (2015) 326e331

compared to the other four C. tetani genomes (Table 2; Supplemental Table 2B). This indicates that the lack of CRISPR/Cas systems renders strain 184.08 more vulnerable to phage infections. The divergence of the CRISPR/Cas regions of C. tetani reflects the distinct bacteriophage infection history of each individual C. tetani strain.

In conclusion, genome comparison of C. tetani highlights flexible, non-core genome content. Larger islands within individual genomes often constitute prophage-like regions. Smaller genomic regions of divergence encode a set of fitness functions, increase competitiveness and protection and guarantee survival under stress conditions and attack, e.g. against bacteriophage infections.

The tetanus-toxin-encoding plasmids of the five isolates have several differences, underlying their highly flexible nature compared to the chromosome. In particular, plasmids of strains 184.08 and 12124569 differ significantly from pE88, the plasmid from the vaccine-producing strains E88 and A.

Conflict of interest

The authors have no conflict of interest.

Appendix A. Supplementary data

Supplementary data related to this article can be found at http://dx.doi.org/10.1016/j.resmic.2015.01.002.