Dafne Musino

Anne-Caroline Genix
email: anne-caroline.genix@umontpellier.fr

Caroline Fayolle

Aurelie Papon

Laurent Guy

Natalia Meissner

Radoslaw Kozak

Pawel Weda

Thomas Bizien

Thomas Chaussee

Aurélie Papon

Radosław Kozak

Thomas Chaussée

Julian Oberdisse

Synergistic Effect of Small Molecules on Large-Scale Structure of Simplified Industrial Nanocomposites

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The composition of modern car tire tread materials is of great complexity, and little is known about the impact of each of the many ingredients on the final structure and performance of such polymer nanocomposites (PNC) [START_REF] Mark | Science and Technology of Rubber[END_REF] , in spite of many fundamental studies of such materials [START_REF] Jancar | Current issues in research on structure-property relationships in polymer nanocomposites[END_REF][START_REF] Schmidt | Properties of polymer-nanoparticle composites[END_REF][START_REF] Kumar | Nanocomposites with Polymer Grafted Nanoparticles[END_REF][START_REF] Jordan | Experimental trends in polymer nanocomposites -a review[END_REF][START_REF] Kim | SBR/Organoclay Nanocomposites for the Application on Tire Tread Compounds[END_REF][START_REF] Chevigny | Tuning the Mechanical Properties in Model Nanocomposites: Influence of the Polymer-Filler Interfacial Interactions[END_REF] . For the industrial application, performance is usually quantified in terms of specific mechanical properties related to, e.g., rolling resistance or wet grip. [START_REF] Heinrich | Wet skid properties of filled rubbers and the rubber-glass transition[END_REF][START_REF] Manoharan | Exploring a highly dispersible silica-elastomer composite for tire applications[END_REF][START_REF] Sarkawi | A Review on Reinforcement of Natural Rubber by Silica Fillers for Use in Low-Rolling Resistance Tyres[END_REF] Since the advent of highly dispersible silica some twenty years ago, performance of such materials has been boosted [START_REF] Wolff | Chemical Aspects of Rubber Reinforcement by Fillers[END_REF][START_REF] Frohlich | The effect of filler-filler and filler-elastomer interaction on rubber reinforcement[END_REF][START_REF] Hasse | Influence of the amount of diand polysulfane silanes on the crosslinking density of silica filled rubber compounds. Kautsch. Gummi[END_REF][START_REF] Heinrich | Why silica technology needs S-SBR in high performance tires? The physics of confined polymers in filled rubbers. Kautsch. Gummi[END_REF] . Silanol surface groups on the silica result in a strong filler-filler interaction though at the detriment of the filler-polymer interaction. [START_REF] Bachmann | Fine Particle Reinforcing Silicas and Silicates in Elastomers[END_REF][START_REF] Wang | Filler-Elastomer Interactions. Part I: Silica Surface Energies and Interactions with Model Compounds[END_REF] The latter is dominant, e.g., for the traditional car tire filler -carbon black -, which display a high dispersive contribution of the surface energy. To overcome this disadvantage of the silica, which is due to the mismatch in polarity between the hydrophilic nanoparticle surfaces and the hydrophobic polymer, an obvious solution was found with the addition of coating and coupling agents. The firsts are small amphiphilic silane molecules playing a role analogous to surfactant molecules in, e.g., microemulsions. They are known to decrease both dispersive and polar contributions of the surface energy, whereas the introduction of a coupling agent compensates the decrease of the dispersive part, as compared to carbon black, by providing covalent links of the filler nanoparticles to the polymer matrix. [START_REF] Wolff | Filler-Elastomer Interactions. Part IV. The Effect of the Surface Energies of Fillers on Elastomer Reinforcement[END_REF] These agents are commonly used in modern car tires, and numerous studies have evidenced their impact on dispersion and macroscopic properties, see e.g., [START_REF] Sternstein | Reinforcement Mechanism of Nanofilled Polymer Melts As Elucidated by Nonlinear Viscoelastic Behavior[END_REF][START_REF] Bouty | Nanofiller Structure and Reinforcement in Model Silica/Rubber Composites: A Quantitative Correlation Driven by Interfacial Agents[END_REF][START_REF] Yan | Effect of silane coupling agent on the polymer-filler interaction and mechanical properties of silica-filled NR[END_REF][START_REF] Wu | The influence of in situ modification of silica on filler network and dynamic mechanical properties of silica-filled solution styrene-butadiene rubber[END_REF] . Note that in some cases, different surface chemistries, for instance with phosphonic acids, have been developed in order to adapt surface modification strategies to, e.g., aqueous environments [START_REF] Schmitt Pauly | Surface modification of alumina-coated silica nanoparticles in aqueous sols with phosphonic acids and impact on nanoparticle interactions[END_REF] , before casting the nanocomposite samples, and its effect on the final nanocomposite structure [START_REF] Schmitt Pauly | Structure of alumina-silica nanoparticles grafted with alkylphosphonic acids in poly(ethylacrylate) nanocomposites[END_REF] .

Several studies have been dedicated to the impact of (in situ) filler-surface modification with silanes on the mechanical properties of filled styrene-butadiene (SB) prepared by solid-phase mixing, sometimes reporting impact on microstructure, sometimes not. [START_REF] Frohlich | The effect of filler-filler and filler-elastomer interaction on rubber reinforcement[END_REF][START_REF] Ramier | Payne effect in silica-filled styrene-butadiene rubber: Influence of surface treatment[END_REF][START_REF] Stockelhuber | Impact of Filler Surface Modification on Large Scale Mechanics of Styrene Butadiene/Silica Rubber Composites[END_REF][START_REF] Yatsuyanagi | Effects of surface chemistry of silica particles on the mechanical properties of silica filled styrene-butadiene rubber systems[END_REF][START_REF] Suzuki | Effects of rubber/filler interactions on deformation behavior of silica filled SBR systems[END_REF][START_REF] Gauthier | Analysis of the non-linear viscoelastic behaviour of silica filled styrene butadiene rubber[END_REF][START_REF] Luginsland | Influence of different silanes on the reinforcement of silica-filled rubber compounds[END_REF][START_REF] Mele | Reinforcement effects in fractal-structure-filled rubber[END_REF] Such systems are close to the industrial application, usually with a crosslinked matrix (rubber, SBR), involving the addition of various additives like, e.g., zinc oxide, sulfur, CBS (cyclohexyl benzothiazole sulfenamide) and DPG (diphenyl guanidine), referred to as 'curing package' [START_REF] Stockelhuber | Impact of Filler Surface Modification on Large Scale Mechanics of Styrene Butadiene/Silica Rubber Composites[END_REF] . Stöckelhuber et al [START_REF] Stockelhuber | Impact of Filler Surface Modification on Large Scale Mechanics of Styrene Butadiene/Silica Rubber Composites[END_REF] have studied SBR filled with different silica nanoparticles and silane grafting, and provided evidence by dynamical mechanical measurements for the impact of the interfacial properties on large-scale mechanical properties and in particular the Payne effect [START_REF] Payne | The dynamic properties of carbon black-loaded natural rubber vulcanizates. Part I[END_REF] . The Payne effect in silica-filled SBR was also studied as a function of surface modification with coating or coupling agents by Ramier et al [START_REF] Ramier | Payne effect in silica-filled styrene-butadiene rubber: Influence of surface treatment[END_REF] . In this case, the dispersion state of the silica was not found to vary, whereas the amplitude of the Payne effect was found to decrease with surface modification, as also observed by Gauthier et al [START_REF] Gauthier | Analysis of the non-linear viscoelastic behaviour of silica filled styrene butadiene rubber[END_REF] . The reversible decrease of elastic modulus associated with the Payne effect was correlated with bound rubber measurements and discussed in terms of chain desorption from the filler surface. Several groups have observed correlations of microscopic structure with mechanical properties. [START_REF] Wu | The influence of in situ modification of silica on filler network and dynamic mechanical properties of silica-filled solution styrene-butadiene rubber[END_REF][START_REF] Yatsuyanagi | Effects of surface chemistry of silica particles on the mechanical properties of silica filled styrene-butadiene rubber systems[END_REF][START_REF] Pérez-Aparicio | Reinforcement in Natural Rubber Elastomer Nanocomposites: Breakdown of Entropic Elasticity[END_REF][START_REF] Conzatti | Morphology and viscoelastic behaviour of a silica filled styrene/butadiene random copolymer[END_REF] Wu et al have studied the influence of surface modification of silica filler with a coupling agent in SBR. [START_REF] Wu | The influence of in situ modification of silica on filler network and dynamic mechanical properties of silica-filled solution styrene-butadiene rubber[END_REF] They have shown that an in situ modification of silica NPs with coupling agent increases the content of bound rubber, and improves the dispersion of silica, the stability of which is ensured by the chemical bonds between filler and rubber. They also provided rheological evidence for a simultaneous improvement of both rolling resistance and wet grip, i.e., the industrially relevant quantities. Conzatti et al have studied the impact of surface modification of silica in SBR on structure by TEM. [START_REF] Conzatti | Morphology and viscoelastic behaviour of a silica filled styrene/butadiene random copolymer[END_REF][START_REF] Castellano | Influence of the silane modifiers on the surface thermodynamic characteristics and dispersion of the silica into elastomer compounds[END_REF] They have shown that the better dispersion obtained with silane agents is correlated with a shift in percolation threshold towards higher filler contents, as well as a higher tensile strength at high deformation. Finally, the use of a coating agent to reduce filler-filler interactions and improve dispersion in SBR was exploited by Mujtaba et al. [START_REF] Mujtaba | Mechanical Properties and Cross-Link Density of Styrene-Butadiene Model Composites Containing Fillers with Bimodal Particle Size Distribution[END_REF] This induces a higher specific surface which contributes strongly to mechanical properties, whereas an additional fraction of micron-sized fillers is found to have little effect. In all these examples, the SBR nanocomposites display a complex formulation -including curing agents as mentioned above -, and the grafting of coating or coupling agents onto the silica surfaces cannot be studied independently.

On the contrary, the impact of individual curing ingredients on rheology has been investigated by Schaal et al. [START_REF] Schaal | The Effects of Certain Recipe Ingredients and Mixing Sequence on the Rheology and Processability of Silica-and Carbon Black-Filled Tire Compounds[END_REF] They followed rheological changes upon annealing and found that DPG leads to a retardation of the presumed re-agglomeration of silica particles. DPG is a basic molecule traditionally used as a secondary accelerator, but it was also found to improve the silanization reaction during mixing. [START_REF] Penot | Rubber composition for a tire comprising a reinforcing inorganic filler and an (inorganic filler/elastomer) coupling system[END_REF] This is beneficial for the curing process as higher coverage of the silica surface allows better efficiency of the accelerator (e.g., CBS), which otherwise tends to adsorb on the silica surface. [START_REF] Kaewsakul | Optimization of mixing conditions for silica-reinforced natural rubber tire tread compounds[END_REF] It follows that the mixing order of the ingredients is of primary importance. If, in many cases, DPG was added together with the curing package -i.e., at the end of the mixing process [START_REF] Yatsuyanagi | Effects of surface chemistry of silica particles on the mechanical properties of silica filled styrene-butadiene rubber systems[END_REF][START_REF] Suzuki | Effects of rubber/filler interactions on deformation behavior of silica filled SBR systems[END_REF][START_REF] Mujtaba | Mechanical Properties and Cross-Link Density of Styrene-Butadiene Model Composites Containing Fillers with Bimodal Particle Size Distribution[END_REF][START_REF] Choi | Difference in bound rubber formation of silica and carbon black with styrenebutadiene rubber[END_REF] or later by using a two-roll mill [START_REF] Ramier | Payne effect in silica-filled styrene-butadiene rubber: Influence of surface treatment[END_REF][START_REF] Schaal | The Effects of Certain Recipe Ingredients and Mixing Sequence on the Rheology and Processability of Silica-and Carbon Black-Filled Tire Compounds[END_REF][START_REF] Castellano | Surface modification of silica: 1. Thermodynamic aspects and effect on elastomer reinforcement[END_REF] -, it was recently introduced in the first stage of mixing together with the coating [START_REF] Baeza | Multiscale Filler Structure in Simplified Industrial Nanocomposite Silica/SBR Systems Studied by SAXS and TEM[END_REF] or coupling [START_REF] Bouty | Nanofiller Structure and Reinforcement in Model Silica/Rubber Composites: A Quantitative Correlation Driven by Interfacial Agents[END_REF] agent and the silica, following a patented innovation by Michelin 37 . To the best of our knowledge, however, no quantitative study of filler dispersion and aggregation in styrenebutadiene of the separate and combined impact of both coating agent and DPG on the final PNC structure exists. This is the objective of the present letter.

In recent contributions, we have developed a combined TEM-SAXS-simulation analysis of nanoparticle and aggregate dispersions in systems of disordered silica in styrene-butadiene copolymer of industrial relevance but with a limited number of ingredients [START_REF] Baeza | Multiscale Filler Structure in Simplified Industrial Nanocomposite Silica/SBR Systems Studied by SAXS and TEM[END_REF][START_REF] Baeza | Effect of Grafting on Rheology and Structure of a Simplified Industrial Nanocomposite Silica/SBR[END_REF][START_REF] Baeza | Studying twin-samples provides evidence for a unique structure-determining parameter in simplified industrial nanocomposites[END_REF][START_REF] Genix | Recent advances in structural and dynamical properties of simplified industrial nanocomposites[END_REF] . Motivated, in particular by the absence of any additive related to curing but DPG, the system with un-crosslinked SB matrix was termed 'simplified industrial'. The objective of the present letter is to use our quantitative analysis to characterize the structure of the filler particles in such systems as a function of the concentration of small molecules, i.e., coating agent and DPG, in order to investigate their -empirically known to be synergetic -impact on nanocomposite structure.

Materials and methods

Sample formulation. Styrene-butadiene random copolymer was purpose-synthesized by Synthos in two highly monodisperse versions (polydispersity index = 1.02): one non functionalized of molecular weight 177 kg/mol (denoted 'NF'), the other (170 kg/mol, denoted 'F') with a single graftable endfunction SiMe 2 -OH, also called D3. The synthesis protocol is outlined elsewhere. [START_REF] Baeza | Effect of Grafting on Rheology and Structure of a Simplified Industrial Nanocomposite Silica/SBR[END_REF] An industrially relevant example of a coating agent is octyltriethoxysilane ("octeo"), which is typically used at 8%wt with respect to silica, and the content of which is varied in the present study from zero to the latter value (0, 4%, 8%). The mass of the short aliphatic chain of the silane was also varied, from octyl-to dodecyl-and octadecyl-triethoxysilane (termed C 8 from Sigma-Aldrich, C 12 and C 18 from Alfa Aesar, of molar mass 276, 332, and 417 g/mol, resp.), keeping the molar concentrations constant with respect to the 8%wt octeo samples. The filler is obtained by fragmentation of precipitated silica micropellets (Zeosil(R) 1165MP from Solvay), which have a nominal specific surface of 160 m 2 /g. The size distribution of the primary nanoparticles obeys a log-normal law in agreement with TEM studies (R 0 = 8.55 nm, σ = 27%, leading to an average bead volume of V Si = 3.6 10 3 nm 3 ). [START_REF] Baeza | Multiscale Filler Structure in Simplified Industrial Nanocomposite Silica/SBR Systems Studied by SAXS and TEM[END_REF] Finally, 1,3-diphenyl guanidine (DPG from Bayer, 211 g/mol, 1%wt with respect to polymer) was also used during the mixing process. As discussed above, its empirical role is to enhance silanization. It is the second small molecule the impact of which is investigated here, on its own and in combination with octeo to characterize the potential synergy. PNC samples were produced by stepwise introduction and mixing in an internal mixer at a temperature reaching 155°C to ensure silanization [START_REF] Reuvekamp | Effects of Time and Temperature on the Reaction of Tespt Silane Coupling Agent During Mixing with Silica Filler and Tire Rubber[END_REF] . The polymer was introduced first at 70°C and masticated for 1 min before introduction in two equal steps of the silica together with the coating agent and DPG. Mixing is followed by two passages between rotating cylinders at 70°C first, then another twelve passages after cooling to room temperature. Finally, samples have been pressed at T = 60°C, P = 60 bar for 3 min, to reach a final sample thickness of about 1.5 mm. Effective grafting after the mixing process was checked by bound rubber (BR) measurements giving the quantity of non-extractable chains. The BR fraction has been determined by extraction of the free chains in xylene and measuring the resulting sample mass after evaporation of the solvent. This indicator for the efficiency of chain grafting on the silica surface is also linked to the structuration of the filler [START_REF] Baeza | Effect of Grafting on Rheology and Structure of a Simplified Industrial Nanocomposite Silica/SBR[END_REF] . Note that due to the absence of antioxidant agent in the formulation, some weak crosslinking of SB occurs during mixing resulting in an incomplete dissolution of the matrix polymer in a good solvent of the polymer (THF). For comparison, some crosslinked samples have also been produced by adding sulfur and CBS (1.5%wt with respect to polymer for each) after the internal mixing step. In order to reduce the number of parameters, a single silica volume fraction of Φ Si ≈ 18%v has been investigated here. For all samples, the volume fractions characterized by thermogravimetric analysis display some small spread of no more than 0.3%v.

Structural analysis.

Small-angle X-ray scattering was performed on beamline SWING at synchrotron SOLEIL (Saint Aubin, France) using standard conditions (sample-to-detector distances 2 m, 5 m and 6.5 m; wavelength 1 Å, giving a q-range from 6.2 10 -4 Å -1 to 5.6 10 -1 Å -1 ). Standard data reduction tools given by Soleil were used (Foxtrot 3.1). Matrix contributions have been measured independently, and subtracted. Transmission electron microscopy (TEM) was performed on samples prepared by ultracryomicrotomy at -120 °C (desired thickness set to 90 nm, copper grids with carbon membrane) using a FEI Tecnai Biotwin at 120 kV. SAXS reveals a multi-scale filler structure: NPs forms aggregates which are themselves concentrated within fractal branches extending over micrometers. In the quantitative structural model that we have developed [START_REF] Baeza | Multiscale Filler Structure in Simplified Industrial Nanocomposite Silica/SBR Systems Studied by SAXS and TEM[END_REF] , the central element is the aggregate defined by its internal characteristics, namely average compacity κ, and aggregation number N agg . The SAXS analysis is performed in two steps: (i) the use of a Kratky plot (q²I(q) versus q) to extract the average aggregate radius R agg , and (ii) the description of the inter-aggregate interaction by a simulated structure factor for polydisperse spheres representing aggregates with hard sphere interactions. As the silica fraction is fixed here, the latter contribution has a minor impact, and changes in the SAXS data are mainly related to a different aggregate form factor (mass and size). Comparison with the SAXS intensity at a low-q value leads to the determination of the average aggregate compacity. [START_REF] Schmitt Pauly | Structure of alumina-silica nanoparticles grafted with alkylphosphonic acids in poly(ethylacrylate) nanocomposites[END_REF][START_REF] Baeza | Multiscale Filler Structure in Simplified Industrial Nanocomposite Silica/SBR Systems Studied by SAXS and TEM[END_REF] Note that quantitative input of TEM in terms of the volume fraction of fractal branches Φ fract (or equivalently the volume fraction of pure polymer zones, given by 1 -Φ fract ) is needed to estimate the internal volume fraction, Φ agg = Φ Si /(κ Φ fract ). [START_REF] Baeza | Multiscale Filler Structure in Simplified Industrial Nanocomposite Silica/SBR Systems Studied by SAXS and TEM[END_REF] On average for the samples with functionalized chains, Φ fract about 80±2% was found. For NF and 50/50, we used Φ fract = 58 and 73%, respectively [START_REF] Baeza | Effect of Grafting on Rheology and Structure of a Simplified Industrial Nanocomposite Silica/SBR[END_REF] .

Results and discussion

Impact of matrix composition on nanocomposite structure. In recent studies of simplifiedindustrial nanocomposites, we have investigated the impact of chain parameters on the filler aggregate structure by SAXS, and in particular chain end-functionalization [START_REF] Baeza | Effect of Grafting on Rheology and Structure of a Simplified Industrial Nanocomposite Silica/SBR[END_REF] and chain mass [START_REF] Baeza | Mechanism of aggregate formation in simplified industrial silica styrene-butadiene nanocomposites: effect of chain mass and grafting on rheology and structure[END_REF] . In order to test the influence of the coating agent and DPG reported in the present letter, a matrix composition had to be chosen. As the chains had been newly synthesized with a slightly different chain mass than before 41, 44 , the impact of end-functionalization was studied first, in order to choose the most promising system in terms of structural evolution. For the same reason, we have chosen to work with a high volume fraction of silica filler (Φ Si = 18%v). In Figure 1, three normalized SAXS intensity curves are shown, corresponding to either a pure non functionalized matrix (NF), a fully functionalized one (F), and a 50/50 mixture. The structural features of these nanocomposites are clearly recognized in the scattered intensities: the high-q scattering is reminiscent of the primary silica NPs, identical in all formulations. A small bump in the range 0.02 -0.03 Å -1 (well visible in a Kratky representation, see later in Figure 3b) characterizes the interactions between particles in contact, and below that q-value aggregate scattering is dominant. A much stronger change in slope around 0.01 Å -1 is indicative of aggregate-aggregate interactions, and this position in q-space is used to extract the average aggregate radius R agg from the scattering [START_REF] Baeza | Multiscale Filler Structure in Simplified Industrial Nanocomposite Silica/SBR Systems Studied by SAXS and TEM[END_REF] . In this intermediate q-range, the strong influence of the matrix composition on the structure is visible. At low-q, finally, a low-q upturn typical for large-scale filler structures is visible. A quantitative analysis confirming the impact of matrix composition has been performed following ref. [START_REF] Baeza | Multiscale Filler Structure in Simplified Industrial Nanocomposite Silica/SBR Systems Studied by SAXS and TEM[END_REF] . It allows the determination of not only R agg , but also the aggregate density (compacity κ), and thus the average aggregation number N agg , all given in Table 1. As in our previous work with a different polymer [START_REF] Baeza | Effect of Grafting on Rheology and Structure of a Simplified Industrial Nanocomposite Silica/SBR[END_REF] , the more functionalized matrix chains are, the lower the aggregate radius determined through the Kratky representation, and the lower the κ. This tendency results in a decrease of the aggregation number by a factor of two. The evolution towards smaller aggregates is accompanied by the formation of a well-defined shoulder around 0.01 Å -1 : nanocomposites formed with the end-functionalized polymer chains thus present the most prominent signature in I(q), indicative of the best dispersion of smallest possible aggregates. Moreover, the effect of the graftable chains can also be evidenced by bound rubber measurements, the results of which are also summarized in Table 1 and in the inset in Figure 1. Even with only un-functionalized polymer (NF), about half of the chains are not desorbed in a BR measurement, most probably due to weak crosslinking of the chains during mixing, or possibly chain adsorption on the filler. The high fraction of bound rubber with graftable chains indicates that grafting took place. We thus choose to study the impact of the small molecules, silane and DPG, only with fully functionalized, graftable matrices, by investigating how the dispersion evolves as the quantity of small molecules is reduced. Impact of hydrophobic character of coating agent on nanocomposite structure: Subsequent to the choice of the polymer matrix, the coating agent needed to be chosen. We have varied the chemical nature of the coating molecule, via the length of its aliphatic chain. The macroscopic bound rubber values are almost identical (96 -99%) and are given in Table 1. This suggests that the chain length of the coating agent has only little influence on polymer grafting, which we have seen in Figure 1 to be an important driver of the PNC microstructure. The scattering results for an 18%-nanocomposite, in a fully functionalized polymer matrix, with DPG and isomolar amounts of C 8 (octeo, 8%wt with respect to silica), C 12 , or C 18 , are shown in Figure 2. Contrarily to what one might expect based on the formation of a more hydrophobic layer on the filler surface, potentially inducing a better compatibilization (and thus dispersion), there is no significant effect of the nature of the coating agent on the aggregate structure. Indeed, the three intensity curves superimpose almost perfectly. This surprising result may be due to a saturation effect, i.e., adding the same number of higher-mass agents would not increase the silica surface coverage. In this picture, C 8 molecules already reach a perfect coverage, and C 18 only modifies the thickness of the interfacial layer. Also, bound rubber measurements seem to show no impact on the efficiency of polymer grafting on the silica surface, further indicating similar surface coverage. In the rest of this article, the surface modification is thus performed with the C 8 coating agent (octeo), which incidentally is also the industrially relevant one.

Impact of amount of coating agent on nanocomposite structure. In Figure 3a, the normalized small-angle intensities scattered by nanocomposites of various octeo contents used during the formulation, 0, 4, and 8%wt with respect to the filler, respectively, are shown, all other parameters being fixed, namely the filler content and DPG. Arrows indicate the q-ranges of particle-and aggregate-scattering, respectively. b) Kratky plot q 2 I(q) vs q highlighting the structural evolution (dashed arrow). Inset: All BR for functionalized polymers versus the specific surface of aggregates.

The scattered intensities in Figure 3a are characteristic of the same type of structure as discussed with Figure 1: obvious changes are induced by the presence of the coating agent, inducing better filler dispersion. A similar result has been reported in the literature for an ungrafted polymer matrix. [START_REF] Bouty | Nanofiller Structure and Reinforcement in Model Silica/Rubber Composites: A Quantitative Correlation Driven by Interfacial Agents[END_REF] In parallel, the fraction of bound rubber shown in the inset is seen to increase with the amount of octeo.

Figure 3a is a key result of our study, as it demonstrates, to our knowledge for the first time, the delicate impact of the coating agent on structure in graftable matrices. The results of our quantitative analysis are given in Table 1. It reveals a decrease by more than 10% in the average aggregate radius R agg , which is also seen in the Kratky representation in Figure 3b, with the shift of the maximum to higher q illustrated by the dashed arrow. This maximum is the direct consequence of the break in slope in I(q), indicative of the aggregate interaction, highlighted by the arrow around 0.01 Å -1 . From the analysis of the height of the signal in the 'aggregate region' up to 0.01 Å -1 , using the volume fraction of fractal branches determined independently from TEM pictures shown in Figure 4 (Φ fract ≈ 80%), the compacity is also found to decrease, as given in Table 1. As a result of the combination of decrease both in size and density, the average aggregation number N agg is seen to diminish by 35%, from 48 to 31. In the inset of Figure 3a, the bound rubber results are seen to slightly increase with increasing octeo content. As this coating agent does not contribute to any grafting of the polymer, this is a priori surprising. Our interpretation is that the decrease of N agg favors polymer-filler interactions, via a higher specific surface S/V estimated by 3*Φ agg /R agg in the inset of Figure 3b. Indeed, both bound rubber and S/V increase by some 10%. It cannot be decided, however, if the higher bound rubber is only adsorbed to the aggregates, possibly entangled with grafted chains making it difficult to desorb, or if the higher S/V favors polymer access to surface and grafting. In parallel, the TEM images of the octeo-series in Figure 4 illustrate the progressive formation of smaller and smaller structures with increasing octeo content. Homogeneity on larger scales has been checked with lower resolution pictures. As expected, compatibilization of the originally hydrophilic filler surface with the hydrophobic polymer matrix through the coating agent induces the formation of smaller aggregates, i.e., a better -more "nano" -dispersion. and DPG, or only with octeo, than in presence of only DPG or both agents. For a quantitative analysis, however, the evolutions are too weak, and SAXS is needed. In Figure 6a, the normalized scattered intensities are plotted for these samples, and the Kratky plots in Figure 6b highlight the structural evolution. All the structures characterized by scattering in Figure 6a are again typical for this type of simplified industrial nanocomposites. The coating agent and DPG have a pronounced impact on the aggregate structure probed at intermediate q, as one can see from the different intensity levels -keeping in mind the logarithmic scale: without any DPG or coating agent, the dispersion is worst, with the highest intensity in Figure 6a, while in presence of both -cf. the lowest, red curve in the figure -, the best dispersion is obtained. This is confirmed by the quantitative analysis of the Kratky plot (Figure 6b) presenting a shift to the right with added molecules, which leads to the average aggregate radius given in Table 1: R agg decreases from 38 (no agents) to 30 nm (both agents). Simultaneously, the compacity and the aggregation number evolve from к = 48% to 44%, and N agg from 67 to 31, respectively, while the bound rubber fraction increases from below 70% to above 95%. The decrease in N agg is indicative of a better compatibility of the filler with the matrix, as there is more polymer in contact with fillers for smaller aggregates, and this leads to the highest bound rubber value in Figure 6a.

One might wonder if octeo and DPG produce a synergetic effect on the microstructure, or if only one of them might be effective. The two intermediate samples with either one of them have also been analyzed; cf. Figure 6 and Table 1. Adding only octeo decreases the aggregate radius to 36 nm and the aggregation number N agg to 61, whereas the bound rubber fraction increases. The latter behavior is consistent with the variation observed for different octeo contents in presence of DPG, cf. Figure 3, and probably has the same causes -a higher specific surface. Finally, adding only DPG has an unexpected effect on the structure, as it leads to a better dispersion than without any agent: the average aggregate radius is about 34 nm, while compacity is 47% and the aggregation number N agg ≈ 48. Such results may be due to the strong affinity of the basic DPG with the silica surface, which displays acidic character. DPG therefore preferentially adsorbs on the (untreated) silica, preventing by this way the formation of hydrogen bonds between hydroxyl groups of different nanoparticles. [START_REF] Zaborski | Activity of fillers in elastomer networks of different structure[END_REF][START_REF] Kosmalska | Adsorption of curatives and activity of silica toward elastomers[END_REF] Note that the number density of silanol groups on silica has been shown to be correlated with the aggregate size in SBR nanocomposites. [START_REF] Yatsuyanagi | Effects of secondary structure of fillers on the mechanical properties of silica filled rubber systems[END_REF] Comparison between the structural results with DPG only and octeo only further indicates that the efficiency of the coating agent alone is not optimal. The addition of DPG together with octeo is required to obtain the better dispersion (see Figure 6 and Table 1), DPG promoting the grafting reaction of the silane [START_REF] Penot | Rubber composition for a tire comprising a reinforcing inorganic filler and an (inorganic filler/elastomer) coupling system[END_REF] . On the other hand, one may wonder if DPG also has a catalytic action on the grafting reaction of the polymer chains on the silica nanoparticles (all chains are expected to carry a Si-OH end-function). This does not seem to be the case as the BR values for DPG only and octeo only are identical. It follows that the 40%-increase of the bound rubber fraction from no agents to both agents in Figure 6a is mostly driven by the structural evolution: the decrease of N agg leads to a higher specific surface (≈ +45%), favoring chain grafting. Confronting these findings with presence and absence of both molecules unambiguously proves that DPG and octeo are necessary to tune the large-scale filler structure efficiently.

In order to extend our discussion towards the influence of small molecules on crosslinked samples as they may be used in applications, we have first investigated the vulcanization kinetics given by the curing curves (torque versus time) of the pure matrix with and without DPG, see inset of Figure 7. As observed by Ramier et al. [START_REF] Ramier | Influence of Silica and its Different Surface Treatments on the Vulcanization Process of Silica Filled SBR[END_REF] , using only CBS (without DPG) shifts the beginning of the vulcanization process to longer times, which means that both accelerators are required for a rapid cure. In the meantime, we found that the presence of DPG impacts the final crosslinking density only by a few percent. We have then formulated four crosslinked samples in analogy to the four samples compared in Figure 6. The experimental difficulty with such samples is that the curing agents are added in the two-roll mill process, inducing slight anisotropy in the samples and thus intensities. Luckily, the average level of these intensities as shown in the next plot, however, varies strongly with the composition, making conclusions robust independently of anisotropy. In Figure 7, the same color code as for Figure 6 applies, comparing crosslinked samples with both DPG and octeo, with only one, or neither. The scattering of the corresponding uncrosslinked samples (Fig. 6a) is superimposed in grey for convenience. In Figure 7, crosslinked and uncrosslinked structures are identical at high q, and at low q intensities tend to converge. In the intermediate q-range, all four intensity curves of crosslinked samples are higher by less than a factor of two than the uncrosslinked ones. This means crosslinking favors slightly the aggregation of silica nanoparticles, similarly to previous results [START_REF] Baeza | Revealing nanocomposite filler structures by swelling and small-angle X-ray scattering[END_REF] . For instance, comparing the intensity levels of the samples having both small molecules, the crosslinked one has bigger aggregates R agg (34 nm, uncrosslinked: 32 nm), about the same compacity (к = 47%, resp. 44%), and a higher N agg = 46 (resp. 31). It is unclear, however, if this increase in aggregation is due to the crosslinker which might influence interactions between particles, or if mechanical action during crosslinking has this effect. The most remarkable fact in Figure 7 is that the order of the curves is unchanged with respect to the uncrosslinked case. Our findings that both small molecules, DPG and octeo, are necessary to tune the large-scale filler structure efficiently thus remain valid also for crosslinked samples.

Conclusion

We have investigated the impact of small molecules -coating agents and DPG added during formulation by solid-phase mixing -on the bound rubber fraction and on the microstructure of polymer nanocomposites made with disordered silica filler of industrial relevance. The bound rubber is found to be correlated with the microstructure (see inset of Figure 3b), suggesting a specific-surface effect, smaller aggregates offering more possibilities to bind polymer. The microstructure is characterized by small-angle scattering and TEM, with a particular focus on aggregates of primary NPs, which we have shown to be the fundamental structural units, cf. recent reviews [START_REF] Genix | Recent advances in structural and dynamical properties of simplified industrial nanocomposites[END_REF][START_REF] Genix | Structure and dynamics of polymer nanocomposites studied by Xray and neutron scattering techniques[END_REF] .

Accordingly, the structure can be quantitatively described with two independent measures, the average aggregate radius R agg , and the average aggregation number N agg , or equivalently, the compacity κ. Both have been shown to evolve with the coating agent content, higher concentrations of this additive giving better dispersion, i.e., smaller aggregates, presumably due to a better compatibility with the matrix. There was surprisingly no impact of the chain length of the coating agent on the nanocomposite microstructure, suggesting that short chains are sufficient to cover the nanoparticle surface. Concerning this latter point, one might nonetheless expect a modification of the local dynamics, caused by a plasticizing effect of the matrix by the longer aliphatic chains, a point which we hope to address in a future study. Along the same lines, one could also investigate if the absence of effect is due to saturation of the surface, e.g., by performing a study of longer molecules as a function of (lower) concentrations.

By confronting the dispersion of nanocomposites with and without the amine DPG, this molecule was unambiguously demonstrated to be necessary to promote the compatibilizing effect of the coating agent: there is a strong synergy between the amine and the coating agent. This conclusion was shown to remain true even for industrially more relevant crosslinked samples. By comparing bound rubber and structure to a sample containing only DPG, it was further shown that DPG does not catalyze the matrix end-grafting. Incidentally, this could be cross-checked in future studies with nanocomposites made only with a non graftable matrix, which was excluded here in order to obtain the strongest possible structural signature of synergy.

Figure 1 :

 1 Figure1: Reduced SAXS intensities I(q)/Φ allowing the structural comparison between nanocomposites containing 18.0±0.3%v of silica, coating agent octeo (8%wt), and DPG, made with different polymer matrices F (squares), NF (circles), and a 50/50 mixture (diamonds). Arrows indicate the q-ranges of nanoparticle-and aggregate-scattering, respectively. Inset: Bound rubber results for the same samples.

Figure 2 :

 2 Figure 2: Reduced SAXS intensities allowing the structural comparison between nanocomposites containing 18.0±0.3%v silica in a fully functionalized matrix with DPG, for three different isomolar silane coating agents: octeo (C 8 , squares), C 12 (circles), and C 18 (triangles). The inset shows the bound rubber of the same PNCs.

1 )Figure 3 :

 13 Figure3: a) Reduced SAXS intensities allowing the structural comparison between nanocomposites containing 18.0±0.3%v silica in a fully functionalized matrix with DPG, made with different amounts of coating agent octeo: 0 (diamonds), 4 (crosses), and 8%wt (squares). The inset shows the bound rubber of the same PNCs. Arrows indicate the q-ranges of particle-and aggregate-scattering, respectively. b) Kratky plot q 2 I(q) vs q highlighting the structural evolution (dashed arrow). Inset: All BR for functionalized polymers versus the specific surface of aggregates.

Figure 4 :

 4 Figure 4: TEM images of nanocomposites containing 18.0±0.3%v silica in a fully functionalized matrix with DPG, made with different amounts of coating agent octeo: 0, 4, and 8%wt with respect to silica.

Figure 5 :

 5 Figure 5: TEM images of nanocomposites containing 18.0±0.3%v silica in a fully functionalized matrix: (a) without small molecules, (b) with octeo only (8%wt with respect to silica, no DPG). Images with only DPG, and with both octeo and DPG are presented in Figure 4.

Figure 6 :

 6 Figure 6: a) Reduced SAXS intensities comparing the structure of nanocomposites containing 18.0±0.3%v silica in a fully functionalized matrix: without small molecules (bare, triangles), with DPG only (diamonds), with octeo only (circles), and with both octeo and DPG (squares). The fraction of octeo (when present) is 8%wt. The inset shows the bound rubber of the same PNCs. b) Kratky representation of the same samples.

Figure 7 :

 7 Figure 7: Reduced SAXS intensities comparing the structure of crosslinked nanocomposites containing 17.5±0.3%v silica in a fully functionalized matrix: without small molecules (bare, triangles), with DPG only (diamonds), with octeo only (circles), and with both octeo and DPG (squares). The fraction of octeo (when present) is 8%wt. The uncrosslinked curves shown in Figure 6 are superimposed in grey. Inset: Curing curves of unfilled polymer: SB+Sulfur+CBS, and SB+Sulfur+CBS+DPG.

Table 1 :

 1 Nanocomposite formulations with silica content Φ Si = 18.0 ± 0.3%v. Average radius of aggregates from Kratky analysis, compacity, aggregation number, and bound rubber of PNCs. All samples contain 8%wt C 8 with respect to silica (resp., isomolar C 12 , C 18 ) and DPG, unless otherwise stated.

	Sample	Coating	Φ Si (%v)	R agg (nm)	κ (%)	N agg	%BR
	NF	C 8	18.3	33	70	67	46
	50/50	C 8	18.2	33	51	47	57
	F	C 8	18.1	30	44	31	96
	F	C 12	18.1	30	45	31	98
	F	C 18	18.3	31	45	34	99
	F	4%wt C 8	17.8	32	46	39	92
	F	-	17.5	34	47	48	88
	F, no DPG	C 8	18.1	36	50	61	88
	F, no DPG	-	17.4	38	48	67	69

Acknowledgements:

The authors are thankful for support by the ANR NANODYN project, grant ANR-14-CE22-0001-01 of the French Agence Nationale de la Recherche. We thank Benoit Moreaux (Solvay Silica) for helping with the preparation of the samples. TEM preparation at Solvay was done by Pauline Grau and Magalie Fontana.

FOR TABLE OF CONTENTS USE ONLY

Synergistic effect of small molecules on large-scale structure of simplified industrial nanocomposites Dafne Musino, Anne-Caroline Genix, Caroline Fayolle, Aurélie Papon, Laurent Guy, Natalia Meissner, Radosław Kozak, Pawel Weda, Thomas Bizien, Thomas Chaussée, Julian Oberdisse


Daniel Borgis

Roland Assaraf

Benjamin Rotenberg

Rodolphe Vuilleumier

Computation of pair distribution functions and threedimensional densities with a reduced variance principle

Introduction

What we all learned from our mentor is that, for a given statistical mechanics problem, it is important (and intellectually rewarding) to first derive the proper, irrefutable, mathematical formulation, to understand as much as possible the formulas in mechanical terms, and, at the end, to be able to express the relevant formulas in computational forms that can be tackled with efficient and accurate simulation algorithms. Such goal can be hardly reached without the help of a chalk and a blackboard, and without questioning every single line. 1We hope that the present contribution is as much as possible in this vein, although it stays away from any complex problems in computational statistical mechanics that have attracted the community recently [START_REF] Ferrario | Computer Simulations in Condensed Matter: From Materials to Chemical Biology[END_REF][START_REF] Ferrario | Computer Simulations in Condensed Matter: From Materials to Chemical Biology[END_REF], in particular quantum-classical statistical mechanics [3,4], statistical mechanics with constraints [5][6][7] or phase-space exploration and characterization for complex systems using accelerated sampling methods [8,9].

We go back here to the computation of the most elementary objects in liquid-state theories, namely the two-particle pair-distribution functions in bulk liquids, or the three-dimensional one-particle liquid density around a tagged molecular object. [START_REF] Mcquarrie | Statistical Mechanics[END_REF][START_REF] Hansen | Theory of Simple Liquids[END_REF] The evaluation of such quantities, as prescribed in simulation textbooks [START_REF] Allen | Computer simulation of liquids[END_REF][START_REF] Frenkel | Understanding Molecular Simulations: From Algorithms to Applications[END_REF], is usually done through histograms, or possibly improved histograms, that rely in any case on the representation of the Dirac δ-function on a discrete grid. With such method, the variance of the results diverges when the grid spacing goes to zero. Elaborating on the work of Assaraf and collaborators [START_REF] Assaraf | [END_REF]15] for the evaluation quantum electronic densities with zero-variance/zero-bias principles, we propose a different approach to the computation of pair distribution functions and one-particle densities that leads to a much reduced variance of the results with respect to straight histogram counting and furthermore achieves a finite variance when the grid spacing tends to zero. As will be seen, smooth curves can be obtained even with very limited statistics.

Their will be a fair amount of explicit rigorous mathematical manipulations in this contribution, as it suits to a rigorous scientist who denies any physical (not to speak of chemical) intuition and thinks that Nature can only speak through well defined mathematical objects and their relations. We must confess, however, that some points of the derivation will not have the degree of rigor that it deserves; we will elude in particular some mathematical proofs of convergence and variance properties -in part because our degree of understanding and current formalization is not yet complete.

Nevertheless, all the basic principles will be set.

The paper is organized as follows. Section 2 is devoted to the computation of pair distribution functions whereas section 4 focuses on three-dimensional densities. Section 4 concludes.

Pair distribution functions

Let us consider a fluid mixture for which microscopic configurations are generated by computer simulations in the canonical ensemble. The pair distribution function (or radial distribution function, RDF) between molecules of type a and b is defined as

g ab (r) = ab 4πr 2 Na i=1 N , b j=1 δ(r -r ij ) , (1)
where ab = (1

-1 2 δ ab ) V NaN b , r ij = r j -r i , and r ij = |r ij |.
The prime in the second sum indicates that i = j should be excluded in the case a = b. This function is generally computed through histograms with bins of of finite width, ∆r, replacing the δ-function by 1 ∆r h ∆r (r), where h ∆r (r) is the characteristic function equal to 1 between r and r + ∆r and 0 otherwise. This statistical estimation is of infinite variance as ∆r → 0 since the instantaneous density in each bead oscillates between 0 and O(1/4πr 2 ∆r).

Accounting for the rotational invariance in the bulk fluid, the pair distribution can be also expressed as

g ab (r) = ab 4π dΩ Na i=1 N b j=1 δ(r -r ij ) , (2)
with Ω = r/r. Following the ideas of Assaraf et al for electron densities [START_REF] Assaraf | [END_REF]15], use can be made of the Poisson equality to replace the three-dimensional delta-function

δ(r -r ij ) = - 1 4π ∆ r i 1 |r -r ij | . (3)
Insertion of the Laplacian with respect to either r i or r j in the canonical average and integration by part yields (after symmetrization)

h ab (r) = - β ab 4π dΩ Na i=1 N b j=1 r ij -r |r ij -r| 3 • 1 2 (F j -F i ) (4)
where h ab = g ab -1.

Using the Gauss theorem for the electric field created by a uniformly charged sphere of radius r at the location

r ij dΩ r ij -rΩ |r ij -rΩ| 3 = r ij r 3 ij H(r ij -r), (5)
with H the Heaviside function, we get from eqn 4

h ab (r) = - β ab 4π Na i=1 N b j=1 1 2 (F j -F i ) • r ij r 3 ij H(r ij -r) (6)
This is a key formula of this work. Compared to the standard histogram procedure, it now involves the force acting on the particles in addition to their positions. It also implies a quite different numerical procedure: Here, for each configuration, every particles pair contributes to all distances r < r ij instead of just to r = r ij . Furthermore application of the formula requires a predefined grid but does not necessarily imply the limit of infinitely small grid separation, ∆r → 0. The only requirement is that the chosen grid, not necessarily regular, is adapted to cope with the variations of the pair distribution at all distances.

The new procedure is illustrated in Figure 1 for the radial distribution function, g(r), of a pure Lennard-Jones fluid composed of 864 particles at a reduced density ρ * = 0.8 and reduced temperature T * = 1.35, computed by molecular dynamics simulation.

We have displayed the RDF computed from ONE SINGLE equilibrated configuration using either histograms or eqn 6 and we compare those "instantaneous" curves to the converged result after 10000 time steps. In both approaches we used the same regular grid with ∆r = 0.005σ. It can be seen that the curve obtained by eqn 6 is already very smooth and quite close to the final converged result. The histogram curve does contain the converged one within its fluctuations but appears very noisy. This is further illustrated in Fig. 2 where we plot the variances

v(r) = 1 T t g t (r) 2 - 1 T t g t (r) 2 (7)
obtained after T = 1000 simulation steps; g t (r) is the "instantaneous" pair distribution function measured at step number t. It can be verified that the variance measured with the "force" approach is indeed much reduced with respect to the histogram approach and appears independent of the chosen grid size; the histogram method leads to a variance that is inversely proportional to ∆r.

0 1 2 3 4 r (σ) 0 1 2 3 4 g(r)
In Fig. 3 is displayed the oxygen-oxygen radial distribution obtained from 100 configurations of a DFT-MD simulation of 128 water molecules at ambient thermodynamic conditions with BLYP functional, after preliminary equilibration. It can be seen that, even for the relatively fine grid chosen and with a very limited number of steps, the "force RDF" is very smooth. It should also be noted that, even if the agreement with the converged RDF is already and overall very good after such a short trajectory, one observes slight discrepancies, in particular in the height of the first peak. At this stage, the force method is able to improve the variance of the RDF, but does not correct for the bias in the statistics. Looking for an improved method with minimum variance AND minimum bias is still another matter. Nonetheless, Fig. 3 is meant to show that

0 1 2 3 4 r (σ) 0 0,05 0,1 v(r)
Figure 2: Variance on the value of the radial distribution function depicted in figure 1 after 1000 simulation steps using either the force approach, eqn 6, or the standard histogram technique, with a grid spacing ∆r = 0.01σ (black and blue lines, respectively)

or ∆r = 0.005σ (cyan and red lines, respectively).

it is certainly in the field of ab-initio studies, where the generation of the trajectories themselves is computationally very expensive, that the force method described here to compute the RDF's reveals its full potential, given that the forces on the nuclei at each time step are readily available from the simulations.

In order to gain some physical insight into the above formula, one can write

H(r ij -r) = ∞ r dr δ(r -r ij ) (8)
Replacing into eq. 6 and inverting integral and canonical average, we get

ρ b h ab (r) = -β ∞ r dr F (r ) (9)
with the mean force density defined by We are back to a histogram procedure, but for F (r) instead of g ab (r) directly.

F (r) = ρ b ab 4πr 2 Na i=1 N b j=1 1 2 (F j -F i ) • r ij r ij δ(r -r ij ) . (10)
Denoting by F (r) the constrained (or conditional) mean force at a given distance, F (r) = F (r)/ρ b g ab (r), we get by differentiation of eqn 9, and division by g ab (r):

1 g ab (r) dg ab dr = F (r) = -β dw ab (r) dr , (11)
which is the definition of the potential of mean force (PMF), g ab (r) = exp(-βw ab (r)).

A fundamental difference of the present approach with respect to standard PMF calculations is the use of the force density F (r) instead of the mean force F (r). Eqn 9 is thus reminescent of the usual PMF formula but not equivalent to it in practice: Here, the integration of the force density F (r) yields g ab (r) directly rather than its logarithm and this quantity is computed in its integrality during the simulation, rather than step by step using constraints or restraints on the a -b distance as in usual PMF calculations.

Note also that eqs 6 and 9-10 are rigorously equivalent only in the limit of infinitely small grid size.

We note in passing that after eqn 4, a second integration by part can be performed, yielding

ρ b h ab (r) = β ab 4π dΩ Na i=1 N b j=1 1 2 ( Φ i + Φ j ) 1 |r -r ij | (12)
with Φ i = βF 2 i -∆ r i U .

Using now the Gauss theorem for the electrostatic potential created by a uniformly charged sphere of radius r at a location r ij

dΩ 1 |r ij -rΩ| = F (r, r ij ) = min( 1 r , 1 r ij ), (13)
eqn 12 can also be re-written as

ρ b h ab (r) = β ab 4π Na i=1 N b j=1 1 2 ( Φ i + Φ j ) F (r, r ij ) . (14)
As for eqn 6, this formula has to be evaluated at all points r of a pre-defined grid by summation over all particle pairs for each microscopic configuration. The key quantities to be evaluated for each particle are now the Φ i 's, that are scalar quantities involving the second derivatives of the potential felt by each particle (minus the mean square force times β).

In analogy to eqs 9-10, one can define also the continuous function Φ(r),

Φ(r) = ρ b ab 4πr 2 Na i=1 N b j=1 1 2 (Φ j + Φ i )δ(r -r ij ) , (15)
that can be evaluated through standard histograms for a given simulation. This function being computed, it is easy to show that the radial distribution can be extracted at the end through the one-dimensional integral relation

ρ b h ab (r) = β ∞ r dr r (1 -r /r)Φ(r ) (16)
This second formula appears much less familiar than the mean force formula of eqn 10.

We found its application, as well as that of the direct formula 14, to be very unstable and thus quite disappointing. We attribute this to the large fluctuation of the Φ i 's that appear as the difference of two large numbers. We do not have a clear and final understanding of this fact yet. For the time being, we limit ourselves to the force formulation, eqs 6 and 9, for which the reduced variance properties was shown to work quite fine.

Three-dimensional densities

The related problem that we address now is the computation of the three-dimensional density of a molecular solvent around a fixed molecule of arbitrary shape located at the origin using molecular simulations such as Monte-Carlo or Molecular Dynamics.

This density is defined by

ρ(r) = N i=1 δ(r -r i ) . ( 17
)
This quantity is usually computed by discretizing space into a 3D-grid and monitoring the filling of 3D-histograms of volume ∆V in the course of the simulation. This statistical process is indeed of infinite variance as ∆V → 0 since the instantaneous density in each bead oscillates betweeen 0 and O(1/∆V ). The exploration of the three-dimensional volume requires much more statistics than that required for radial functions and it is indeed known that three-dimensional densities are hard to get accurately that way. To bypass that problem, the Dirac function in eqn 17 may be again transformed according to the Poisson's identity of eqn 3. Substitution in the average above and integration by part with respect to the Boltzmann weight gives

∆ρ(r) = ρ(r) -ρ 0 = - β 4π N i=1 r -r i |r -r i | 3 • F i ( 18
)
F i = -∇ r i U representing the force on molecule i, U is the total interaction potential.

A second integration by part yields

∆ρ(r) = - β 4π N i=1 Φ i |r -r i | (19) with Φ i = βF 2 i -∆ r i U .
Those two formulas have a different content that eqn 17 since every particle now contributes instantaneously to every grid point in space, even those at long distances, and not only to the closest ones.

The straight application of those formulas is time consuming however, with a cost of order N × N g , where N g is the number of grid points, for each solvent configuration. This might not be a limitation in ab-initio molecular dynamics where the cost of a time step is already very high. In molecular mechanics however, the cost should be compared to N 2 or even N ln(N ) for advanced algorithms. A compromise has therefore to be found.

Inspired by the developments of the previous section, we can introduce a deltafunction in r into eqn ( 18)

∆ρ(r) = - β 4π N i=1 dr δ(r -r i ) r -r |r -r | 3 • F i (20)
and with the definition of the force density at point r

F(r ) = N i=1 δ(r -r i )F i (21) one gets ∆ρ(r) = - β 4π dr r -r |r -r | 3 • F(r ). (22)
In practice, F(r ) can be computed by histograms during a simulation and at the end Fourier-transformed to F(k) using a discrete three-dimensional fast Fourier transform

(FFT). ρ(r) is then obtained by inverse FFT of ∆ρ(k) = iβ k 2 k • F(k). (23)
Defining in the same way

Φ(r ) = N i=1 δ(r -r i )(βF 2 i -∆ r i U ) (24)
one has equivalently generated with CP2K code using the conventional BLYP functional. It is seen that the integrated force method described by eqs. 22-23 yields smooth densities that are basically independent of the chosen histogram bin size, whereas the straight histogram method yields noisy densities with fluctuations increasing with the fineness of the grid. Restricting the sampling, which is rather high here, to a smaller number of configurations yields analogous results with coarser grids.

∆ρ(r) = - β 4π dr Φ(r ) |r -r | (25) or ∆ρ(k) = -β Φ(k) k 2 . ( 26
) Noting that in fact Φ(r) = ∇ r • F(r), or in k-space, Φ(k) = -ik • F(k),
We note also that, as already found for the RDF's, the use of the second derivative approach using the function φ(r) in eqs. 24-26 instead of the above first-derivative approach with the forces yields numerical instabilities in the computation that we attribute to the very large statistical fluctuations of the molecular quantities Φ i 's; we still do not fully understand this behaviour.

Conclusions

We have derived in this paper two rather unusual formulas to compute the pair distribution functions in a bulk liquid mixture, namely eqs. 6 and 14. They require at each time step and for each pair distance realization to accumulate statistics over a whole pre-defined grid rather than just updating one histogram bin corresponding to that distance, as it is usually done. The first formula was shown to be related to the computation of the potential of mean force between two particles. It leads in practice to a computation on-the-fly that is quite different from the constraints/restraints method that is usually prescribed. The variance of the results is reduced by a large factor with respect to histogram techniques and was shown to be insensitive to the grid size. The second formula, that was associated to a zero-variance principle in the electronic density case [START_REF] Assaraf | [END_REF]15], seems impaired by the large fluctuations of the second-order quantities Φ i and did not prove useful at this stage. We think that the force approach is already very valuable in the context of ab-initio molecular dynamics simulations to improve the usually poor statistics due to short runs and small number of particles. It benefits moreover from the fact that in such simulations the forces on individual nuclei are already available.

The formalism has been extended to the computation of three-dimensional liquid densities around given molecular objects, quantities that are submitted to large statistical fluctuations when tackled by histogram techniques. Mere application of the Poisson identity for the three-dimensional Dirac function and successive integrations by parts, as proposed by Assaraf et al for electronic densities [START_REF] Assaraf | [END_REF]15], yield two formal expressions of the density, eqs. 18 and 19, which, to our knowledge, are unexploited in a statistical mechanics context. Straight application of the formulas requires the definition of a pre-determined 3D-grid and, for all grid points at each time step, an appropriate summation over all particles. The associated computation cost is largely affordable in an AIMD context. For classical simulations, in order not to impair the N 2 cost of the simulations, we have proposed an alternative, approximate scheme in which the mean forces (or mean second-order derivatives) are accumulated locally by histograms, and a final convolution is performed at the end with the help of 3D Fourier transforms. The first order scheme, using the forces, was shown to produce smooth densities with a much reduced variance compared to histograms. We hope to be able to clarify why the second-order scheme, aimed at giving a zero variance, turns out not to function so well in practice.

In any case, we think that the first order formalism described above will already prove very useful for e.g. the characterization of water densities around complex solutes such as protein, in order to estimate solvation free-energies, and to map local properties such as local entropies and enthalpies [16]. The introduction of the angular variables in the formalism in order to describe the molecules orientations in addition to their positions will be a necessary generalization in that context. Last but not least, extension to densities along general collective variables seems possible by using an auxiliary vector field as introduced in Ref. [7].

Figure 1 :

 1 Figure 1: Radial distribution function obtained for a SINGLE equilibrated configuration of a Lennard-Jones liquid composed of 864 particles using either the force approach, eqn 6, or the standard histogram technique, with a grid spacing ∆r = 0.005σ The dashed blue line indicates the converged result after 10000 simulation steps.

Figure 3 :

 3 Figure 3: Oxygen-oxygen radial distribution function averaged over 100 configurations extracted from a DFT-MD trajectory with 128 water molecules at ambient liquid conditions. The dashed blue line indicates the converged result obtained by averaging over 38600 configurations.

Figure 4 :

 4 Figure 4: Three-dimensional water density around a fixed water molecule from 36800 DFT-MD configurations of 128 molecules, represented in three orthogonal directions: z along the H 1 OH 2 bissector, y perpendicular to the molecular plane, x = z × y. The red and black lines are the results of the histogram and force methods, respectively. The grid spacing is ∆r = 0.06 Å for the top figure and ∆r = 0.04 Å for the bottom one.

This contribution is dedicated to Prof. Giovanni Ciccotti, in honour of his 70th birthday


‡ Rs2e

Supercapacitors based on an ionic liquid electrolyte and graphite or nanoporous carbon electrodes are simulated using molecular dynamics. We compare a simplified electrode model in which a constant, uniform charge is assigned to each carbon atom, with a realistic model in which a constant potential is applied between the electrodes (the carbon charges are allowed to fluctuate). We show that the simulations performed with the simplified model do not provide a correct description of the properties of the system. Firstly, the structure of the adsorbed electrolyte is partly modified. Secondly, dramatic differences are observed for the dynamics of the system during transient regimes. In particular, upon application of a constant applied potential difference, the increase in the temperature, due to Joule effect, associated with the creation of an electric current across the cell follows Ohm's law, while unphysically high temperatures are rapidly observed when constant charges are assigned to each carbon atom.

Cations

Supercapacitors are used for the electrochemical storage of energy. Unlike Li-ion batteries, there is no redox process involved in the bulk of an active material, and the charge is accumulated through the electrosorption of ions from an electrolyte on the electrode surface, which confers on supercapacitors very specific properties. [START_REF] Simon | Materials for Electrochemical Capacitors[END_REF] First, they show very high power densities and are therefore used in applications where there is a need for a large charge uptake or delivery in a limited amount of time (typically a few seconds). Second, only a slight volume change (< 3 %) occurs in the material during operation and the electrolyte is not subject to electrochemical degradation, and thus most devices show an outstanding cycle life (> 1,000,000 cycles). They have nevertheless an inherently moderate energy density, and there is an increased need for higher-energy devices. To this end a better understanding of the electrolyte/electrode interface is required, for which molecular simulations appear as a promising tool, [START_REF] Kondrat | A Superionic State in Nano-Porous Double-Layer Capacitors: Insights from Monte Carlo Simulations[END_REF][START_REF] Merlet | On the Molecular Origin of Supercapacitance in Nanoporous Carbon Electrodes[END_REF][START_REF] Wu | Voltage Dependent Charge Storage Modes and Capacity in Sub-nanometer Pores[END_REF][START_REF] Vatamanu | Molecular Insights into the Potential and Temperature Dependences of the Differential Capacitance of a Room-Temperature Ionic Liquid at Graphite Electrodes[END_REF][6] complementary to electrochemical [START_REF] Chmiola | Anomalous Increase in Carbon Capacitance at Pore Sizes Less Than 1 Nanometer[END_REF][START_REF] Chmiola | Desolvation of Ions in Subnanometer Pores and Its Effect on Capacitance and Double-Layer Theory[END_REF] or spectroscopic [START_REF] Wang | Real-Time NMR Studies of Electrochemical Double-Layer Capacitors[END_REF][START_REF] Richey | In Situ Molecular Level Measurements of Ion Dynamics in an Electrochemical Capacitor[END_REF] techniques.

At such an electrified interface, electrostatic interactions play a predominant role. The objective of this work is to study the importance of accounting realistically for electrode polarisation effects in the presence of an electrolyte. To this end, we perform a first series of simulations in which we impose that the potential on each electrode atom at each molecular dynamics step is constant and equal to a specified value 11 V i 0 , which is the electrode potential in a conventional electrochemistry experiment [START_REF] Reed | Electrochemical Interface between an Ionic Liquid and a Model Metallic Electrode[END_REF] (V i 0 = V + or V -for atoms in the positive and negative electrodes, respectively, with the potential difference equal to V + -V -). In practice, this is done by treating the partial charges carried by the electrode carbon atoms as additional degrees of freedom, which fluctuate during the simulation. Their values are determined from a self-consistent calculation ; this approach has now been used in several simulation studies [START_REF] Merlet | On the Molecular Origin of Supercapacitance in Nanoporous Carbon Electrodes[END_REF][START_REF] Vatamanu | Molecular Insights into the Potential and Temperature Dependences of the Differential Capacitance of a Room-Temperature Ionic Liquid at Graphite Electrodes[END_REF]6,[START_REF] Reed | Electrochemical Interface between an Ionic Liquid and a Model Metallic Electrode[END_REF][START_REF] Pounds | Ion Adsorption at a Metallic Electrode: an Ab Initio Based Simulation Study[END_REF][START_REF] Vatamanu | Molecular Simulations of the Electric Double Layer Structure, Differential Capacitance, and Charging Kinetics for N-Methyl-Npropylpyrrolidinium Bis(fluorosulfonyl)imide at Graphite Electrodes[END_REF][START_REF] Merlet | Imidazolium Ionic Liquid Interfaces with Vapor and Graphite: Interfacial Tension and Capacitance from Coarse-Grained Molecular Simulations[END_REF][START_REF] Xing | On the Atomistic Nature of Capacitance Enhancement Generated by Ionic Liquid Electrolyte Confined in Subnanometer Pores[END_REF] . Were this electrostatic response found negligible, a significant gain of simulation time would be obtained. In order to assess the importance of these effects, simulations performed at constant applied potential are compared here with a second series of simpler, computationally less expensive simulations where fixed partial charges are assigned on each atom, as reported in a majority of simulation studies. [START_REF] Lanning | Screening at a Charged Surface by a Molten Salt[END_REF][START_REF] Pinilla | Polarization Relaxation in an Ionic Liquid Confined between Electrified Walls[END_REF][START_REF] Lynden-Bell | Simulation of Ionic Liquids, Solutions and Surfaces[END_REF][START_REF] Fedorov | Ionic Liquid near a Charged Wall: Structure and Capacitance of Electrical Double Layer[END_REF][START_REF] Feng | Microstructure and Capacitance of the Electrical Double Layers at the Interface of Ionic Liquids and Planar Electrodes[END_REF][START_REF] Kislenko | Molecular Dynamics Simulation of the Electrochemical Interface between a Graphite Surface and the Ionic Liquid [BMIM][PF 6[END_REF][START_REF] Shim | Nanoporous Carbon Supercapacitors in an Ionic Liquid: A Computer Simulation Study[END_REF][START_REF] Si | Effect of Specific Adsorption on the Differential Capacitance of Imidazolium-Based Ionic Liquid Electrolytes[END_REF][START_REF] Li | Molecular Dynamics Simulation Study of the Capacitive Performance of a Binary Mixture of Ionic Liquids near an Onion-like Carbon Electrode[END_REF][START_REF] Paek | A Computational Study of the Interfacial Structure and Capacitance of Graphene in[END_REF] In both cases, the same coarsegrained model is used to describe the electrolyte. [START_REF] Merlet | Imidazolium Ionic Liquid Interfaces with Vapor and Graphite: Interfacial Tension and Capacitance from Coarse-Grained Molecular Simulations[END_REF][START_REF] Merlet | New Coarse-Grained Models of Imidazolium Ionic Liquids for Bulk and Interfacial Molecular Simulations[END_REF] We show that the simplified simulation scheme does not reproduce well a realistic electrode, both from the structural and the dynamic point of view.

In a first step we have focused on a supercapacitor simulation cell based on a 1-butyl-3methylimidazolium hexafluorophosphate (BMIM-PF 6 ) ionic liquid electrolyte and graphite electrodes. This system was chosen because it is well known [START_REF] Merlet | Imidazolium Ionic Liquid Interfaces with Vapor and Graphite: Interfacial Tension and Capacitance from Coarse-Grained Molecular Simulations[END_REF][START_REF] Merlet | New Coarse-Grained Models of Imidazolium Ionic Liquids for Bulk and Interfacial Molecular Simulations[END_REF] and because the planar surface of the electrodes facilitates the analysis of the adsorbed electrolyte structure. In order to perform a fair comparison, the constant charges were assigned based on the average total charge measured during constant potential simulations for various applied potential differences (see the Supplementary Information). This total charge was equally distributed among the carbon atoms, and the selfconsistency of the procedure was checked by calculating the potential difference arising from the so-built constant charge simulations: The initial values were well reproduced.

The ionic density profiles along the normal to the graphite surfaces are shown on Figure 1. On each electrode, only the profiles corresponding to the oppositely charged ions are shown for clarity.

The region corresponding to z values between 3.0 and 9.3 nm is not shown because it corresponds to the bulk for which ionic densities remain flat under all conditions. Although the constant poten- tial simulations (top panel) profiles look very similar to the constant charge ones (bottom panel) at first glance, some differences are observed for both electrodes. Indeed, for the positive electrode (left panels), the second peak of the PF - 6 ions density profiles, located at distances between 1.5 and 2.0 nm, shows a different variation depending on the simulation scheme. In the first case, the intensity remains constant whatever the applied potential difference, but the peak is shifted toward larger distances upon increase of the electrode polarization, switching from 1.58 nm for 0 V to 1.73 nm for 9 V. When constant charges simulations are employed, not only is this shift more pronounced, but the peak also becomes much flatter. In the case of the BMIM + cations (right panels), the changes concern the first intense peak, which also becomes less intense and flatter in the case of the constant charge simulations. Finally, we observe that the differences in the structural properties become more visible when the electrodes are highly charged. A comparison of the full ionic density profiles for an applied potential difference of 9.0 V or a constant charge density of 1 e nm -2 (Figure 2) clearly shows that most of the oscillatory profile close to the interfaces is lost in the second case. The difference between these density profiles is due to the fact that in the second case, the constant charge on the electrode constrains the liquid to adopt a given structure. When a constant potential is applied, there is an interplay between the charge on the electrode and the local structure of the electrolyte layer facing it, which is therefore allowed to change more smoothly. It is worth mentioning that the variation of the charge with respect to the potential drop at the electrode/electrolyte interface obtained with the constant charge simulations agrees very well with that obtained with the constant potential simulation ones, resulting in similar values for the differential capacitance (data provided in the Supplementary Information). This result may however be partly fortuitous, and should not be taken as a general rule.

The performance of supercapacitors in terms of power (i.e how fast the charge can be delivered) can be assessed by modelling the polarization relaxation during transient regimes. This is done by performing simulations where the applied potential difference or the charge is suddenly changed and by characterizing the processes which occur before the new steady-state is reached [START_REF] Lanning | Screening at a Charged Surface by a Molten Salt[END_REF][START_REF] Pinilla | Polarization Relaxation in an Ionic Liquid Confined between Electrified Walls[END_REF] (note that the objective of such simulations is not to mimic the cyclic voltammetry setups, in which the potential is applied with a slow scan rate, but rather to study the relaxation phenomena which occur in a supercapacitor under operation). The difference between constant potential and constant charge simulations is then much more spectacular. Although it is already visible in the case of graphite electrode systems, we have picked a more illustrative example by simulating a supercapacitor with nanoporous carbide-derived carbon electrodes. [START_REF] Merlet | On the Molecular Origin of Supercapacitance in Nanoporous Carbon Electrodes[END_REF] We have performed 3 different simulations, starting from the same initial configuration. The latter was extracted from a trajectory obtained with uncharged electrodes, and presents the particularity to have some empty pores.

The two first simulations were performed with constant applied potential differences of 1.0 V and 5.0 V respectively, while the third one was performed with a constant charge of 0.01 e (respectively -0.01 e) assigned to each carbon atom of the left (right) electrode: This value corresponds to the average of the fluctuating charges carried by the carbon atoms for a 1 V applied potential simulation. [START_REF] Merlet | On the Molecular Origin of Supercapacitance in Nanoporous Carbon Electrodes[END_REF] As shown on Figure 3, after only 5 ps, the empty pore is already filled for the constant charge simulation. This is not the case for any of the constant potential ones, for which 100 ps is not sufficient to obtain a full charge (only the 1.0 V situation is shown on Figure 3 for planar electrodes. These authors concluded that constant potential simulations provided much more reasonable double layer relaxation times. [START_REF] Vatamanu | Molecular Simulations of the Electric Double Layer Structure, Differential Capacitance, and Charging Kinetics for N-Methyl-Npropylpyrrolidinium Bis(fluorosulfonyl)imide at Graphite Electrodes[END_REF] Here, even more striking is the increase in the temperature due to the Joule effect, as shown in Figure 4. In the constant potential simulations, the calculated heating rate is of 0.08 K ps -1 (respectively 2.1 K ps -1 ) for a constant applied potential of 1.0 V (respectively 5.0 V), in very good agreement with a basic estimation based on Ohm's law (see the Supplementary Information for details) of 0.10 K ps -1 (respectively 2.5 K ps -1 ). In contrast, in the constant charge simulations the temperature reaches 25,000 K in 1.0 ps, reflecting an extreme energy dissipation and the driving of the system far from equilibrium, with a strong macroscopic electric field across the cell. The relaxation processes observed during the latter simulations are therefore unrealistic and this procedure should be avoided for the study of transient regimes occurring in supercapacitors at the molecular scale. Nevertheless, the constant potential simulations with nanoporous carbon electrodes are very expensive from the computational point of view, so that constant charge simulations remain very useful in the first stages of equilibration of the simulation cells. In the case of constant applied potential simulations, the charge on the electrode atoms varies in response to the ionic liquid; as a consequence the charging of the initially empty pores occurs gradually while they become wetted, resulting in much slower relaxation times. In conclusion, we have showed that the use of constant charge simulation setups for modelling supercapacitors based on ionic liquid electrolytes and carbon electrodes alter both the structure of the adsorbed fluid at the interface and the timescales over which relaxation phenomena occur. The situation is much more critical for the latter, which are modified by orders of magnitude, but the changes in the structure are also visible and they increase with the charge of the electrode, which may lead to modifications in the thermodynamic properties of the systems. This situation can be paralleled with the role played by ion polarization effects in bulk ionic liquids and molten salts:

T (K) V + -V -= 1.0 V 5.0 V 0 1 2 3 4 5 t (ps)
In the latter case, it is now well known that the simulations performed in the framework of rigid ion models provide a correct description of the structure of singly-charged ions based liquids, but, as soon as multivalent ions are added, such models are not able to yield the complex structural patterns which are observed and it is necessary to resort to a polarizable ion model. [START_REF] Borodin | Polarizable Force Field Development and Molecular Dynamics Simulations of Ionic Liquids[END_REF][START_REF] Bedrov | Influence of Polarization on Structural, Thermodynamic, and Dynamic Properties of Ionic Liquids Obtained from Molecular Dynamics Simulations[END_REF][START_REF] Schroder | Comparing Reduced Partial Charge Models with Polarizable Simulations of Ionic Liquids[END_REF][START_REF] Salanne | Polarization Effects in Ionic Solids and Melts[END_REF] Note that here the use of more complex force fields for the ionic liquids, which include polarization effects either explicitly or effectively [START_REF] Chaban | A New Force Field Model for the Simulation of Transport Properties of Imidazolium-Based Ionic Liquids[END_REF][START_REF] Chaban | A New Force Field Model of 1-Butyl-3-Methylimidazolium Tetrafluoroborate Ionic Liquid and Acetonitrile Mixtures[END_REF] (and even the electron density breathing effects [START_REF] Madden | Covalent' Effects in 'Ionic' Systems[END_REF] ), would also enhance the description of the system. Finally the dynamics is highly affected, and transport properties are usually underestimated by one order of magnitude with rigid ion models. [START_REF] Bedrov | Influence of Polarization on Structural, Thermodynamic, and Dynamic Properties of Ionic Liquids Obtained from Molecular Dynamics Simulations[END_REF][START_REF] Schroder | Comparing Reduced Partial Charge Models with Polarizable Simulations of Ionic Liquids[END_REF][START_REF] Salanne | Polarization Effects in Ionic Solids and Melts[END_REF] Here the effect is opposite, since constant charge simulations greatly enhance the dynamics of the system.

It is therefore important to use constant applied potential conditions when modelling electrodes, despite the associated increase in computational cost. This recommendation is likely to hold for the modelling of systems other than supercapacitors which involve a metal / ionic liquid interface; 35 for example the role of polarization on the stabilization of metallic nanoparticles in ionic liquids has recently been demonstrated. [START_REF] Pensado | Solvation and Stabilization of Metallic Nanoparticles in Ionic Liquids[END_REF] It is worth noting that schemes alternative to the selfconsistent calculation of the charges have been proposed, [START_REF] Petersen | A Computationally Efficient Treatment of Polarizable Electrochemical Cells held at a Constant Potential[END_REF] which allow to decrease significantly this computational time, but these methods have yet to be adapted to the case of porous electrodes.

Computer simulation details

Following our previous works, the coarse-grained model of Roy and Maroncelli, which yields the correct structural, thermodynamic and dynamic properties of BMIM-PF 6 , was used for the ionic liquid. [START_REF] Roy | An Improved Four-Site Ionic Liquid Model[END_REF] In all the simulations, the Coulomb interactions are calculated through a two-dimensional Ewald summation [START_REF] Reed | Electrochemical Interface between an Ionic Liquid and a Model Metallic Electrode[END_REF][START_REF] Gingrich | On the Ewald Summation of Gaussian Charges for the Simulation of Metallic Surfaces[END_REF] and a time step of 2 fs, which ensures a good conservation of the total Hamiltonian, is used to integrate the equations of motion. For the ionic liquid-graphite simulations, the system consists of 320 ion pairs sandwiched between two groups of three (respectively one)

graphite layers for the constant potential (respectively charge) simulations. In the latter case, one plane only was used in order to match with previous studies. The lengths in x and y directions are chosen to match with the lattice parameters of graphite and are equal to 3.22 and 3.44 nm, respectively. The length in the z direction is adjusted to recover the bulk density of the ionic liquid at 400 K leading to L z = 12.32 nm. For each point, the system is equilibrated during 200 ps before gathering data during 1 ns in the NVT ensemble. For the ionic liquid-porous carbon simulations, the system consists of 600 ion pairs sandwiched between two porous electrodes [START_REF] Palmer | Modeling the Structural Evolution of Carbide-Derived Carbons Using Quenched Molecular Dynamics[END_REF] containing 3649 atoms each, as explained in reference. [START_REF] Merlet | On the Molecular Origin of Supercapacitance in Nanoporous Carbon Electrodes[END_REF] The length in x and y directions are both equal to 4.37 nm, while the length in the z direction is of 18.64 nm. The systems are simulated during 100 ps in the NVE ensemble, starting from an initial configuration as explained in the main text. Note that the number of simulated ions was enough to keep a constant density in the bulk region during all the simulations.

Figure 1 :

 1 Figure 1: Ionic density profiles for the constant applied potential (top) and constant charges (bottom) simulations. For each electrode, only the oppositely charged ions profiles are shown for clarity. The arrows highlight the main differences between the two sets of data, which are discussed in the text.

2 Figure 2 :

 22 Figure 2: Ionic density profiles for the 9.0 V constant applied potential (V + = -V -= 4.5 V) and the 1 e nm -2 constant charge simulations. Left: positive electrode, right: negative electrode. Similar plots for other electrode potentials are provided in the Supplementary Information.

-Figure 3 :

 3 Figure 3: Evolution of the structure inside the positive electrode upon application of a constant charge (top) or a constant applied potential (bottom) simulation. Only part of the simulation cell is shown here. The initially empty pore region is highlighted with a red circle. Green spheres: anions, red spheres: cations.

Figure 4 :

 4 Figure 4: Temperature drift when a constant potential difference is suddenly applied (insert: when a constant charge is assigned to each carbon atom).

Acknowledgements

We acknowledge the support of the French Agence Nationale de la Recherche (ANR) under grant ANR-2010-BLAN-0933-02 ('Modeling the Ion Adsorption in Carbon Micropores'). The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement n. 102539. We are grateful for the computing resources on JADE (CINES, French National HPC) obtained through the project x2012096728. This work made use of the facilities of HECToR, the UK's national high-performance computing service, which is provided by UoE HPCx Ltd at the University of Edinburgh, Cray Inc and NAG Ltd, and funded by the Office of Science and Technology through EPSRC's High End Computing Programme.

Supplementary information

Additional methodology details and supplementary figures are all given in the Supplementary Information. Two movies showing the trajectories used for Figure 3 are also provided. Finally, in a third movie we show the full trajectory (100 ps) for the 1 V constant applied potential simulation.


Benjamin Rotenberg

Richard Taïeb

Valérie Véniard

Alfred Maquet

H2+ in intense laser field pulses: ionization versus dissociation within moving nucleus simulations

Keywords: numbers: 36, 40, -c, 32, 80, Rm

come

Introduction

The interaction of intense laser pulses with either atoms, molecules or clusters has motivated a large number of studies describing interesting non-pertubative effects such as high-order harmonic or X-rays generation [START_REF] Huillier | Atoms in Strong Fields[END_REF], above threshold ionisation -dissociation- [START_REF] Muller | Atoms in Strong Fields[END_REF][START_REF]Molecules in Laser Fields[END_REF][START_REF] Giusti-Suzor | [END_REF] or very energetic ions emerging from Coulomb explosion [5] . The major difficulty with molecules is to deal with both the electronic and the nuclear degrees of freedom. Regarding the latter, one has often neglected it, through the use of the Born-Oppenheimer approximation, keeping the nuclei fixed, [6,7], or considered only the lowest, or essential, electronic states, [START_REF] Giusti-Suzor | [END_REF], ignoring the role of ionisation. More refined approaches are based on a full quantum description of the system in 1-D, [8,9] or 3-D, [10,11]. Here we propose to use a semi-classical approach, somewhat similar to the one used routinely in collision theory, [12] i.e. we describe the electron quantum mechanically while the motion of the two nuclei is described classically. In this letter, we restrict ourself to a 1-D configuration, using soft-core potentials, since such a model describes reasonably well both the atomic electron, [START_REF] Eberly | Atoms in Strong Fields[END_REF], and the molecular dynamics, [START_REF] Pegarkov | [END_REF], in linearly polarized laser fields. We note that this approach can be straightforwardly extended to 3-D. This method gives an insight on what parameters are essential in the evolution of the molecular ion in the presence of a laser pulse and it sheds light on the competition between the different processes like ionisation, dissociation or excitation of the vibrational levels. We first present the model, then we compare the results obtained with either frozen or moving nuclei, exhibiting, in both cases, charge-resonanceenhanced-ionisation (CREI) already seen in previous works, [6,7,8,9,10]. We will use atomic units throughout this letter.

Model

The electronic state of the H + 2 ion is described by its wavefunction, the time evolution of which is derived by solving the following 1-D time dependent Schrödinger Equation (TDSE), in the velocity gauge:

i ∂ψ(x, t) ∂t = [- 1 2
∂ 2 ∂x 2 + V (x, t) + i c A(t). ∂ ∂x ]ψ(x, t) (1)
V (x, t) = -1

a 2 + (x -X 1 (t)) 2 - 1 a 2 + (x -X 2 (t)) 2 ,
where X 1 (t), X 2 (t) are the positions of the two nuclei. The electric field E(t) was chosen as E 0 f (t)sin(ωt), with f (t) representing the field envelope. It is related to the vector potential A(t) through the relation

E(t) = -1 c . ∂A(t)
∂t . The smoothing parameter a was set equal to √ 2 in order to get -0.5 a.u. as an asymptotic energy for the isolated atom. The motion of the protons is treated classically via the equations:

M d 2 X i dt 2 = F j→i (t) + F e -→i (t) + F laser (t), (2)
with F j→i (t), F e -→i (t), F laser (t) represent the repulsive internuclear, attractive electronic and laser forces, respectively. Their expressions are given by:

F j→i (t) = X i (t) -X j (t) (q + R 2 (t)) 3/2 (3)
F e -→i (t) = - +L -L (X i (t) -x) | ψ(x, t) | 2 (a 2 + (X i (t) -x) 2 ) 3/2 dx (4)
F laser (t) = E(t), (5)
where

R(t) = | X 1 (t) -X 2 (t)
| is the internuclear distance, and M ≈ 1836 a.u. the mass of the proton. The smoothing parameter q in the internuclear term was chosen equal to 0.03, in accordance with previous works, [8,9]. The time evolution of the different particles was obtained by the Crank-Nicholson method for Eq.( 1), [START_REF] Press | Numerical Recipes[END_REF] and the Verlet method for Eq.( 2), [START_REF] Verlet | [END_REF]. In order to avoid reflexions of the electronic density from the boundaries of the box we have used an absorber, [17]. Therefore, the norm of the wavefunction can decrease in time, providing a measure of the ionisation. We can also compute the positions of the two nuclei as well as the electronic probability density P i around each proton (in a box of 5 a.u.). The variations in time of the latter quantity can indicate dissociation or vibrational excitation of the molecular ion.

Results

We have first computed the energies of the two first molecular states, 1σ g and 1σ u , of the ion as determined by diagonalizing the laser-free hamiltonian for fixed internuclear distance, R. With the help of this simple model we have obtained the two first molecular states i.e. a pure dissociative state 1σ u , while 1σ g has a minimum located at R = 3.5 a.u. with a dissociation energy of 0.07 a.u.. These values, which are typical of 1-D simulation, are in reasonable agreement with results for full 3-D simulation (R e ≈ 2 a.u.).

Fixed nuclei ionisation rates

We have first determined the ionisation rates for a Nd:YAG laser (ω = 0.043 a.u., with a cycle duration of T L ≈ 3.5 f s). To this end we have used a trapezoidal laser pulse with two cycles linear turn-on and turn-off and 6 cycles flat part. The ionisation rates were deduced from the decrease of the norm of the wavefunction absorbed at the boundaries. The variations of the ionisation rate with the internuclear distance R for two intensities, I = 5 × 10 13 and 10 14 W/cm 2 are shown in Fig. 1. There exists an intensity dependent critical value for R c around 5 a.u. where ionisation starts to be enhanced. The value of R c increases when the intensity decreases. For larger values of R one observes local minima. For distances longer than 20 a.u., we recover the ionisation rate for an isolated 1-D hydrogen atom. This phenomenon can be explained in terms of Charge-Resonance states, using the quasi-static approximation and the associated tunelling rates, [18,19]. Briefly, in the presence of the low frequency electric field the HOMO 1σ g and the LUMO 1σ u states are transformed in localised atomic states with an energy separation of the order of E 0 R. Each state has its own ionisation rate, that is strongly dependent on its energy in the "laser-modified" potential barrier. Therefore, for particuliar values of R the electron can be localised in one of these states and can ionise more rapidly.

Moving nuclei simulation

As shown above, ionisation yields are very small for femtosecond laser pulses if the internuclear distance is smaller than R c ≈ 5 a.u.. However, as ionisation takes place, as the electronic density starts to expand, the bond is softened and, therefore, R increases. Then, if R c is reached while the laser is on, enhanced ionisation can take place, thus leading to the Coulomb explosion of the molecule. This is the case when, starting from the equilibrium internuclear distance R = 3.5 a.u., we take the same pulse shape as above with an intensity of I = 10 14 W/cm 2 . One can see, in Fig. 2, that the critical distance is reached for t ≈ 20 f s. Ionisation is the predominant phenomenon in this case with a final ionisation yield ≥ 99.5%. It should be noted that, for fixed nuclei, ionisation would only amount to less than 0.1% If we now consider more realistic gaussian pulse shapes, f (t) = e -t 2 τ 2 , we can easily identify three regimes -excitation of the molecular ion, dissociation and ionisationcorresponding to -very short, short and long -pulse durations as measured by the Full Width Half Maximum (FWHM = 2τ √ ln2) illustrated in Figs. 3. As said before, the value of R c increases when the intensity decreases. Therefore, the final state reached by the system will mainly depend on what time, in the pulse, this critical distance is reached. If we are in the long pulse regime, Fig. 3a,e.g. if R c is reached for t ≤ τ , we are in a situation close to the one of Fig. 2 followed by Coulomb explosion of the two bare nuclei, is dominating the dynamics. In the very short pulse regime, the ion is not stretched enough while the laser is on and there is almost no ionisation. One can see, in Fig. 3c, the two nuclei oscillating as H + 2 is vibrationaly excited. The period of oscillation (≈ 25 f s) is in agreement with the one we would expect from the potential curve of the 1σ g state. More interesting is the intermediate case (short pulse regime, Fig. 3b) corresponding to the situation when R c is reached for t ≈ τ , i.e. when the electric field intensity is already decreasing. Then, the system can evolve in two different channels leading either to excitation, as previoulsy described, or to dissociation. The latter corresponds to the case when there is separation of one ion, the electron staying close to the other nucleus. It turns out that the choice between either final state depends critically on the phase of the laser field, see Paulus et al, [20] for a similar situation in atoms. Interestingly, our simulation indicates that, when there is dissociation, the electron is in a superposition of the two states (1σ u and 1σ g ) and its probability of presence oscillates between the two nuclei. Then the nuclei move apart from each other, the internuclear potential barrier increases and the oscillating electron will end eventually being trapped in one of the potential holes. Again, the choice of which nuclei, right or left, depends on the laser phase.

Conclusions

In the present letter we have presented an analysis based on a 1-D model of the response of a simple molecular ion to an intense infrared laser pulse. Our findings show that it is necessary to include the motion of the nuclei in order to properly address the question of the competition between the different processes: excitation, ionisation (leading to Coulomb explosion) and dissociation. We have shown also that which final state is reached depends critically on the duration of the laser pulse and also on the absolute phase of the laser. It would be most interesting to extend this model to a 3-D treatement in order to analyse the effects of the alignment of the molecule on its behaviour.

Figure 1 .

 1 Figure 1. Ionisation rate as a function of the interatomic distance R for a YAG:Nd laser (ω = 0.043 a.u.). Solid line: I = 10 14 W/cm 2 ; Dashed line: I = 5 × 10 13 W/cm 2 . The critical internuclear distance R c (see text) is indicated by an arrow.

2 Figure 2 .

 22 Figure2. Ionisation yield (solid line), nuclei positions X i (dashed lines) and electronic probabilities P i around the two nuclei (dotted lines) for trapezoidal laser pulse (thick line) with I = 10 14 W/cm 2 and ω = 0.043a.u.

2 Figure 3 .

 23 Figure 3. Ionisation yield (solid line), nuclei positions X i (dashed lines) and electronic probabilities P i around the two nuclei (dotted lines) for a gaussian laser pulse (thick line) different durations. Here I = 10 14 W/cm 2 and ω = 0.043a.u.. (a) : F W HM = 50 f s, (b): F W HM = 14.5 f s, (c): F W HM = 7 f s.

in intense laser field pulses: Ionisation vs. Dissociation within moving nuclei simulations.5

Acknowledgements

The Laboratoire de Chimie Physique-Matière et Rayonnement is a Unité Mixte de Recherche, Associée au CNRS, UMR 7614, and is Laboratoire de Recherche Correspondant du CEA, LRC # DSM-98-16. Parts of the computations have been performed at the Centre de Calcul pour la Recherche (CCR, Jussieu, Paris) and at the Institut du Développement et des Ressources en Informatique Scientifique (IDRIS).