Brownian cluster dynamics with short range patchy interactions. Its application to polymers and step-growth polymerization

We present a novel simulation technique derived from Brownian cluster dynamics used so far to study the isotropic colloidal aggregation. It now implements the classical Kern-Frenkel potential to describe patchy interactions between particles. This technique gives access to static properties, dynamics and kinetics of the system, even far from the equilibrium. Particle thermal motions are modeled using billions of independent small random translations and rotations, constrained by the excluded volume and the connectivity. This algorithm, applied to a single polymer chain leads to correct static and dynamic properties, in the framework where hydrodynamic interactions are ignored. By varying patch angles, various chain flexibilities can be obtained. We have used this new algorithm to model step-growth polymerization under various solvent qualities. The polymerization reaction is modeled by an irreversible aggregation between patches while an isotropic finite square-well potential is superimposed to mimic the solvent quality. In bad solvent conditions, a competition between a phase separation (due to the isotropic interaction) and polymerization (due to patches) occurs. Surprisingly, an arrested network with a very peculiar structure appears. It is made of strands and nodes. Strands gather few stretched chains that dip into entangled globular nodes. These nodes act as reticulation points between the strands. The system is kinetically driven and we observe a trapped arrested structure. That demonstrates one of the strengths of this new simulation technique. It can give valuable insights about mechanisms that could be involved in the formation of stranded gels.


I. INTRODUCTION
The structure and dynamics of a wide range of complex liquids is determined by the aggregation of small particles in solution such as colloids [1][2][3][4], proteins [5][6][7], micelles [8,9] or oil droplets [10]. Depending on the concentration, the range and strength of the attraction, stable cluster dispersions, transient gels, glasses, phase separated systems can be formed (see for example recent reviews [11,12]). In order to better understand these processes computer simulations have been done on model systems. They are important tools in colloidal physics giving relevant insights about the structure and dynamical properties, even far from the equilibrium.
Systems made of attractive particles in solution can be considered from two different viewpoints. Either we speak of aggregates, which form or breakup continuously, or we speak of concentration fluctuations of individual particles. Which of the two approaches appears most natural depends on the system, e.g., in the case of short range attractions one might tend to speak of reversible aggregation, while in the case of long range attractions it may appear more appropriate to speak of concentration fluctuations. Of course, both approaches are strictly equivalent and it is a matter of semantics whether we consider two neighboring particles whose movement is temporarily correlated as belonging to a single transient cluster or as two different particles whose position is correlated by the influence of an attractive potential.
Classical simulations are based on the particle viewpoint. Newtonian dynamic simulations are numerically solving the equations of motion for a set of particles. This algorithm was first proposed for a simple hard sphere fluid in 1957 [13]. It was then extended to other potential like square-well and Lennard-Jones. Most of the time the solvent is ignored and particles have linear trajectories between collisions. Taking the solvent into consideration by representing each of its molecules would be far too expensive in computation time and memory. There are often billions of solvent molecules per particle. To mimic the particle Brownian motion resulting from the random collisions with the solvent molecules, a random force and a friction term are introduced in the equations of motion; leading to either dissipative particle dynamics (DPD) if hydrodynamic interactions are taken into account or Brownian dynamics (BD) if not (see for example [14,15]). In the case of DPD the friction and the random force felt by a given particle depend on the position of the others. Nowadays BD techniques can deal with a set of 10 4 particles integrating the equations of the motion over a physical time up to few seconds for micrometric colloidal particles in water at 20°C. This approach is very useful to investigate systems that rapidly reach their equilibrium. Static and dynamical quantities can be computed. For colloids, this is generally the case in the one phase domain of the phase diagram (i.e. for weak attractions, low volume fractions of colloids). But for stronger attractions, when phase separation occurs, the time needed to reach the thermal equilibrium is out of the accessible range using molecular dynamic approaches. In that case Monte-Carlo (MC) simulations based on the Metropolis algorithm [16] are employed to determine the equilibrium state in particular to determine the phase coexistence line [17,18].
In the past years, we have developed a novel simulation technique called Brownian Cluster Dynamics (BCD), based on the cluster approach. It can handle up to 10 6 particles during several hours of physical time for micrometric colloidal particles in water at 20°C. It is based on the algorithm originally developed in 1983 by Meakin [19] and Kolb et al. [20] to mimic the irreversible aggregation of Brownian spherical particles which leads to the formation of fractal aggregates. We have extended it to reversible association by allowing already built clusters to break and reform [21,22]. The algorithm does not involve any resolution of the equations of motion and is only based on a probabilistic approach (see below). It consists of chaining a huge number of times the three following steps: (i) clusters are randomly formed among interacting pairs of particles; (ii) particles undergo random small translational displacements maintaining cluster integrity. This step is very similar to an off-lattice version of the bond fluctuation model [23]; (iii) time is incremented. The algorithm has been described in details elsewhere [24] but will be recalled in the next section to include the patchy interaction.
It gives same dynamics as BD simulations and predicts same phase diagrams and static properties as MC simulations [24][25][26]. Hydrodynamic interactions are not taken into account and hence dynamics is obtained in the free-draining case also called the Rouse limit [27].
Hard sphere fluids interacting via an isotropic short range potential are part of self-assembled materials but they often lead to polydisperse and irregular structures.
More precise control over their association to produce predictable structures at the length scale of several particles is an important challenge for novel materials. In comparison, most biological particles assemble into highly monodisperse and precise structures due to the presence of localized and specific interactions at their surface.
Following a "biomimetic" approach, an idea was to model those interactions by decorating the surface of the particles with various sticky patches [28] conferring them an anisotropic potential. Nowadays, patchy particles are the subject of growing interest in the self assembled material community and various ingenious particle synthesis techniques have been developed (see for example the review by Pawar and Kretzschmar [29] and recent articles [30][31][32]). Concomitantly, numerical models have been tested and the influence of that kind of anisotropy on the system behavior has been studied (see for example [28,[33][34][35][36][37][38][39][40]). Sciortino's group has studied the reversible self assembly of particles with two monovalent patches into polymer chains [41,42]. Their simulation results were remarkably well described by the Wertheim theory [43]. Very recently Marshall and Chapman also developed a new theoretical model to describe self assembling mixtures of single [44] or double [45] patch colloids with spherically symmetric ones.
Here, we would like to present an extended version of the BCD algorithm which takes into account patchy interactions localized on the surface of the particle. In addition to the translational Brownian motion, we have to mimic the rotational one in order to relax both particle positions and patch orientations. The potential we use is similar to a square-well limited to some angular sectors introduced for the first time by Jackson et al. [46] and later studied by Kern and Frenkel [33] and Sciortino's group [36,42,47] using computer simulations. To interact, pairs of particles have to be both in range and correctly oriented.
Step (ii) of BCD is modified to take into account random small rotations of particles in addition to translational displacements while still keeping bond lengths and orientations in their tolerance domain (cluster integrity). In near future we aim to employ this new algorithm, hereafter denoted patchy Brownian cluster dynamics (PBCD), to study many relevant experimental systems combining strong or irreversible patchy attractions with a superimposed weak attractive isotropic potential (colloidal aggregation of cementitious materials [48], formation of polymeric proteins like gluten [49] or spider silk [50,51], aggregation of Janus particles [52,53]). In the present paper we centre our attention in linear polymer chains. On the one hand, the linear architecture of polymers represents one of the most stringent scenarios for illustrating the ability of our computational scheme to capture the underlying physics of anisotropic clusters. On the other hand polymers provide a particularly convenient case of study because many theoretical predictions, simulations and experiments have been made on the subject.
Firstly, using our model, we will revisit static and dynamic properties of linear chains with or without excluded volume effects. In particular we will illustrate how the effects of the chain flexibility (or the lack thereof, stiffness) can be tuned in the algorithm by playing with the cone angle of the patchy interactions, and how the PBCD model can be employed to get insights on the equilibrium configuration of semi-flexible polymers, along with their collective relaxation properties.
Next we will apply the PBCD algorithm to describe the step-growth polymerization; i.e. a type of chemical reactions in which bi-functional monomers react, irreversibly, to form first dimers, then trimers and so on… and eventually long polymer chains. As many naturally occurring and synthetic polymers are produced by stepgrowth polymerization, e.g. polyesters, polyamides, polyurethanes, etc [54], it is clear that the physics involved in these process is of great practical use. Most of the existing theoretical studies rely on population balance models [55,56] which keep track of the concentration of every chain length from the knowledge of the starting populations and the effective constant rates [57][58][59]. In practice, some difficulties appear when assigning values to the effective step-growth chemical rates, as it is still not fully understood how they are influenced by chain-length, diffusion effects and ring formation. In this scenario it is clear that the herein presented PBCD algorithm can serve as a valuable tool for getting a complementary insight, since this new computational scheme naturally account for these mentioned effects. Let us note that in terms of the PBCD algorithm the step-polymerization is merely the irreversible aggregation process of interacting particles with two accessible monovalent patches. Tuning the strength of an additional isotropic potential makes it possible to mimic polymerization in solvent with various thermodynamic qualities. Here, we will present some preliminary results in that sense where polymerization and phase separation are in competition.
The remainder of this paper is organized as follows: Section II describes computational details and applied approximations. Section III presents the results obtained on single polymer chains. Section IV gives some preliminary results about the influence of local flexibility and solvent quality on step-growth polymerization. We will particularly focus on the formation of out of equilibrium arrest networks, very similar to stranded gels obtained with biopolymers (see for example the very recent review by Nicolai and Durand [60] and references therein). Section V summarizes the main findings and gives some conclusions.

A General case
The patchy potential is modeled using a square-well (SW) potential restricted to some given orientations between both particles [46]. Particles i and j are in interaction if their distance, r i,j , is within the range and the vector, r i,j , linking both particle centers, crosses both patch surfaces. Assuming d is the particle diameter,  is the relative range of the SW potential and a patch is delimited by a cone of axis v i and half opening angle  (see figure 1) then the potential, V(r i,j , v i , v j ), between both particles is given by: where , are unit vectors, V 1 is an isotropic SW potential of relative width  and depth u 0 < 0: and V 2 given by: if cos and coŝˆ, , 0 else For a given patch,  i is the angle between r i,j and v i .
The simulation starts at simulation time t sim = 0 with an ensemble of N randomly distributed and oriented spheres in a box of size L box with periodic boundary conditions.
The occupied volume fraction is: . A simulation step is made of the three following procedures that are repeated till we reach the desired physical time. To be in interaction with potential -u0, ri,j has to cross both patch surfaces (with solid angle 2π⋅(1-cosω)) and d < ri,j < d⋅(1+ε). vi defines the patch orientation and γi is the angle between ri,j and vi.

Cluster formation procedure
Spheres, also named monomers, are considered to be in contact when they are in interaction condition (with correct range and orientation), we call Nc the number of such contacts. In the so-called cluster formation step, monomers in contact are bound with probability P. Alternatively, bonds are formed with probability α and broken with probability β, so that the P = α /(α + β). In the latter case, one can vary the kinetics of the aggregation from diffusion limited (α = 1) to reaction limited (α → 0) with the same P that P = Nb/Nc, we can express P in terms of u0 [22]:

Movement procedure
The cluster construction procedure is followed by a movement procedure, where 2N times a sphere is randomly selected and an attempt is made to move it either by an elementary random rotation or translation (with 50% probability each). After 2N trials, each sphere has been in average entitled to a rotation step and a translation step in an

Physical Time
After a cluster construction and movement procedure the simulation time (t sim ) is incremented by one. Calling t the physical time, relationships between t, t sim , s T and s R are obtained considering a free diffusing sphere with the following translational ( ) and rotational ( ) diffusion coefficients [61]: with k B the Boltzmann constant, T the absolute temperature and  the solvent viscosity.
For a free particle mean square displacement (MSD), R 2 (t), is given by [24]: At short time ( R 1 1 D t   ), the orientation vector v performs a two-dimensional random walk on the surface of the sphere and we have [61]: The characteristic time, t 0 , is defined as the time needed for the centre of mass of a free sphere to diffuse a distance equal to its squared diameter: t 0 =d 2 /(6D 1 T ). Using t 0 , equation (7) gives the following relationship between the physical time, t, and the simulation time, t sim : The rotational diffusion coefficient of a sphere is simply given by: Taking and combining equations (5), (6), (7) and (8) we obtain a simple relationship between s T and s R : Equation (9) tells us the simulation time needed to reach a given physical time is inversely proportional to the square of the Brownian step size. Too large steps give too many rejections and slow down artificially dynamics but too small ones restrain the physical time accessible in a reasonable simulation time. A good compromise has to be found between the duration of the simulation and a correct dynamics as we will see in part II.C.
In this simulation study the well-width will be set to 10% of the diameter of a sphere ( = 0.1) and  will be given in radian.

B. Modeling step-polymerization and polymer chains
As a test case to our model we have done simulation on step-growth polymerization. In this case, two patches are located in opposite directions at the surface of the monomers. As polymerization is a chemical reaction involving covalent bonds the chemical groups form bonds which can be considered irreversible, so u 0 is set to infinite and bonds are limited to one per patch. In the cluster construction step, we use  = 1 and  = 0 instead of P = 1 to ensure that a previously formed bond is not broken at the In order to study in detail static and dynamic properties of a polymer chain, single chains with various m and  have been generated in the box. This was realized by randomly choosing positions for the consecutive monomers within bond constraints with or without excluded volume. The next bond vector was chosen with a uniform probability density within the interaction volume (a conic shell). In order to generate unbiased self avoiding cases when excluded volume effects are considered, when a random position leads to an overlap the entire chain is rejected and restarted from the beginning. Obviously this limits the maximum size of chains that can be generated with excluded volume effects, especially when  is too high due to an exponential attrition [62].
The measurement of static properties requires only the chain generation step, while dynamics requires many additional movement steps until the chain conformation and position are sufficiently relaxed.
The use of a grid for small patch angles ( < 0.2 using q = 80, see appendix A) leads to some finite size effect biases when generating the next monomer direction.
Limiting our choice to only discrete positions makes the occurrence of the exact forward direction (along v i ) too high and increases artificially the amount of bond angles with  =  exactly. To overcome this problem we had to choose the next direction among a continuous set of values which implied floating point calculations of sines and cosines (while they were pre-calculated on the grid). This method is computationally very expensive and prevents the study of chain dynamics for too small patchs. So static properties were measured down to  = 0.1 using a continuous set of directions while dynamical properties were measured down to  = 0.2 using the grid method.
With our patchy polymer chain model (PPC), the bond length of the polymer can fluctuate between d and d(1+) whose average value is defined by l b  = r i 2  1/2 , then the average contour length is given by L= (m -1)l b . Similarly the angle between a patch and a bond can fluctuate between 0 and  with its average orientation given by cos , and the angle between two consecutive bonds will fluctuate from 0 to 2 with its average value given by cos . For polymer chains without excluded volume constraints (called ideal chains) all these averages can easily be obtained by considering a uniform distribution of bonds within conic shells. For the completion of this discussion, we are giving the formula of l b , cos  and cos  in the appendix B. These local averages will be noted with a star in exponent to indicate ideal quantities.
As given in figure 2, we define G as the chain center of mass and R e as the end-toend vector. R cm 2 (t) is the MSD of the polymer chain center of mass of and R 2 (t) the MSD of an average monomer in the chain. In this work the relaxation of the chain orientation is followed by monitoring the temporal evolution of the normalized correlation function of the end-to-end vector, defined as: Monitoring the MSD of the center of mass of the polymer chain gives access to the translational diffusion coefficient of the chain, D m T : In all cases, angular brackets denote averages over all configurations and all possible evolution of the chain. They are obtained by generating at least 10 5 independent chains for static properties and at least 10 3 independent configurations and temporal realizations for dynamics.  (15) and with  the Flory exponent. = 1/2 for ideal chains and a mean field approach gives   3/5 for self-avoiding ones. A more sophisticated derivation leads to  = 0.588 for self avoiding chains (see [66] for more details). p is the mode index and only odd modes contribute to the relaxation of C(t). The first mode (p = 1) has the longest relaxation time and is called  max . The subsequent relaxation modes (p = 3, 5…) have smaller intensities and relaxation times by a factor 1 2 p   as written in equations (15) and (16).  (14), (15) and (16).  Table 1 gives the expected relative contribution of each odd mode, A( p ), together with the ratio  max / p for both values of . We clearly see that the main contribution arises from  max that is responsible for more than 80% of the relaxation of C(t). The 2 nd relaxation time (p = 3) is around ten times smaller than  max with a contribution less than 10%. This means that for t >  max , C(t) behave mainly as a single exponential: The rotational diffusion coefficient of a polymer chain, D m R , is calculated using the longest relaxation time  max : For rigid objects, C(t) relaxes as a single exponential with  max  m 3 [63].
The distribution of relaxation times can be obtained by a regularized Inverse Laplace Transform (ILT) of C(t) using a constrained regularization calculation algorithm called REPES [67]. ILT is notoriously mathematically ill-conditioned [68] and is very sensitive to noise in the data [69]. But considering the expected shape of the distribution (the first mode is the main contribution and is well separated from the third one) this technique is sufficiently accurate to extract a correct value for  max if C(t) is relaxed around 1%.

C. Finite size effects
As we discussed previously, the Brownian step size, s R or s T (both being related by equation (11)) has to be as small as possible to give correct dynamics. But a too small value of the step size slows down the calculation by increasing the number of simulation steps needed to reach the same physical time (see equation (9)). In figure 4(a), we have to the chain have to respect two constraints in order to move. They should not break any bond nor should overlap. If such an event occurs their motion is rejected but the simulation time is still incremented. As the step size increases, monomers face more and more potential rejections and dynamics is artificially slowed down. In the range of step sizes used, we observe a linear dependency of this phenomenon and that the diffusion coefficient converges toward a finite value at s T = 0, which we define to be the real translational diffusion coefficient of the chain. A similar behavior is observed for rotational diffusion (see figure 4(b)). In the following, all results given for dynamics of a single chain have been extrapolated to zero step size. But using such an extrapolation to study the complete polymerization reaction would be computationally too expensive and we have to define what we consider to be an acceptable step size. We arbitrarily state that to be acceptable apparent dynamics should be within 10% of the real one. We already had to face that problem before in the case of BCD, but with isotropic potential and no rotational motion. We found that to be acceptable, the translational step size had to be small enough compared to any characteristic length scale in the system (see [25,70] for details). This means the step size has to be small compared to the well width but also compared to the average distance to the first neighbor () when the simulation starts with a random distribution of hard spheres. Introducing the rotational motion adds another characteristic length that is the size of the patch which delimits the accessible positions to the tip of the orientation vector. Combining all these criteria with the constraint to be within 10% of real dynamics leads to the following condition: s R /d < /10 and s T /d < /5 and s T < (-1)/3.
It is well known that introducing excluded volume interactions lead to the swelling of the chain which leads to increasing characteristic sizes of the polymer chain (end to end distance, radius of gyration…). This effect is not a simple proportionality factor but also changes the Flory exponent involved in the description of the asymptotic scaling of R e 2  with the number of monomer in the chain [72]: Despite average local quantities are very little influenced by that effect when  is small, the bond correlation function, cos (n) is strongly modified. cos (n) is the average cosine of the bond angle between bond r i and r i+n . For an ideal chain, it is a single relaxing exponential with a relaxation length l p given in appendix C. But for large self avoiding chains, cos (n) must scale as: so the correct scaling from equation (19) is recovered [73].
In figure 9 we have plotted cos (n) for a self avoiding chain with m = 210 4 and  = 0.2. We clearly see that up to a given n, the relaxation is the same as the one of an ideal chain with a persistent length calculated using equation 35 and 45 (see appendix C): Then a crossover takes place and equation (20) is recovered. For a chain with finite length, as n approaches its maximum value (n = m-2) we observe a deviation from the power law and a final cut-off occurs. In figure 10 we see that the complete behavior of the bond correlation can be remarkably adjusted all along using the following phenomenological equation: where A, n 1 and n 2 are fitting parameters. Figure 10 shows the fit made on a self avoiding chain with m = 810 3 and  = 0.4 using a nonlinear regression with the Marquardt-Levenberg algorithm (see for example [74]). For a given , values obtained for n 1 are independent on m as soon as m >> n 1 , see figure 11(a). Figure 11(b) shows that n 1 is indeed equal to l p /l b  where l p has been calculated using cos: cannot be used to characterize the local flexibility for self avoiding chains (see [75] for details). In the following, we have fixed n 1 = l p /l b  in the fitting procedure. This leads to more accurate determinations for A and n 2 . Figure 12    Finally, figure 13 shows l p , calculated with equation (23), as a function of  for PPC with excluded volume interactions together with the theoretical prediction for ideal chains. In the range 0.1    , the persistent length of the self-avoiding PPC model varies over 2 decades. In the limit of small , l p is not influenced anymore by excluded volume effects and scales like: FIG. 13. Evolution of the ratio l p /l b  as a function of  for excluded volume chains. The solid line represents the predicted ideal behavior (see equation (21)). Error bars are smaller than the symbol size

D. Deviation from ideality, influence of the local flexibility on the thermal blob size
In order to study the deviation from ideality, R e 2  has been normalized by its ideal value, R e 2 * (see appendix C), and plotted for various m and . Figure 14  Using scaling arguments for semi flexible chains with excluded volume interactions, Schaefer et al. [77] have proposed an explanation to describe the length scale, l T , at which that cross-over takes place. They use the thermal blob concept and assume the chain is made of blobs, with characteristic squared length scale  T 2 , freely jointed in a self avoiding manner. Inside each blob, containing m T monomers, the chain is considered as ideal with persistent length, l p , and bond length l b . l T is simply the contour length inside a blob, given by: l T = (m T -1)l b . On one hand, considering that kind of structure leads to the following relationship: On the other hand, taking a Flory approach to calculate the free energy of the semi flexible chain and minimizing it gives for very long chains in good solvent (high temperature limit see equation (7) in [77]): Taking  = 3/5 and combining equations (25) and (26) Schaefer et al. [77] predict: Unfortunately, for the smallest , we cannot generate sufficiently large chains such that the swollen behavior in unambiguously recovered. Chains are still in the crossover regime and the limiting scaling law, with exponent 2-1, is not yet valid (see figure   14). However, assuming that for  sufficiently small the ratio R e 2 /R e 2  * is only a function of L/l T we have built a master curve by horizontally shifting curves on a log scale for various  on the top of a reference curve ( = 1.00), trying only to superimpose values at larger L. The shift factor enable the calculation of l T for each  assuming l T at  = 1.00 being arbitrary defined as the intersection of the limiting power law with the horizontal line R e 2 /R e 2  * = 1 (l T /d  20, see figure 14(a)). Results are presented in figure 15 where data of figure 14 are plotted but as a function of L/l T . Except for the highest , data superimposes reasonably well and the cross-over can be described with the following phenomenological equation: with a = 0.3. This cross-over is large and takes place over about four decades (10 -2 << L/l T << 10 2 ). Figure 16 shows the evolution of l T /d as a function of l p /d. The theoretical behavior predicted by Schaeffer et al. [77] is well recovered (see equation (27)).  The PPC model appears to be in very good agreement with theoretical predictions for static properties of self avoiding polymer chains. It is an off-lattice model that enables a very convenient tuning of the local chain rigidity by playing on the patch angle. We did study other static properties like the local persistent length, l p (k), defined as the projection of the local bond vector k on the end-to-end vector. This quantity has been introduced by Yamakawa [78] and theoretically investigated by Schäfer and Elsner [79]: It has received much attention recently [75,80,81]. Our results concerning l p (k) are beyond the scope of that paper and will be published elsewhere.  (13)). This can be done very accurately even after short running times. In   (14)).

OUT OF EQUILIBRIUM ARRESTED STRANDED GELS
This part is dedicated to some preliminary results about the competition between polymerization and phase separation. It demonstrates the efficiency of our algorithm to investigate arrested out of equilibrium structures.
Here, we introduce a secondary isotropic SW potential coupled to the patches.
Equilibrium properties of such a system have been investigated by Liu et al. [82] using MC simulations and a thermodynamics perturbation theory approach. They limit their study to one possible bond per patch using geometrical constraints to be consistent with the theoretical assumptions [83]. They focused on the influence of the number of patches and their coverage on the resulting phase diagram and showed that the second virial coefficient is a scaling parameter for a generalized law of corresponding states [36,84]. In this study we concentrate on the competition between an irreversible patchy t/ max 0.001 0.01 0. polymerization (highly directional) coupled with an isotropic reversible interaction that mimics the quality of the solvent. For simplicity, the range of the isotropic SW is the same as that for patches ( = 0.1). In a bad solvent, patch free particles will have a tendency to self associate and even phase separate. The resulting SW fluid has already been extensively studied using BCD [24][25][26]. Its adhesiveness can be characterized by the attractive part (B att ) of its reduced second virial coefficient (B 2 ). We speak about a reduced quantity that has been divided by the particle volume to obtain dimensionless numbers. For a SW fluid we have: B att combines in a single parameter the range and the energy of the interaction. It has been shown that in the limit of small ranges ( << 1), the phase diagram is no more influenced by the precise shape of the potential and a unique phase diagram is obtained in the plane (B att , ) [84]. Previously we have shown that systems with  = 0.1 can be considered in that limit [25,85]. When B att = 0, the system is in good solvent conditions but as B att increases, its thermodynamic quality decreases and favors isotropic reversible aggregation in the solution. When B att = 4, the hard core repulsion is balanced by the attractive part of the potential and B 2 = 0. This corresponds to the Boyle temperature of the fluid. To some point, it can be considered to be close to theta conditions for the polymer chain [86]. Figure 25 (adapted from [26]) shows the phase diagram of a SW fluid. To trigger a spontaneous crystallization it is necessary to work below the metastable liquid-liquid binodal. Phase separation above that line is nearly impossible without introducing crystalline seeds in the solution. Snapshots in figure 25 compare a system with  =10% at the same time, above (B att = 4) and below (B att = 12) the liquid-crystal binodal. In the first case the system is homogeneous while in the second one, we observe crystallization that takes place in dense liquid droplets. We introduce now patchy irreversible aggregation to take into account the stepgrowth polymerization process in both solvent conditions. An infinite patchy potential is simply added to the finite isotropic one. For some convenience we have used patches with angle  = 0.8179, so that polymerization kinetics is not too slow and takes place within similar time scales as the phase separation. In that case, patches cover a total amount of 32% of the particle surface. In good solvent (B att = 0, self avoiding chains) or close to the theta condition (B att = 4, ideal chains), isolated chains built with that cone angle would be rather flexible with a persistent length around 3 bonds (see figure 14).  To verify whether it is arrested or not, we have first monitored the degree of reaction of the polymerization, x, which is simply the fraction of reacted patches.
Plotting 1-x or m n as a function of the time on a log-log scale reveals clearly the difference of kinetics between B att = 12 or 4 (see figure 28). The reaction appears to have stopped for B att = 12 while it continues to slowly evolve at B att = 4. The system in figure   27 is "chemically" arrested. We have also monitored dynamics in both cases by plotting the MSD of an average monomer during the polymerization reaction for two different waiting times (t w ). Figure 29 shows the MSD of an average monomer for both solvent conditions with t w /t 0 = 0 and t w /t 0 = 5378. For t w /t 0 = 0 monomers exhibit free diffusing dynamics at short times which becomes sub-diffusive as the reaction proceeds. At large times, for B att = 12 a plateau is reached while for B att = 4 dynamics is more compatible with the one of a Rouse polymer chain at a time scale smaller than the Rouse time (see for example [72]). Restarting to monitor the MSD at t w /t 0 = 5378 we also notice a clear difference of behavior. At B att = 12 a monomer endures in average a displacement of its diameter over a period of time equal to around 5000t 0 . At B att = 4 monomer dynamics is much faster and seems to converge toward the one of a monomer in a Rouse chain at a time scale smaller than the Rouse time. For B att = 12, the reduction of dynamics at larger time scales supports the argument of a dynamically arrested system. Only "breathing motions" of the stranded gel remain visible (see both movies in supplementary materials [87], they illustrate both kind of dynamics) and the system can reasonably be considered chemically and dynamically arrested.
We have checked the effect of the size of the simulation box using a bigger box and found that the resulting structure does not depend on the system size. Figure 30 shows two snapshots taken at the same time (t/t 0 = 603) for L box /d = 100 (containing 190985 monomers) and L box /d = 37 (containing 9674 monomers) which corresponds to  = 0.1 and at B att = 12. We do not observe appreciable differences in the system and both systems reach the same polymerization degree at the same time. However, the computation time needed for the bigger system is about thirty times bigger and does not allow us to go much further. Figure 31 shows the effect of concentration on the obtained network structure.
When the concentration decreases, the correlation length of the system seems to increase and the strands become more stretched and contain a lesser number of chains.
Whether there exists a critical concentration below which the network structure is not

V. CONCLUSION
Finally to conclude we have seen that the PBCD algorithm applied to a single linear polymer chain in good solvent gives expected static and dynamical properties.
Billions of small local random motions (translations and rotations), constrained by excluded volume interactions and the connectivity of the chain, lead to the correct collective dynamics in the framework where hydrodynamic interactions are ignored. It was also an opportunity, using the PPC model, to revisit the influence of the local flexibility on the bond correlation function and on the transition from ideal chain to swollen chain for self avoiding semi flexible polymers.
We would like to note here that PBCD is not specially designed to study dynamics of polymer chains but rather to focus on complex aggregation reactions involving both isotropic and directional potentials. Neither it is optimized to determine accurately equilibrium states if they exist and a Metropolis approach has to be preferred for that purpose. However PBCD is particularly suited to monitor the temporal evolution of particulate systems far from their equilibrium (like phase separating system with strong interactions) or out of equilibrium when they are kinetically driven by chemical reactions involving the formation of irreversible bonds. The example given in Sec. IV illustrates the ability of our algorithm to give relevant insights where classical tools (molecular dynamics or theoretical approaches) are inoperative. It enables us to study the competition between polymerization and phase separation that leads to some very interesting out of equilibrium structures like stranded gels. These structures can be viewed as arrested microphase separation resulting from a balanced competition between both effects. The solidity of the resulting 3d network results from a combination of irreversible bonds (by polymerization) and highly cooperative reversible interactions between polymer chains that form gel strands and nodes.
With this model we can study static, kinetic and dynamic properties of irreversible as well as reversible aggregation processes with varying the patch size, the strength of the anisotropic potential adding or not an isotropic potential for competition effects. This could provide more insights into the connection between gels and gas-liquid behavior of colloidal systems. We would like to extend our model to study physical systems like Janus particles but also the formation of calcium-silicate-hydrate (CSH) gels formed during cement hydration [88]. In the later case, it could be interesting to mix an aggregation process with the nucleation and growth of the basic building blocks of CSH gels [89].