Tim Hupfeld

Gaétan Laurens

Samy Merabia

Stephan Barcikowski

Bilal Gökce
email: bilal.goekce@uni-due.de

David Amans
email: david.amans@univ-lyon1.fr

Dynamics of laser-induced cavitation bubbles at a solid-liquid interface in high viscosity and high capillary number regimes

Keywords: Laser ablation in Liquids, oil, gold, yttrium iron garnet, contact angle, bubble shape, 1. Introduction

No unified model is available yet to explain the dynamics of laser-induced cavitation bubbles during laser ablation of solid targets in liquids, when an extremely high capillary number is achieved (>100), i.e. when the viscous forces strongly contribute to the friction. By investigating laser-induced bubbles on gold and yttrium iron garnet targets as a function of the liquid viscosity, using a nanosecond laser and an ultrafast shadowgraph imaging setup, we give a deeper insight into what determines the bubble dynamics. We find that the competition between the viscous forces and the surface tension (capillary number Ca), on the one hand, and the competition between the viscous forces and inertia (Reynolds Number Re), on the other hand, are both key factors. Increasing the viscous forces and hereby Ca up to 100 has an impact on the bubble shape and result in a very pronounced rim, which separates the bubble in a spherical cap driven by inertia and an interlayer. The temporal evolution of the footprint radius of the interlayer can be addressed in the framework of the inertiocapillary regime. For an intermediate viscosity, the thickness of the interlayer is consistent with a boundary layer equation.

Interestingly, our data cannot be interpreted with simplified hydrodynamic (Cox-Voinov) or molecularkinetic theory models, highlighting the originality of the dynamics reported when extremely high capillary numbers are achieved. Upon bubble collapse, spherical persistent microbubbles are created and partly dispersed in water, whereas the high-viscous polyalphaolefines lead to long-standing oblate persistent bubbles sticking to the target´s surface, independent of the ablated target. Overall, liquid´s viscosity determines laser ablation-induced cavitation.

of accessible materials [1][2][3][4][5] . In contrast to chemical routes, LAL works without using surfactants or chemical precursors in a one-step process 6 with a wide variety of liquids [7][8][9][10][11][12][13][14] . Although LAL is scalable 15 , economically feasible 16 and has a high potential for applications in biomedicine 17,18 , catalysis 2,19 , optics [20][21][22] , additive manufacturing 23,24 and nanoparticle-polymer composites 25,26 , there is still a lack of understanding for the basic mechanism involved during nanoparticle synthesis. The dynamics after the impact of a short laser pulse on a target in a liquid environment can be separated into the following stages. An early laser-induced plasma is quickly quenched a few microseconds after the impact, leading to a vapour bubble [27][28][29][30][31][32] containing nanoparticles [33][34][35][36] . The bubble grows and collapses for first within a few hundreds of microseconds and finally leads to the release of the nanoparticles into the liquid 34,35 when the bubble fully collapses and disappears after a few oscillations 26,[37][38][39] . Persistent microbubbles are released during the collapse phase, which may shield a large portion of subsequent laser pulse in particular at high liquid viscosities 40 , and consist of chemical reaction products of the solvent and target, such as volatile carbohydrates 40 , hydrogen, oxygen and hydrogen peroxide 41 .

The dynamics of the cavitation bubble and its influence on the nanoparticle formation has been investigated by many researchers 28,[37][38][39][42][43][44] . Surprisingly, there is a lack of literature investigating cavitation bubbles in highly viscous liquids, although viscous liquids are becoming more and more interesting for LAL applications in oils 45,46 or monomers 47 through offering the possibility for one-step synthesis of colloids. For example, nanoparticles as additives in engine oil are relevant for the automotive industry, where nanoparticles are often used to enhance the tribological properties of the lubricant 46,48 . The one-step synthesis via LAL could replace the process of dispersing chemical synthesised particles in oils which often leads to agglomeration and thus to lower functionality. However, it is still mostly unknown how viscosity and surface tension influence shape, lifetime and size of the laser-generated bubble, and how they influence the nanoparticle growth and agglomeration inside the bubble [START_REF] Tomko | Colloids Surfaces A Physicochem[END_REF]50 . A deeper understanding could offer new ways to ultimately improve the overall LAL process, i.e. repeatability and productivity. A perfect hemispherical bubble is expected in an inviscid flow 51 , and most of the above-mentioned studies on laser ablation in low-viscosity liquids (water, ethanol …) assume a hemispherical shape. Only a few reports deal with the complete description of the bubble shape (rim, interlayer, contact angle). Tomko et al. reported an asymmetric shrinking 35 . Contact angle hysteresis and bell-shaped bubbles in water have been observed 39 . 36 which agrees with a numerical study by Lechner et al. using finite-volume method applied to Navier-Stokes equations for a compressible fluid 51,52 . From studies of cavitation bubbles close to a solid-liquid interface it is further known that the viscosity has a significant influence on the bubble shape 53 . Furthermore, the very fast and thin jet observed in the cavitation-collapse-simulations by Lechner et al. is unlikely to be observed at liquid viscosity 40 times higher than water as the annular inflow is no longer fast enough. At such high viscosities, a jet forms which is much wider and much slower (on the order of 100 m/s) as in water (about 1000 m/s) [personal communication with Christiane Lechner]. Still, there is no model available yet to fully describe the shape of the LAL-induced bubbles for any viscosity. Moreover, laser-induced bubbles close to a surface are characterised by fast moving vapour-liquid interface. By experimentally investigating the bubble dynamics in viscous liquids, we could get access to an unusual high value for the capillary number and address original and unprecedented condition in the literature. In this study, we investigate the relevant parameters (e.g. bubble size, velocity, lifetime, contact angle, contact area) for laser ablation of both gold and yttrium-iron-garnet (YIG) in water and in highly viscous polyalphaolefin (PAO), and we take a deeper look on the origin of the deviations from the semi-spherical shape.

Ibrahimkutty et al. reported an inward jet

Experimental

The experiments were carried out in commercially available Spectrasync PAO6 (density ρ=0.827 g/cm 3 , kinematic viscosity 𝜈=80.8 mm 2 /s @ 293°K) and Spectrasync PAO40 (density ρ=0.850 g/cm 3 , kinematic viscosity 𝜈=764 mm 2 /s @ 293°K) supplied by Exxon Mobile. The viscosity of PAO is influenced by its chain length. Pure water (18.2 MΩ/cm 2 ) was used as a reference fluid (density ρ=0.997 g/cm 3 , kinematic viscosity 𝜈=1.00 mm 2 /s @ 293°K, surface tension 𝛾=72.8 mN/m). The surface tension γ also changes to 29.7 mN/m for PAO6 and 31.5 mN/m for PAO40 (@ 297°K). As ablation targets we have chosen a commercially available gold foils (99.99 %) and yttrium-iron-garnet (YIG, Y3Fe5O12) wafers to experience different wettability. The uncleaned gold foil is assumed more hydrophobic than the hydrophilic oxides. Note that the hydrophilicity does not only depend on the material, but can also be altered by impurities, microstructures, and adsorbed species on the surfaces 54 . The first bubble expansion and its collapse are very fast with respect to the heat conduction characteristic time.

Therefore, the bubble growth is assumed to be adiabatic as it was shown by Lam et al. 28 and we can neglect the heat transfer from the inner bubble to the liquid. As a result, we consider the kinematic viscosity and the surface tension constant during the first bubble expansion and collapse.

For the fast imaging shadowgraphy experiments, a target is placed in a cubic vessel with 6 ml of liquid.

The liquid layer thickness above the target is 6.5 mm. For the sake of comparability, the liquid layer thickness above the target is fixed throughout all experiments. The third harmonic of a pulsed Nd-YAGlaser (355 nm, 5 ns, 9 Hz, 7.2 mJ/pulse) is focused on a spot with a diameter of approximately 260 µm, resulting in a fluence of 13.6 J/cm 2 . In our setup geometry, the liquid´s transmission at 355 nm for a 6.5 mm liquid thickness is 99% for water, 97% for PAO6, and 96% for PAO40, ensuring similar fluence for the three liquids (see supplementary Fig. S1). Ablation is observed at 210 000 frames per second, which corresponds to an integration time of 4.75 µs, by an ultra-fast camera (Phantom v711, Vision Research). On the camera, a Zoom 6000 from Navitar is mounted. The magnification of the optical system is 0.65. According to its performance specification, the resolution limit of the Zoom 6000 optical assembly is then 31 µm, which also corresponds to the image of a single-pixel (20 µm in size divided by the magnification). It leads to an overall resolution of 44 µm (√2 x 31 µm). A pattern of ultra-bright LEDs is used for illumination. To avoid distorting the results by inhomogeneities of the target surface, or through changes in the liquid parameters causing changes in the target illumination, each experiment is performed on a fresh target spot in fresh liquid while the vessel is cleaned between each experiment. The camera is synchronised with a laser beam shutter to ensure that the camera starts imaging right after the first laser shot on a fresh surface. Fig. 1a shows the dynamics of the bubbles in three liquids for a gold target (videos for both the gold and YIG targets are available in the supplementary materials). Starting with the first image 4.75 µs after the laser impact, the bubble starts to expand, reaching its maximum expansion after 100 -150 µs.

Results

Interestingly, the bubbles in PAO6 and PAO40 continue expanding in the vertical direction, even after the expansion in the horizontal direction stopped, highlighting the strength of the friction. The footprint of the bubble in PAO40 stays almost constant during the collapse which could affect the redeposition of ablated material on the target surface (see Fig. 1b). Fig. 1c depicts the vertical bubble radius (height) for the different liquids as a function of time.

Assuming a rotational symmetry of the bubble shape, the maximum volume is deduced from the bubble shape at its maximum vertical radius (Fig. 1d). The average maximum volume for each liquid is shown in Fig. 1e. After the maximum expansion, the shrinking starts, ending with the collapse of the cavitation bubble followed by bubble oscillations (videos of the whole process including oscillation are provided in the supplementary information). The total lifetime of the first bubble is between 190 and 280 µs, depending on the liquid. For each bubble, the lifetime is consistent with the Rayleigh collapse time deduced from the maximum vertical radius (See supplementary Fig. S2). Compared to water, the number of oscillations is smaller in the PAOs. Obviously, the highly viscous PAOs provide higher damping of the bubble oscillation through viscous energy dissipation, resulting in only one oscillation for PAO40. After oscillations, persistent bubbles 55 are formed from the cavitation bubble and their mobility decreases dramatically for high viscosities. The collapse in water releases microbubbles into the liquid that only partly adhere to the target surface. In contrary, LAL in PAO40 results in persistent oblate-shaped bubbles sticking to the target, independent of the ablated target material (see Fig. 2c and the videos in the supplementary information).

There are changes when increasing the viscosity of the liquid, not only in the lifetime and the maximal size of the bubble, but also in the shape and contact angle between the target and the bubble. The physics of spherical and hemispherical bubble oscillation has extensively been studied in the past 56 with focus on bubbles in water. Various studies are also available for cavitation bubbles at different distances from a solid boundary 51,57,58 To discuss the differences in the bubble dynamics in the present case, we have developed a Python code to deduce from the movie the time evolution of the relevant geometrical parameters characterizing the bubble shape. For each image, the bubble shape is fitted using a circular arc for the top of the bubble (spherical shape) and a continuous line for the interlayer (conical shape). The shape of the interlayer is assumed to be conical to measure the apparent contact angle θ. The geometrical quantities deduced from the fit are the height of the bubble R, which remains almost equal to the radius of curvature during the bubble expansion, the footprint radius X, the distance L between the bubble boundary (edge of the interlayer) and the centre of the bubble, the thickness e of the interlayer, and the contact angle θ (see Fig. 2a-b). From the measurement of R, X and L as a function of time, we can compute velocities:

t cl e dR dX dL V = ; V = ; V = dt dt dt (1)
Vt is the velocity of the vapour/liquid interface at the top of the bubble, Vcl is the velocity of the apparent three-phase contact line, and Ve is the velocity of the vapour/liquid interface at the edge of the bubble. account an uncertainty of ±30% on the viscosity, an uncertainty of ±3% on the surface tension 59,60 , and an uncertainty of ±22 µm on the geometrical parameters (R and X). and 3 orders of magnitude, respectively, and is close to unity for the more viscous oil. This trend is confirmed by the evolution of Ca. In water, Ca is less than unity in the first µs before reaching 10 -2 at the maximum bubble size. In PAO6 and PAO40, the same evolution is shifted by 2 and 3 orders of magnitude towards high capillary numbers. Contrary to the case of water, the viscous force cannot be neglected anymore with respect to inertia in polyolefin. The relatively high polyolefin viscosity leads to contact-line friction [START_REF] Bertrand | Colloids Surfaces A Physicochem[END_REF] which competes with inertia and then drives the bubble motion close to the surface.

The main differences in the bubble dynamics between different viscosities can be observed in the behaviour of the interlayer defined in fig 2 . The spreading of the bubble on the solid target can be compared to the advancing and receding sessile droplet [START_REF] Lavi | Colloids Surfaces A Physicochem[END_REF][63][64] or a spreading bubble 65,66 . During the early stage of bubble expansion (first quarter of the bubble lifetime), the time-dependent footprint radius X(t) of the bubble follows the power-law X(𝑡) ∝ 𝑡 𝑛 (see supplementary Fig. S3 andS4). The corresponding exponent n which characterizes the spreading regimes is shown in Table 1 and indicates only a weak dependence of the bubble expansion on the target material. The exponent n is around 0.39 for water. Such value is consistent with regimes where the dominant resistance is liquid inertia.

For a bubble characterized by a constant volume, these regimes lead to theoretical time evolution between t 1/2 and t 1/3 65 depending on the driving force, the gravity or the capillarity, respectively, which counterbalances the inertia. If the bubble size is larger or smaller than the capillary length, this leads to a flattened or a spherical shape, respectively. By considering our bubble shape as hemispherical, an inertiocapillary regime with t 1/3 is in good agreement with the data in Table 1 for water and PAO6.

However, for laser-generated bubbles, the driving force during the early expansion of the bubble is clearly different and corresponds to the bubble inner pressure. Moreover, the bubble volume is not constant. In water, for large Weber and Reynolds numbers, a simplified Rayleigh-Plesset equation applies and reflects the balance between inertia and inner pressure (PB) 28 :

𝜌 (𝑅𝑅 ̈+ 3 2 𝑅 ̇2) = 𝑃 𝐵 (𝑡) -𝑃 𝑙 (5)
where Pl denotes the surrounding liquid pressure (<<PB), and the inner pressure PB follows the isentropic relation 𝑃 𝐵 ∝ 𝐴 𝑅 -3𝛼 with 𝛼 the heat capacity ratio and a constant A 28 . It leads to the balance 𝜌 𝑅 2 𝑡 2 = 𝐴 𝑅 -3𝛼 , and then to the power law coefficient 𝑛 = 2 3𝛼+2

. For a laser-generated bubble, we can expect 𝑛 = 1 3 ⁄ for water (𝛼 water =1.33) and 𝑛 = 2 5 ⁄ for large molecules (𝛼 tends to 1). These values are consistent with the measured ones for water (X(𝑡) and R(𝑡) are identical for water) and PAO6. In contrast, n is significantly smaller for the most viscous PAO40 (n≈0.2) which indicates a change of regime to domination of viscous forces. By increasing the viscosity, the Reynolds number decreases, and the Capillary number drastically increases. Close to the surface, friction has to be included in the resistance forces along with inertia. For a droplet characterized by a constant volume, when inertia is marginal, the balance between the surface tension and the viscous dissipation leading to friction at small scale is described by Tanner's law with t 1/10 . A decrease of the n value with increasing viscosity makes sense, even if Tanner's law assumes again a constant volume, as well as a low capillary number (Ca<<1), which only applies for the bubble's dynamics close to their maximum radius (see figure 3). To our best knowledge, there is no study available in the literature describing a droplet or bubble spreading in a high capillary number regime such the ones reached for high viscous PAOs, i.e. up to 100. However, large Ca numbers are achieved in the context of high-speed coating, up to 1000 for the coating of optical fibres (see 67,68 and references inside). The velocity of the apparent three-phase contact line is drawn as a function of the apparent dynamic contact angle in Fig. 4 for PAO6, PAO40 and water. In the case of water, the expanding bubble is hemispherical, and the apparent contact angle is constant around 90° during most of the bubble lifetime (also see Fig. 1a) due to the inertial regime with a low capillary number (see Fig. 3c). There is only a deviation from 90° during the collapse phase (see Fig. 4) characterized by a sharp acceleration of the bubble. When increasing the liquid viscosity and hereby the capillary number Ca by several orders of magnitude (see Fig. 3c), we observe the appearance of hysteresis between expanding and shrinking bubble for the apparent contact angle between target and bubble. The hysteresis is similar to the behaviour observed for a moving liquid droplet or a dynamic sessile drop at small Ca numbers 63,69,70

. In the limit of the small Ca number, the apparent contact angle is described by the Cox-Voinov hydrodynamic model of dynamic wetting 71 and has been extensively studied 63,69 . But this relation does not apply here. A molecular-kinetic theory has been also developed to describe the velocitydependence of the dynamic contact 58,[START_REF] Bertrand | Colloids Surfaces A Physicochem[END_REF] Blake et al. introduced activation energies of adsorption and desorption from the surface of the molecules of the liquid. The velocity-dependence of the dynamic contact angle is then due to the perturbation of the adsorption equilibrium due to motion of the contact line. Yet again, the model didn't succeed to fit the data for the largest capillary numbers, highlighting the originality of the dynamics reported here.

Fig. 4: Velocity of the apparent contact line as a function of the apparent contact angle plotted for a few laser-induced bubbles on a gold target in PAO6 (green squares), PAO40 (blue triangles) and in water (red circles). Pictures correspond to a laser-generated bubble in PAO6.

The interlayer does not only come with a time-dependent angle, but also with a time-dependent height.

𝑉 𝑒 = 𝑑𝐿 𝑑𝑡 ⁄ is the velocity of the liquid/vapour interface at the edge of the bubble, ν is the kinematic viscosity and L is the distance between the bubble boundary at the edge of the interlayer and the centre of the bubble. The slope of the regression lines in Fig. 5

(b) represent the dimensionless constant

Zb. The slope is almost independent from the target material with Zb,YIG=1.1 and Zb,Gold=0.9. A similar behaviour has been observed by Van Ouwerkerk for bubbles on a hot surface during boiling of the liquid 72,73 . Assuming no change in the shape during the bubble growth, Van Ouwerkerk deduced Zb=0.9

from a self-similar approach and reported experimental values for Zb between 0.3 and 0.8.

Conclusion

Laser-induced cavitation in liquids is a key method applied to understand cavitation near solid boundaries, but also highly relevant for laser synthesis of colloids. Indeed, the cavitation bubble but also the resulting persistent microbubbles may shield subsequent laser pulses, affecting productivity and reproducibility of the laser generation of nanoparticles. While it is intuitive that viscosity plays a role, experimental investigations under variation of both the target and liquid types have been missing.

In this work we showed that an increase of the kinematic viscosity by 2 orders of magnitude compared to water does affect the lifetime, the shape, the size and the damping of the oscillation of laser-induced cavitation bubbles. Bubbles can be divided into two geometrical sections, separated by a rim: A spherical cap with a dynamic driven by inertial forces and an interlayer close to the target where the shape differs from a sphere when the contribution of the viscous forces to the friction drastically increases. Indeed, looking at the bubble footprint dynamics close the target, an inertiocapillary regime is observed for intermediate viscosity and the thickness of the interlayer is in accordance with a boundary layer equation. Moreover, a contact angle hysteresis appears with increasing viscosity. These findings poorly depend on the target material which suggests that the shape and the dynamics of the bubble are driven by the contribution of the viscous forces to the friction, rather than by static hydrophilicity of the target. However, by increasing the viscosity, the behaviour of the contact line cannot be interpreted with simplified hydrodynamic (Cox-Voinov) or molecular-kinetic theory models, highlighting the originality of the dynamics reported when extremely high capillary numbers are achieved. In addition to viscosity effects on the first cavitation of laser-induced bubbles, persistent bubbles show noteworthy phenomenological differences caused by the liquid and appear independent of the target type. These bubbles are spherical and partly dispersed for laser ablation in water, whereas in PAO40 a larger, flat, persistent bubble adheres to the target.

Fig. 1 :

 1 Fig. 1: (a) Shadowgrams of the cavitation bubbles after the impact of a ns-laser pulse on a gold target in water (top), PAO6 (middle) and PAO40 (bottom). Videos of both gold and YIG LAL in water and PAOs are available in the supplementary materials. (b) Zoom on the collapse in PAO40 with a focus on the footprint diameter, as well as the persistent oblate bubble after oscillation. (c) Vertical radius (height)

Fig. 2 :

 2 Fig. 2: (a) Image of a laser-generated bubble on gold in PAO6 showing the relevant geometrical parameters. (b)The red and black curves correspond to the bubble shape returned by the image processing (raw data) for the expanding and receding phases, respectively. Inset: the black curve is the bubble contour obtained from raw data and the blue curve corresponds to its fit. c) Images of the cavitation bubbles at their maximal expansion (left panel) and persistent bubbles after oscillation (right panel) depending on target (gold, YIG) and liquid (water, PAO6, PAO40) combinations.

Fig. 3 :

 3 Fig. 3: (a) Weber number We and (b) Reynolds Re number calculated from the velocity 𝑉 t at the top of the bubble. (c) Capillary number Ca was calculated from the velocity 𝑉 cl of the apparent contact line. The normalized time corresponds to the duration of the first expansion and shrinking of the bubble. No movement is observed (𝑉 cl = 0) for the apparent contact line during the period delimited by the arrows. Error bars are deduced from the propagation of uncertainty (see supplementary) and take into

Fig. 3 (

 3 Fig.3 (a-b) show the Reynolds number Re and the Weber number We computed from the velocity 𝑉 t of the top of the bubble, while figure3 (c)shows the capillary number Ca computed from the velocity Vcl of the apparent contact line.

Fig. 5 (

 5 a) shows the time-dependent height of the interlayer during the expansion for the two targets in PAO6. The error bars are based on the calculated temporal and special resolution of the experimental setup. The interlayer for PAO40 is not addressed because the round shape of the bubble would lead to an arbitrary definition of the parameters L and e. For PAO6, Fig. 5 (a) shows no differences in the height of the boundary layer between a gold and an oxide target. The thickness of the interlayer increases with time and fits the following boundary layer equation with the dimensionless constant Zb: 𝑒 = 𝑍 𝑏 √ νL 𝑉 𝑒

Fig. 5 :

 5 Fig. 5: Measured thickness e of the interlayer of bubbles on gold and YIG in PAO6 during the expansion phase of the bubble (a) as a function of time and (b) as a function of the ratio √𝜈 * 𝐿/𝑉 𝑒 with the

Table 1 :

 1 Value of the parameter n from the fit with X(𝑡) ∝ 𝑡 𝑛 , where X(t) is the time-dependent of the footprint radius of the bubble during its early expansion, i.e. the first quarter of the bubble lifetime (see supplementary Fig.S3 and S4). The values in brackets correspond to the 80% confidence interval according to the Student's t-distribution.

	n	Gold	YIG
	Water	0.389 (±0.051) 0.393 (±0.028)
	PAO6	0.332 (±0.073)	0.334 (±0.05)
	PAO40	0.194 (±0.035) 0.200 (±0.053)

Acknowledgement

This work has been financially supported by the European Cooperation in Science and Technology through COST Action MP1302 Nanospectroscopy. Tim Hupfeld thanks Evonik Industries for financial support. We would like to thank Carsten Dubs from INNOVENT e.V. Technologieentwicklung for providing the YIG wafer. BG and SB thank the Deutsche Forschungsgemeinschaft DFG for funding.

Supplementary Material

See supplementary material for more information on bubble height, bubble shrinking and the automated data processing for shadowgraph imaging. The supplementary material also provides information on the time dependence of the footprint radius and on the Rayleigh collapse time.


Yang Wang

Gaëtan Maurel

Marc Couty

François Detcheverry

Samy Merabia

Implicit medium model for fractal aggregate polymer nanocomposites: Linear viscoelastic properties

Dispersing solid fillers into a polymer matrix is a common strategy to enhance and tailor its properties. The polymer nanocomposites so obtained with fractal-like aggregates have exceptional rheological behavior that have long been exploited in the tire industry. However, due to the disparity of time and length scales, our understanding of the relation between nanocomposites structure and rheology remains far from complete. We propose in this work a mesoscopic model to address the dynamics and the rheology of aggregate nanocomposites. While aggregates are represented explicitly as groups of interacting particles, we use for the polymeric matrix an implicit description based on generalized Langevin and Stokes equations, that captures the average effect of a viscoelastic medium. These two-level modelling allows us to simulate large systems containing dozens of aggregates. Focusing here on the linear viscoelastic properties of PNCs, we characterize the influence of aggregate size, volume fraction, rigidity and polydispersity. We demonstrate that rigid aggregates systems display salient features of the phenomenology of PNCs, namely slow relaxation in the stress relaxation response, non-Maxwell elastic response at low frequency and supralinear dependence of their storage modulus with volume fraction. We also point out the critical role played by aggregate rigidity, and show that for either flexible aggregates or well dispersed nanofillers, the effects are far less spectacular. Our findings should help in designing nanocomposites with enhanced mechanical properties.

I. INTRODUCTION

Dispersing a tiny amount of solid nanofiller within a polymer matrix can yield a new material with strikingly distinct properties compared to the pure system [START_REF] Vilgis | Reinforcement of Polymer Nano-Composites: Theory, Experiments and Applications[END_REF][START_REF] Nielsen | Mechanical properties of polymers and composites[END_REF][START_REF] Mark | The Science and Technology of Rubber[END_REF]. The relentless interest in designing such polymer nanocomposites (PNCs) rests on the ever growing variety of fillers available [START_REF] Winey | Polymer nanocomposites[END_REF][START_REF] Kumar | 50th Anniversary Perspective: Are Polymer Nanocomposites Practical for Applications?[END_REF]. They come in many shapes (spherical nanoparticle, plate-like clays, graphene sheets, carbon nanotubes, or metallic nanorods) and they now extend over three decades in size, from the subnanometric OAPS additive [START_REF] Cheng | Big Effect of Small Nanoparticles: A Shift in Paradigm for Polymer Nanocomposites[END_REF], to the silica or carbon particles typically a few dozen of nanometers in diameter, which are the building blocks for aggregates and agglomerates reaching microns [START_REF] Vilgis | Reinforcement of Polymer Nano-Composites: Theory, Experiments and Applications[END_REF]. By tuning the functionalizations and adding grafted chains with well-chosen length and density, one may also control the state of dispersion, and generate large-scale, self-assembled structures such as sheets or strings [START_REF] Akcora | Anisotropic self-assembly of spherical polymer-grafted nanoparticles[END_REF], thus opening the way to a hierarchical arrangement of nanofillers [START_REF] Vaia | Polymer nanocomposites with prescribed morphology: going beyond nanoparticle-filled polymers[END_REF]. Given those countless possibilities, PNCs are being sought for a variety of applications, ranging from gas separations to electronic devices [START_REF] Winey | Polymer nanocomposites[END_REF][START_REF] Kumar | 50th Anniversary Perspective: Are Polymer Nanocomposites Practical for Applications?[END_REF][START_REF] Godovsky | Device Applications of Polymer-Nanocomposites[END_REF]. To fully come to fruition, the design of PNCs should ideally rest on a fundamental understanding of their properties. Yet, this endeavor has proved a recurrent challenge, and many puzzles remain [START_REF] Kumar | Perspective: Outstanding theoretical questions in polymernanoparticle hybrids[END_REF].

Though one of the oldest application, PNCs based on solid, fractal-like aggregates are no exception. Carbon black has been studied as an additive as early as the 1940s, later followed in the 1990s by silica aggregates. Both aggregates are built of up to hundreds of primary particles organized in ramified structure with size in the range 100 -500 nm. When introduced in a polymer or rubber matrix, even in moderate amount, they are known to yield mechanical reinforcement, enhance dissipative properties and augment shock and wear resistance [START_REF] Vilgis | Reinforcement of Polymer Nano-Composites: Theory, Experiments and Applications[END_REF][START_REF] Nielsen | Mechanical properties of polymers and composites[END_REF]. Such outstanding viscoelastic properties [START_REF] Khan | Dynamic rheological behavior of flocculated fumed silica suspensions[END_REF][START_REF] Mélé | Reinforcement effects in fractal-structure-filled rubber[END_REF][START_REF] Jouault | Well dispersed fractal aggregates as filler in polymer-silica nanocomposites : long range effects in rheology[END_REF][START_REF] Baeza | Multiscale filler structure in simplified industrial nanocomposite silica/SBR systems studied by SAXS and TEM[END_REF][START_REF] Cassagnau | Linear viscoelasticity and dynamics of suspensions and molten polymers filled with nanoparticles of different aspect ratios[END_REF] have long been exploited in tire applications [START_REF] Nielsen | Mechanical properties of polymers and composites[END_REF][START_REF] Mark | The Science and Technology of Rubber[END_REF]. Despite heavy commercial use and decades of research [START_REF] Heinrich | Recent Advances in the Theory of Filler Networking in Elastomers[END_REF][START_REF] Jancar | Current issues in research on structure-property relationships in polymer nanocomposites[END_REF][START_REF] Hall | Molecular theories of polymer nanocomposites[END_REF][START_REF] Song | Concepts and conflicts in nanoparticles reinforcement to polymers beyond hydrodynamics[END_REF][START_REF] Vogiatzis | Multiscale Molecular Simulations of Polymer-Matrix Nanocomposites[END_REF], our account of the relation between the PNC structure and its rheology is still incomplete. This lack of understanding is to be related to the competing mechanisms whereby nanofillers may affect the material behavior and contribute to reinforcement. Schematically, four classes of mechanisms may be distinguished. First, even at low volume fraction, the presence of solid additives can change the polymer dynamics, by shifting the spectrum of relaxation times [START_REF] Pryamitsyn | Origins of linear viscoelastic behavior of polymer -nanoparticle composites[END_REF], reducing the segment mobility, or increasing the entanglement density near the particle surface [START_REF] Jouault | Multiscale characterization of filler dispersion and origins of mechanical reinforcement in model nanocomposites[END_REF]. Second, chains that are surface-immobilized by adsorption, bonding, or entanglement may connect the particles into a long-lived, percolating filler-polymer network [START_REF] Sternstein | Reinforcement mechanism of nanofilled polymer melts as elucidated by nonlinear viscoelastic behavior[END_REF][START_REF] Moll | Mechanical reinforcement in polymer melts filled with polymer grafted nanoparticles[END_REF][START_REF] Thomin | Network Effects on the Nonlinear Rheology of Polymer Nanocomposites[END_REF][START_REF] Sen | Molecular underpinnings of the mechanical reinforcement in polymer nanocomposites[END_REF][START_REF] Akcora | Gel-like" Mechanical Reinforcement in Polymer Nanocomposite Melts[END_REF][START_REF] Masnada | Entanglement-induced reinforcement in polymer nanocomposites[END_REF]. Third, particles at high volume fractions might percolate or jam, thus ensuring efficient stress transmission [START_REF] Heinrich | Recent Advances in the Theory of Filler Networking in Elastomers[END_REF]. Such filler-filler network might involve the formation of a gel [START_REF] Chen | Mechanical reinforcement of polymer nanocomposites from percolation of a nanoparticle network[END_REF]. Finally, an alternative scenario focuses on filler-induced shift in the glass transition temperature [START_REF] Berriot | Gradient of glass transition temperature in filled elastomers[END_REF][START_REF] Montes | Influence of the Glass Transition Temperature Gradient on the Nonlinear Viscoelastic Behavior in Reinforced Elastomers[END_REF][START_REF] Merabia | A Microscopic Model for the Reinforcement and the Nonlinear Behavior of Filled Elastomers and Thermoplastic Elastomers (Payne and Mullins Effects)[END_REF]. Immobilized layers may form a network of glassy bridges that provides preferential stress-sustaining paths, and enhanced dissipation under strain-induced yielding. Which mechanism is relevant will heavily depend on the experimental conditions, physical chemistry details and composition of of PNCs. For aggregates, the parameter space is partic-ularly rich. The volume fraction, the morphology, size and polydispersity are all expected to play an important role, not to mention the polymer matrix itself.

Given the multiplicity of phenomena and the large parameter space, disentangling the contribution of each effect from experiments alone is fraught with difficulties. Accordingly, a large body of simulation work [START_REF] Jancar | Current issues in research on structure-property relationships in polymer nanocomposites[END_REF][START_REF] Vogiatzis | Multiscale Molecular Simulations of Polymer-Matrix Nanocomposites[END_REF][START_REF] Wang | Influence of Morphology on the Mechanical Properties of Polymer Nanocomposites Filled with Uniform or Patchy Nanoparticles[END_REF], has sought to obtain a microscopic view of the processes at play in PNCs. Molecular dynamics can offer a detailed description of the chain structure and dynamics in the vicinity of the nanofiller [START_REF] Pandey | Polymer/Nanoparticle Interactions: Bridging the Gap[END_REF]. By resorting to more coarse-grained approach, such as DPD simulations [START_REF] Lahmar | Onset of entanglements revisited. Dynamical analysis[END_REF] or slip-links models [START_REF] Masnada | Entanglement-induced reinforcement in polymer nanocomposites[END_REF], one can address PNCs systems comprising a small number of nanoparticles while keeping an explicit description of the polymeric chains.

Not surprisingly, the case of PNCs comprising large aggregates has not been tackled so far, as the length scale required remains out of reach by microscopic methods. Simulations at the scale of a single aggregate 100 nm in radius can already be quite demanding in computation resources. Moreover, a faithful representation of the PNCs requires that many aggregates be included, so that structures appearing at a higher level, such as a network, can also be captured. Only in large systems containing dozens of aggregates can the role of morphology, size and polydispersity be investigated. Given the wide disparity of time and length scale present in the system, an explicit representation of the chains is currently impossible.

We propose in this work a mesoscopic model capable of reaching the length scale relevant to fractal aggregate PNCs. The idea is to climb one more step in the level of coarse-graining, by resorting to an implicit description of the polymer medium, for which no degrees of freedom are allocated explicitly. Rather, its presence is accounted in an effective manner through a non-Markovian friction force. This force depends on the history of the velocity of the filler and on a memory kernel that directly embodies the viscoelastic behavior of the polymeric matrix. In contrast with the surrounding medium, the aggregates are represented explicitly in our simulation. We are thus free to chose their size, interaction and morphology. Our approach is in line with recent efforts to model the behavior of jack-shaped and fractal aggregates [START_REF] Petersen | Shear rheology of extended nanoparticles[END_REF][START_REF] Heine | Effect of particle shape and charge on bulk rheology of nanoparticle suspensions[END_REF][START_REF] Seto | Hydrodynamic stress on small colloidal aggregates in shear flow using Stokesian dynamics[END_REF] suspended in a purely viscous fluid, but addressing the case where the matrix is viscoelastic.

With this model in hand, we investigate the reinforcement and viscoelastic properties of aggregate PNCs. We find that even moderate volume fractions of aggregates strongly enhance the stress response of the model nanocomposite. We examine the influence of aggregates volume fraction, size, rigidity and polydispersity. We also describe the slow relaxation in stress response and provide an interpretation based on the aggregate rotation. We focus exclusively on linear rheology. A subsequent study will be devoted to non-linear effects.

The remainder of this article is structured as follows. We introduce in Sec. II the generalized Langevin equation and generalized Stokes relation that underlie our approach, and explain the numerical method employed. Section III describes our choice of aggregates, their morphology and interactions, together with two other fillersflexible and individual nanoparticles-as reference points. Results and discussion are presented in Sec. IV.

II. MODEL AND METHOD

The central idea of our model is an implicit description of the polymer matrix. Given the length scale targeted -not only a single aggregate but a whole distribution of them -, the explicit representation of polymeric chains is challenging, to say the least. This is obviously the case with molecular dynamics, which would require immense computational resources, but also with more coarse-grained approaches such as slip-links [START_REF] Biondo | Numerical study of a slip-link model for polymer melts and nanocomposites[END_REF] or DPD [START_REF] Lahmar | Onset of entanglements revisited. Dynamical analysis[END_REF]. Therefore, there is no choice but to go one step up in the ladder of coarse-graining levels. One possibility to do so is to forget all degrees of freedom of the medium, by resorting to an implicit description. The effect of the medium on an embedded particle is then reduced to two distinct forces: a friction force that represents the average drag, and a random force that originates in thermal fluctuation. From Brownian dynamics to DPD, this method has been widely used when the medium is a viscous liquid. As we detail in this section, the same approach can be used when the medium is a viscoelastic polymer. For simplicity, we first focus on the motion of a single particle. Extension to our model aggregates immediately follows since they are represented as a group of particles.

A. Implicit description of medium by GLE

The generic tool to describe the motion of a particle in an implicit medium is the generalized Langevin equation (GLE). For a particle of mass m, position r and velocity v, it reads as

m dv(t) dt = F (r(t)) - t -∞ Γ(t -t ′ )v(t ′ ) dt ′ + F r (t). ( 1
)
The first term F in the right-hand side is the conservative force accounting for the interactions with other particles. The second term in the RHS of Eq. ( 1) is the drag force F d exerted by the medium. Importantly, it does not only depend on the instantaneous velocity but involves all values taken in the past, as weighted by the memory kernel Γ(t). Such non-locality in time makes the evolution equation non-Markovian. Finally, the last term F r is the random force, whose time correlation function is nontrivial, since it is directly related to the memory kernel by the fluctuation-dissipation theorem

F r (t) F r (t ′ ) = k B T Γ(t -t ′ )δ, (2)
with k B the Boltzmann's constant, T the temperature and δ the identity matrix. The simple Langevin equation is recovered as a particular case for a kernel instantaneous in time Γ(t) ∼ δ(t).

The influence of the medium on particle dynamics is entirely encoded in the memory kernel, which should therefore depends on the material properties. For a viscoelastic polymer matrix, the connection is provided by a generalized Stokes law [40]. The kernel Γ(t) is thus directly proportional to the stress relaxation modulus G p (t) of the polymer matrix. For a spherical particle of radius R,

Γ(t) = 6πR G p (t). ( 3
)
The case of a viscous fluid with viscosity η leads back to the usual Stokes law Γ(t) = 6πRηδ(t). Equation (3) derives from a correspondence principle [START_REF] Indei | Treating inertia in passive microbead rheology[END_REF]. To find the flow of a viscoelastic medium, it suffices to consider the equivalent problem with a purely viscous fluid, and working in the frequency domain (ω), to substitute the viscosity η with G * p (ω)/iω, where G * p (ω) is the complex modulus of the material. Because G * p (ω)/iω is the Fourier transform of G p (t), one immediately gets Eq. (3).

Though its derivation is exact, the generalized Stokes equation is based on a number of approximations [START_REF] Indei | Treating inertia in passive microbead rheology[END_REF]: (i) The medium is described as a continuum. (ii) The medium is incompressible. (iii) Only viscoelastic properties in the linear regime are relevant, implying that the medium remains in the vicinity of equilibrium. (iv) The effect of bead and medium inertia are both neglected. In particular, this includes the Basset force that would be relevant at very high frequencies. For the nanocomposite systems considered in this work, assumptions (ii) to (iv) are well justified. As regards requirement (i), it certainly applies for large aggregates whose size is above entanglement length, but it would break down for a nanoparticle [START_REF] Yamamoto | Theory of nanoparticle diffusion in unentangled and entangled polymer melts[END_REF][START_REF] Yamamoto | Spatially dependent relative diffusion of nanoparticles in polymer melts[END_REF][START_REF] Yamamoto | Microscopic Theory of the Long-Time Diffusivity and Intermediate-Time Anomalous Transport of a Nanoparticle in Polymer Melts[END_REF]. For small aggregates, Eq. (3) should be seen as a zero-order approximation. As a minor point, we also note that a no-slip boundary condition is assumed at the particle surface, which does not necessarily hold for polymers (see discussion in [START_REF] Pryamitsyn | Mechanisms of steady-shear rheology in polymer-nanoparticle composites[END_REF]), but this assumption is easily relaxed [55].

B. Dynamics and numerical method

Having defined our model for the medium, we now explain how it can be put to work in practice. To describe the particle motion in the viscoelastic matrix, it is necessary to solve the GLE numerically. With respect to the simple Langevin equation used in Brownian dynamics, Eq. (1) differs in two ways. First, instead of taking independent values at each time step, the noise is correlated, thus requiring a sequence of correlated random numbers. Second, instead of involving only the current velocity, the drag force is obtained by convoluting the memory kernel with all past values of the velocity. For a large number of particles, the storage of velocity histories and the convolution will become demanding in memory and computation time. Fortunately, Baczewski and Bond have recently proposed a method to circumvent this issue [START_REF] Baczewski | Numerical integration of the extended variable generalized Langevin equation with a positive Prony representable memory kernel[END_REF]. The principle is to introduce new variables to rewrite the non-Markovian GLE in an expanded state space where the evolution becomes Markovian.

The method is applicable on the condition that the memory kernel can be developed as a Prony series, i.e. a sum of exponentially decaying terms [56]. The stress relaxation modulus is thus written as

G p (t) = M m=1 G p,m exp - t τ m , (4)
where G p,m and τ m are the amplitude and decay time of mode m, and M the total number of modes. The coefficients Γ m = 6πRG p,m of the memory kernel are denoted as c m /τ m for convenience. For simplicity, let us consider in the GLE equation only one vectorial component of the velocity. For each mode, two new variables are introduced. The first is

Z m (t) = - t 0 c m τ m exp - t -t ′ τ m v(t ′ ) dt ′ , (5)
which satisfies

τ m dZ m (t) = -Z m (t) dt -c m v(t) dt. ( 6
)
The second intermediate variable F r,m is the contribution of mode m to the random force, defined from

τ m dF r,m (t) = -F r,m (t) dt + 2k B T c m dW m (t), (7)
where W m (t) is a standard Wiener process. Importantly, since F r,m obeys an Ornstein-Uhlenbeck process, its time correlation function is an exponential, with characteristic time τ m :

F r,m (t) F r,m (t ′ ) = k B T c m τ m exp - t -t ′ τ m . (8)
Finally, using the total random force

F r (t) = M m=1 F r,m (t), one finds that the variable S m (t) = Z m (t) + F r,m (t) satisfies m dv(t) = F c (t) dt + M m=1 S m (t) dt, (9)
τ m dS m (t) = -S m (t) dt -c m v(t) dt + 2k B T c m dW m (t). ( 10
)
In contrast with the original GLE, those equations no longer involve time convolution: they are local in time, since they only require the current value of each variable.

In practice, the Baczewksi and Bond method has been implemented in the molecular dynamics open code LAMMPS [START_REF] Plimpton | Fast parallel algorithms for short-range molecular dynamics[END_REF], which can be used to solve the GLE dynamics.

To summarize the approach so far, the implicit description of the medium eliminates all associated solvent degrees of freedom but to solve the resulting non-Markovian GLE, it is convenient to make the evolution equation Markovian by introducing additional variables. If there are N p particles in a three-dimensional system and M modes in the memory kernel, their total number is 3M N p , which in general is much smaller than the number of variables required by an explicit description of the chains, hence the interest of the method.

It remains to discuss the choice of the memory kernel. For a viscoelastic matrix, there are actually several routes possible. The first is to rely on an analytical model of the polymeric melt. In this regards, we note that both the Rouse model and reptation model lead to relaxation modulus that are Prony series [START_REF] Rubinstein | Polymer physics[END_REF]. The amplitude and characteristic time of each mode are then known exactly [57]. The second possibility to specify the modes in the Prony series is simply to fit a given relaxation modulus, with a number of modes determined by the accuracy required. The target G p (t) may come directly from experimental data, thus allowing to describe a specific polymer matrix. Alternatively, G p (t) may also be obtained from simulations where the polymer chains are represented explicitly (such as molecular dynamics, dissipative particle dynamics or slip-links models). By measuring the correlation of the random force acting on an embedded particle, one immediately gets the memory kernel of the medium, thanks to the fluctuation-dissipation theorem Eq. ( 2). The ability to integrate information obtained at a lower level of coarse-graining makes the implicit medium approach an attractive scheme for a multi-scale strategy.

III. AGGREGATES

The second ingredient of our model is an explicit description for the aggregates. They are represented as groups of interacting particles, whose size, morphologies and interaction can be chosen at will. There is thus considerable freedom in this modelling step. Here we describe our particular choice for the ramified aggregates used in this work.

A. Aggregate model

We model an aggregate as a group of repulsive particles connected through a spring network that maintains its shape and integrity. Specifically, all particles interact through a repulsive potential U R defined as

U R (r) = U LJ (r) -U LJ (r m ) if r < r m , (11)
= 0 otherwise,
where

U LJ (r) = 4ǫ (σ/r) 12 -(σ/r) 6 (12)
is the Lennard-Jones (LJ) potential, r is the distance between the particles, ǫ is the depth of the potential, and σ (resp. r m ) is the distance at which the potential is zero (resp. minimum). Particles are seen as spheres of diameter σ that are thus in contact when r = σ. As shown in Fig. 1, an aggregate is represented as a collection of particles at contact, endowed with two sets of springs. On the one hand, the connective springs join neighboring particles, with a quadratic potential

U S (r) = k 2 (r -l c ) 2 , ( 13
)
where k is the spring constant and the rest length l c is chosen so as to keep spheres at contact [58], thus preserving the aggregate connectivity. On the other hand, a set of virtual springs is added to preserve the aggregate geometrical structure. Whereas we are primarily interested in the rigid aggregates defined above, it is useful for comparison purpose to introduce two other kinds of systems. First, the individual particles, devoided of any springs, will provide a reference point. Second, we also define the flexible aggregates, which include only connective springs. Because their shape can fluctuate, they are reminiscent of star polymers. They will help to assess the role of aggregate rigidity in the rheology of the nanocomposite. For all type of filler, the volume fraction φ of filler is defined as

φ = πN p σ 3 6V , (14)
where N p is the total number of particles and V the volume of the system.

B. DLA morphology

We have chosen to consider aggregates generated by a diffusion-limited aggregation (DLA) algorithm [START_REF] Sander | Diffusion-limited aggregation: a kinetic critical phenomenon?[END_REF]. Briefly, the process starts with an aggregate of one particle at the system center. A particle released far away diffuses until it touches the aggregate and irreversibly becomes part of it. The process is then repeated many times. In the limit of a large number of particles, the aggregates have a branching, disordered structure [START_REF] Sander | Diffusion-limited aggregation: a kinetic critical phenomenon?[END_REF][START_REF] Seager | Slippery diffusion-limited aggregation[END_REF], whose average coordination number is 2.0, and with fractal dimension d f = 2.5. Our aggregates are generated by stopping the DLA process when it reaches a prescribed number of particles N . Because only moderate sizes are considered (N ≤ 100), the aggregates are not genuine fractal objects, but nonetheless exhibit a loose, ramified structure. To better characterize their properties, it is useful to in-troduce the gyration tensor

S = 1 N N n=1 x n x T n , (15)
Here, x T n is the transpose of x n = r nr cm , with r n the position of the n th particle, and r cm the position of the center of mass. The spatial extent of an aggregate is given by the radius of gyration R g defined from

R 2 g = Tr S = 1 N N n=1 (r n -r cm ) 2 . ( 16
)
The eigenvectors and eigenvalues (λ

2 1 > λ 2 2 > λ 2
3 ) of the gyration tensor define the orientation and the dimension of an equivalent ellipsoid. Shown in Fig. 2 are the radius of gyration for various aggregates sizes. The distribution indicates that most of the aggregates share a size close to the mean value, which exhibits a modest increase, from R g = 2.2 to 3.3 for N = 20 to 50. Besides, the aggregates are anisotropic objects as revealed by the aspect ratio κ = λ 1 /λ 3 , whose distribution shows significant departure from unity.

C. Polydispersity

We will examine both monodisperse and polydisperse ensembles of aggregates. In the former, all aggregates have the same size N [60]. In the latter, which is the relevant situation in experiment, the size distribution is governed by a probability distribution P (N ), which can be quite broad [START_REF] Baeza | Multiscale filler structure in simplified industrial nanocomposite silica/SBR systems studied by SAXS and TEM[END_REF]. For the sake of simplicity, we use a gamma distribution

P (N ) = N α Γ(α + 1)β α+1 exp - N β , ( 17
)
where Γ is the gamma function and α and β are free parameters. The most frequent aggregate size is then N * = αβ and the mean is N = (α + 1)β. As an illustration, the distribution given by Eq. ( 17) is plotted in Fig. 3, with α = 2 and β = 10, leading to N * = 20 and N = 30 and a size polydispersity around 30%. This is the default case, unless otherwise mentioned. Note that in practice, the largest aggregates generated have size N = 100.

D. Initial configurations

Once the aggregates are generated and their size distribution chosen, it remains to obtain the initial configurations by placing them in the simulation box. We used two different methods to do so. In the insertion method, the aggregates are placed one after another. Random positions and orientations are attempted for the new aggregate until it has no overlap with any neighbors. For polydisperse aggregates, it is more efficient to place the aggregates in order of decreasing size. In the compression method, the aggregates are initially placed in a very large box, whose dimensions are then gradually reduced. At each compression step, the system can relax so as to minimize overlap between aggregates. Applying this process to rigid aggregates may change their shape. To make sure that they keep their original structure, aggregates are treated as undeformable objects during the compression process [START_REF]The springs forces thus do not play any role[END_REF]. Their relaxation thus involve only rigid motion, through translation and rotation induced by the total force and torque. In both initialization methods, the aggregates once placed are subject to an additional relaxation process governed by Brownian dynamics, which is stopped when the global energy and pressure of the system have converged to constant values.

For individual particles and flexible aggregates, which can deform to accommodate constraints, it is easy to generate initial configurations with high volume fractions. This is much more difficult with rigid aggregates. It seems that the compression method works best with monodisperse aggregates and the insertion method with polydisperse systems. Shown in Fig. 4 are the volume fractions accessible for monodisperse aggregates as a function of their size N and the number N a of aggregates in the box. For small aggregates (N = 20), the highest volume fraction is φ = 17%, whereas for larger aggregates (N = 50), it decreases to φ = 10%. Except at the highest volume fraction accessible, the total number of aggregate in the box N a can be varied from 20 to 50, suggesting that finite size effects are weak. For polydisperse systems, the densest available systems have φ in the range 18 -15%, as the aggregate number N a increases from 200 to 1000. Thus, our simulated systems do not reach the highest nominal volume fractions (up to 30%) that are reported in experiments. Note however that our volume fractions are not directly comparable to that used with experimental systems [62]. Accordingly one should be cautious in comparing absolute values and should focus primarily in the trends induced by changing the volume fraction. φ (%)

N a =20 N a =30 N a =40 N a =50
FIG. 4: Volume fraction of initial configuration accessible by the compression method, as a function of aggregate size N for different number Na of aggregates in the system.

IV. RESULTS AND DISCUSSION

In this section, we will examine how the linear rheology is influenced by the filler and matrix parameters. We first briefly describe our reference system and how the results were obtained.

A. Preliminaries

System studied and parameters

As explained in Sec. II A, the viscoelasticity of the matrix can be chosen so as to match experimental data for particular polymer, which allow us to examine a specific nanocomposite. The corresponding Prony series would involve multiple modes. Here, our goal is to understand, generically, the role of the polymer viscoelasticity and the interplay with filler characteristics. Besides, the parameter space is already quite large. For these two reasons, we have preferred to start with a simpler situation, where the medium is a Maxwell fluid. The polymer relaxation is thus assumed mono-mode

G p (t) = G 0 p exp - t τ 0 p , (18)
where G 0 p is the plateau stress modulus and τ 0 p the terminal relaxation time of the polymer matrix. For a typical system -a polybutadiene of molecular weight 40 K [START_REF] Zhu | Investigating Linear and Nonlinear Viscoelastic Behavior Using Model Silica-Particle-Filled Polybutadiene[END_REF] -, one gets G 0 p = 0.5 MPa and τ 0 p = 1 ms. It remains to specify our choice of units and numerical parameters. Unless mentioned otherwise, the unit length is the particle diameter σ, the unit time is τ 0 p and the unit energy is k B T . Accordingly, the pressure unit is P u = k B T /σ 3 . Throughout the study, we consider particles with diameter 9.4 nm. For the system considered above, and always working at room temperature, this implies G 0 p = 100 P u . We fix ǫ = 1 throughout the study, a choice expected to be inconsequential since only repulsive interaction are taken into account. Besides, we note that the particle mass enters the dynamics only through its combination with a time step (see Eq. ( 10)). For the sake of numerical efficiency, the mass is not set to a realistic value (we always set m = 1 in reduced units) but serves only to set the value of the time step. We use ∆t = 0.01 τ 0 p . The units and some typical parameters are summarized in Tab. I. Finally, our simulated systems contains typically 10 3 particles, thus including dozens of aggregates (N a in the range 10-50 depending on the aggregate size). We have checked that finite size effects, though sometimes detectable, would not change significantly the results presented below.

Stress relaxation modulus

The stress tensor is computed from the position and conservative forces between particles

σ αβ = - 1 V i<j F ij,α R ij,β . (19)
Here, F ij,α is the α component of the interaction force applied on particle j by particle i, R ij,β is the β coordinate of the vector joining particle i to particle j, and indices i and j run over all N p particles in the system. It is important to note here that the interaction force which appears in the expression of the stress is the conservative force, which derives from a pair potential

F ij,α = - ∂U ∂R ij,α (20)
where U (R) is the interparticle potential that includes the repulsion U R (r) and the spring potentials U S (r) relevant to the case of aggregates (Eqs. ( 11) and ( 13) respectively).

In the linear regime, the stress relaxation modulus G(t) can be obtained from the time correlation of the shear stress, via the linear response theory [START_REF] Chandler | Introduction to modern statistical mechanics[END_REF],

G(t) = V k B T 1 3 3 α=1 3 β>α σ αβ (t)σ αβ (0) , (21)
where . . . denotes an ensemble average. The complex modulus G * (ω) = G ′ (ω) + iG ′′ (ω), where G ′ and G ′′ are the storage and loss modulus, is obtained from the G(t) by the Fourier transform of G(t) as

G * (ω) = iω ∞ 0 G(t) exp (-iωt) dt. ( 22
)
Two others remarks apply to all the data shown below. i) The stress relaxation modulus G(t) and other computed quantities all result from an average over at least 10 independent simulations. This is necessary to improve the statistics, especially for the behavior of the observables at long times. ii) In all curves, including the stress relaxation modulus or the dynamic moduli, only the contribution of the filler is shown. The contribution of the polymer matrix is not taken into account.

B. Influence of aggregate properties

Effect of filler type

The stress relaxation modulus G(t) is shown in Fig. 5 (left) for the three types of solid fillers and a volume fraction φ = 10%. One common feature is the existence of an oscillatory response in the early regime. Even though it may be influenced by the springs properties, it is already present in the nanoparticles alone and seems to originate primarily from the interplay between inertia and the elastic behavior of the polymer matrix at short time. Whatever the filler type, oscillations occurs at a time t ≃ 0.1 significantly smaller than the polymer relaxation time [63] and therefore do not affect the intermediate and long-time behavior which are the focus of this work.

Nanoparticle and aggregate systems exhibit striking differences in their stress relaxation modulus, both in terms of relaxation time and level of reinforcement. Let us consider the former first. In the case of nanoparticle, the decay of G(t) is well described by a single exponential, whose characteristic time is of the order of τ 0 p . In contrast, both the flexible and the rigid aggregates display a mechanical response whose decay is not purely exponential at long times but better described by a stretched exponential G(t) ∼ exp(-(t/τ ) β ) with β ≃ 0.4, in line with the idea that the response involves a distribution of relaxation times, some of them much larger than τ 0 p . The origin of such large relaxation times will be discussed in Sec. IV C in relation with the microscopic motion of the particles. As regards reinforcement, the effect of filler type is also clearly seen. At short and intermediate times, the presence of flexible and rigid aggregates induces a mechanical response which exceeds by respectively one and two orders of magnitude that observed with nanoparticles. Perhaps surprisingly, the difference of G(t) between flexible and rigid aggregates subsides at long times, suggesting that in this regime the effect of connectivity tends to dominate over rigidity.

The effect of the filler type may be equally evidenced by looking at the storage and loss moduli G ′ and G ′′ , as displayed in Fig. 5 (right). Because those quantities are often measured in experiments, and carry the same physical content as the stress relaxation modulus, they are used from now on to present our results. Consistent with the Maxwell model, the simple exponential decay in the stress relaxation of the nanoparticle system now translates into a power-law dependence ω 2 and ω 1 at low frequency for G ′ and G ′′ respectively. In the aggregate systems, power-law dependence are also seen at low frequencies, but with different exponents. Though numerical noise and a limited range of data prevent an accurate estimate, the exponents are unambiguously below their Maxwell values, and seem to be rather similar both for G ′ and G ′′ and for the two types of aggregates. As above, it is clear that at all frequencies the aggregates systems display a large reinforcement with respect to individual nanoparticles. Below, we focus primarily on the behavior of rigid aggregates, since they are more relevant to experimental systems.

Two remarks are in order. i) It is worth pointing our that in our model, the only way particles pertaining to different aggregate may interact is by pure repulsion. The absence of any polymer-mediated interactions that could induce polymer-bridging between particles does not prevent the existence of relaxation times much longer than that of the polymer matrix. ii) Long relaxation times have been observed experimentally in the stress relaxation modulus of polybutadienes filled with individual silica nanoparticles (see Fig. 2 of Ref. [START_REF] Zhu | Investigating Linear and Nonlinear Viscoelastic Behavior Using Model Silica-Particle-Filled Polybutadiene[END_REF]). The existence of long relaxation times is often attributed to polymer bridging between neighbouring particles. Alternatively, our results suggest that aggregated particles may be the cause of slow stress relaxation, even in the absence of polymer bridging.

Effect of volume fraction

The filler volume fraction is known to have a major influence on the mechanical and rheological properties of PNCs [START_REF] Baeza | Multiscale filler structure in simplified industrial nanocomposite silica/SBR systems studied by SAXS and TEM[END_REF][START_REF] Jancar | Current issues in research on structure-property relationships in polymer nanocomposites[END_REF][START_REF] Akcora | Gel-like" Mechanical Reinforcement in Polymer Nanocomposite Melts[END_REF]. Figure 6 (top) shows the storage modulus for various volume fractions up to 17% in nanoparticles and rigid aggregates systems. The stronger reinforcement induced by the latter is again visible: the response of nanoparticle PNC remains below that of aggregates, even when the loading is threefold higher (17% vs 5% respectively). Besides, in the case of nanoparticles, the change in volume fraction induces only a shift of G ′ (ω) over the whole range of frequencies. With rigid aggregates, the behavior is modified at low frequencies. While it remains power law, G ′ ∼ ω ν , the exponent decreases from approximately 2 to 1/2, when the volume fraction increases from 5 to 17%.

The relative enhancement of the storage modulus with the volume fraction is further illustrated in Fig. 6 (bottom). Two distinct behaviors may be distinguished depending on the frequency. At relatively high values ω > 0.1, the increases in modulus is almost linear with the volume fraction: G ′ ∝ φ. In contrast, for small ω, because the exponent changes with φ, the increase of G ′ is clearly superlinear. For instance, when ω = 0.01, one finds G ′ ∝ φ m , with m ≃ 2.4. One can expect an even sharper variation, and a higher exponent at lower frequencies. This phenomenology is typical of PNCs [START_REF] Song | Concepts and conflicts in nanoparticles reinforcement to polymers beyond hydrodynamics[END_REF].

The experimental observation of superlinear behavior is commonly interpreted as the signature of percolation [START_REF] Song | Concepts and conflicts in nanoparticles reinforcement to polymers beyond hydrodynamics[END_REF], and the existence of mechanical paths pertaining to transmit the stress across the system. This interpretation is in line with the possible formation of a gel mediated by physical interactions between fillers. Our simulations show that this phenomenon occurs for fractal repulsive aggregates. Conversely, for individual nanoparticles no superlinear behavior is seen.

Effect of aggregate size

The number N of particles in an aggregate controls its spatial extent, as illustrated in Fig. 2 with the radius of gyration and aspect ratio. How the aggregate size af- fects the storage modulus is shown in Fig. 7. Clearly, the larger the aggregates, the stronger the mechanical response. However, we note that the low frequency exponent G ′ ∼ ω ν seems to remains relatively independent of size, at least when N > 5. In this regime, increasing the aggregate size or volume fraction has thus qualitatively different consequences. For high and intermediate frequencies, the situation is different. The level of reinforcement is best visualized by representing G ′ as a function of the volume fraction, as done in Fig. 8 for various aggregate size. Two behaviors are observed, depending on the frequency ω analyzed. At high frequency ω = 1 (top), the storage modulus of all the systems considered increases almost linearly with volume fraction :

G ′ (φ) ∝ φ.

Effect of polydispersity

So far, we have considered only monodisperse aggregates characterized by a fixed number of particles N . However, aggregates in real polymer nanocomposites may be polydisperse, which may affect the reinforcement. To address this point, we show in Fig. 9 the storage modulus of a polydisperse system with polydispersity around 30%. The corresponding size distribution is represented in Fig. 3, and has mean N = 20 and most frequent size N * = 20. The storage moduli of monodisperse system with size N and N * are also reported for comparison. Whatever the frequency, the mechanical response of the polydisperse system turns out to be slightly stronger. This effect may be ascribed to the existence in the polydisperse system of a few large aggregates (N > 30), which may have a significant effect on the rheology of the model PNC.

C. Slow relaxation times

Structure factors

We gave evidence in Sec. IV B of the long relaxation times which emerge in the stress relaxation function in the presence of aggregates (Fig. 5). Here, we aim at providing a microscopic analysis of the slow dynamics of fillers. To this end, we probe the particle dynamics through several observables. Specifically, we will consider the standard, self and intra-aggregate dynamic structure factors, that are based on the time correlation of different microscopic density fields.

The standard dynamic structure factor characterizes the collective dynamics of the particles S(q, t) = 1 N p ρ(q, t)ρ * (q, 0) ,

since it involves the density field of the whole system ρ(q, t) = Np i=1 exp(-iq•r i (t)), with N p the total number of particles, q the wave vector and r i the position vector of particle i. The self dynamic structure factor probes the motion of individual particles

F (q, t) = 1 N p Np i=1 ρ i (q, t)ρ * i (q, 0) , (24)
where ρ i (q, t) = exp (-iq • r i ) is the local density associated to particle i. Last, the intra-aggregates structure factor, is useful to probe the relationship between the local density field and the motion of the aggregates

S intra (q, t) = 1 N a Na k=1 ρ k (q, t)ρ * k (q, 0) . (25)
Here, N a is the total number of aggregates, and

ρ k (q, t) = N k n=1 exp (-iq • r n (t)
) is the local density field of aggregate k, which contains a number N k of particles. The observable S intra is useful to probe the change in the dynamics due to the interaction between the aggregates.

The three dynamic structure factors are shown in Fig. 10 for the different types of fillers. Two values of the wave vectors have been considered: q = 2π corresponds to the primary particle scale, and q = 2π/3 provide hints on the dynamics at the scale of a single aggregate (whose radius of gyration is R g ≃ 3, see Fig. 2). Focusing first on S(q, t), we see that at the particle scale, individual nanoparticles and flexible aggregate systems relax their density with comparable decay times on the order of τ 0 p . The local dynamics of rigid aggregates is slower, and displays long time tails that are absent from other systems. Thus, the existence of these long relaxation times does not derive from the connectivity but primarily from the rigidity of the aggregates. When considering the scale of the aggregates q ≃ 2π/R g , the difference between flexible and rigid aggregates is even more pronounced: rigid aggregates relax their density field after a time approximately one order of magnitude larger than the relaxation time characterizing the flexible aggregates. This observation is again in line with the conclusion that the rigidity of the aggregates is key to slow dynamics. Finally, it is clear that, whatever the scale analyzed, the relative difference between the density relaxation of flexible and rigid aggregates appears to be independent of 23), [START_REF] Moll | Mechanical reinforcement in polymer melts filled with polymer grafted nanoparticles[END_REF], and (25), and computed for wave vector q = 2π and q = 2π/3 ≃ 2π/Rg, where Rg is the aggregate radius of gyration.

the dynamic structure factor considered. In particular, for rigid aggregates, the relaxation probed by the standard, self and intra-aggregate structure factor is almost the same. As the intra-aggregate structure factor is a measure of the relaxation of the aggregate density field, we shall conclude that interactions between aggregates are the primary cause of the slow dynamics highlighted in the dynamic structure factor, and intra-aggregate interactions play a subdominant role.

Orientational dynamics

The structure factor considered so far to probe the system dynamics all rely on densities. To investigate further the microscopic internal dynamics of the aggregates, we now examine their individual rotational motion. Whereas flexible aggregates would also feature deformation, the motion of rigid aggregates (that are nearly undeformable) may be decomposed in two contributions: diffusion of the center of mass and rotation with respect to it. To measure the latter, we have used the gyration tensor given in Eq. ( 15) to define for each aggregate an instantaneous orientation vector u as the eigenvector with the largest eigenvalue and computed its correlation function C(t) = u(t) • u(0) [64]. The result averaged over all aggregates is shown in Fig. 11. In all cases considered (aggregate size of N = 20 and 50), rigid aggregates have a slower orientational dynamics than flexible aggregates having the same size N . Note that the typical relaxation times associated to orientation are more than one order of magnitude larger than those associated to stress relaxation modulus. The orientation decorrelation times are also larger than the decay times of the dynamic structure factors probed at the aggregate length scale. Thus the densities field at the particle, aggregate or global scales that all relax much faster than the orientation. In conclusion, we have highlighted the very slow internal rotation of the aggregates. The corresponding kinetics is slower than the relaxation of G(t), implying than the slow internal motion does not participate to the stress relaxation of the system.

V. CONCLUSION

We introduced a mesoscopic model aimed at describing the dynamics and the rheology of polymer composites containing fractal-like aggregates. Resorting to an implicit description of the polymer matrix through the generalized Langevin and Stokes equations allows to reach the relevant time and length scales, while at the same time keeping an explicit representation of fillers. We have shown how the filler type, loading, size, rigidity and polydispersity can all affect the mechanical response. Our simulations reproduce salient features of the phenomenology of PNCs: existence of long relaxation times in G(t), power law behavior G ′ ∼ ω a with an exponent a < 1 and superlinear increase of reinforcement with loading: G ′ ∝ φ m , with m in the range 1.7 -2.5. These features are all observed experimentally [START_REF] Song | Concepts and conflicts in nanoparticles reinforcement to polymers beyond hydrodynamics[END_REF][START_REF] Zhu | Investigating Linear and Nonlinear Viscoelastic Behavior Using Model Silica-Particle-Filled Polybutadiene[END_REF]. Importantly, our study shows the importance of the aggregation and rigidity of the aggregates: for individual well dispersed nanoparticles, the suspension linear rheology is more akin to a Maxwell like model described by the polymer matrix terminal time. Superlinear behavior is primarily dependent on the aggregate rigidity: the storage modulus of flexible aggregates increases linearly with loading. We have finally interpreted the emergence of slow dynamics in rigid aggregates systems as resulting from interaggregate interactions.

It is worth pointing out that our model is minimal in the sense that it includes only two basic ingredients: i) the effect of viscoelastic polymer matrix on the particle motion and ii) the repulsive interactions between aggregates. Accordingly, among the mechanisms put forward to rationalize the PNCs behavior, several are absent: no polymer-mediated network between fillers or glassy bridges can occur. Yet, for the fractal-like rigid aggregates considered here, the filler-filler repulsive interactions and the presence of a viscoelastic medium are sufficient to generate a rich phenomenology typical of PNCs.

There are several possible extensions and directions that remain to be explored within the present framework. First, we considered a single exponential memory kernel as a simple proxy to the rheology of a polymer melt. It is clearly possible, through the memory kernel, to obtain the kernel from more microscopic calculations, where the polymer degrees of freedom are explicitly taken into account. This would open the way to build a multiscale description of the rheology of PNCs with aggregated fillers. Second, we focused on DLA-generated aggregates but it is straightforward to examine any type of morphology [START_REF] Heine | Effect of particle shape and charge on bulk rheology of nanoparticle suspensions[END_REF][START_REF] Knauert | The effect of nanoparticle shape on polymer-nanocomposite rheology and tensile strength[END_REF]. Likewise, filler interactions were repulsive only but the effect of attraction between particles, originating in physical forces or short-range interactions, deserves to be investigated. This is work in progress. Finally, as in Brownian dynamics, hydrodynamic interactions are neglected. Accounting for such generalized (viscoelastic) interactions is an interesting perspective that is left for future work. We anticipate though that aggregates would interact with each other before contact and the presence of aggregates would modify the mechanical response at even lower volume fraction.

the volume fraction. To make contact with experimental value, one should probably choose a value in between.

[63] Remember that τ 0 p is taken as the unit time.

[64] Remember that we consider only the linear regime and that our simulations are at equilibrium.

FIG. 2 :

 2 FIG. 2: Characteristics of our DLA aggregates. (Top) Aggregate of size N = 20 and 50. (Middle) Distribution of radius of gyration Rg for several aggregate sizes. The average Rg is shown in the inset. (Bottom) Distribution of aspect ratio κ, as defined in the text.

FIG. 3 :

 3 FIG.3: Aggregate size distribution P (N ) given by Eq. (17) with α = 2 and β = 10. The most frequent aggregate size is N * = 20 and the mean size is N = 30.

2 FIG. 5 :

 25 FIG.5: Effect of the filler type on the stress relaxation modulus (left) and storage and loss moduli (right, filled and empty symbols respectively). The aggregates are monodisperse with size N = 50 and the volume fraction is φ = 10%.

FIG. 6 :

 6 FIG. 6: Effect of volume fraction on the storage modulus. (Top) G ′ (ω) for rigid aggregates of size N = 20 and for nanoparticles. (Bottom) Storage modulus renormalized by its value at 10% volume fraction, Ḡ′ = G(ω, φ)/G(ω, 0.1), for various frequencies.

FIG. 7 :

 7 FIG.7: Effect of the aggregate size on the storage modulus of rigid aggregate PNCs. The case of nanoparticle is also shown (circle). The volume fraction is φ = 10%.

FIG. 8 :

 8 FIG. 8: Storage modulus as a function of the volume fraction, for frequency ω = 1 (top) and ω = 0.01 (bottom). We compare flexible aggregates (lozenge) and rigid aggregates (squares) for size N = 20, 30, 40, 50. The case of individual nanoparticles is also shown (circles).

FIG. 9 :

 9 FIG.9: Comparison between the monodisperse and the polydisperse rigid aggregates systems, see text for details. The volume fraction is φ = 10%.

FIG. 10 :

 10 FIG. 10: Dynamic structure factors computed for individual nanoparticles, flexible and rigid aggregates, with aggregate size N = 50 and volume fraction φ = 10%. From left to right, shown are the standard, self and intra-aggregate dynamic structure factors, as defined by Eqs. (23),[START_REF] Moll | Mechanical reinforcement in polymer melts filled with polymer grafted nanoparticles[END_REF], and[START_REF] Thomin | Network Effects on the Nonlinear Rheology of Polymer Nanocomposites[END_REF], and computed for wave vector q = 2π and q = 2π/3 ≃ 2π/Rg, where Rg is the aggregate radius of gyration.

FIG. 11 :

 11 FIG. 11: Correlation of the orientation vector for flexible and rigid aggregates of size N = 20 and 50. The volume fraction is φ = 10%.

  Each particle is linked by three virtual springs to particles randomly chosen among those that are are not neighbors. The virtual springs have the same spring constant k and their rest length is individually fixed by the structure of the aggregate[59]. As long as the spring constant is sufficiently high, the rigidity of the spring network will ensure that the aggregate structure is maintained, and represents a good approximation to the solid aggregates used in real nanocomposites.

FIG.

1

: Schematic (bidimensional projection) model of our aggregates. A group of repulsive particles are linked by two spring networks. The connective springs (black solid lines) between neighbors ensure the connectivity. The virtual springs (red dashed lines) maintain the geometry. For clarity, the latter are shown for only one particle.

TABLE I :

 I Choice of units and parameters for the typical polymer nanocomposite considered in simulation.

	Quantity	Symbol Simulation Typical value
	Energy	kBT	1	4.1 10 -21 J
	Nanoparticle diameter	σ	1	9.4 nm
	Polymer terminal time	τ 0 p	1	1 ms
	Polymer plateau modulus G 0 p	100 Pu	0.5 MPa

  For rigid aggregates, the modulus increases with aggregate loading and size. Both dependence are approximately linear over the range considered. As a result, the reinforcement is left unchanged if aggregates are twice larger but half less numerous (compare N = 20 and 40, at φ = 16% and 8 for instance). As a side note, one can see that somewhat unexpectedly, in flexible aggregates, reinforcement at ω = 1 is almost independent of size. At lower frequency ω = 0.01 (bottom), a differ-ent scenario appears. First, the storage modulus of rigid aggregates no longer increases linearly with volume fraction. Rather, we observe a superlinear behavior which can be rationalized under the form : G ′ (φ) ∝ φ m where m is an exponent which depends slightly on the aggregate size and takes values between 1.8 and 2.5 for the cases shown here. The behavior of flexible aggregates and individual nanoparticles is again clearly different from rigid aggregates. Both systems are characterized by a linear dependence of G ′ with volume fraction, and lower levels of reinforcement as compared to rigid systems. For flexible aggregates, a small size dependence appears. Overall, the level of reinforcement depends sensitively on aggregate type, size and content and on frequency regime. It is also clearly evident that aggregate rigidity is key to observing superlinear behavior, as observed experimentally in the same volume fraction range[START_REF] Zhu | Investigating Linear and Nonlinear Viscoelastic Behavior Using Model Silica-Particle-Filled Polybutadiene[END_REF].