Hind el Hadri Alice Pradel

Cloé Desmet

Jessica Ponti

Stéphanie Reynaud

Bruno Grassl

Julien Gigault
email: julien.gigault@univ-rennes1.fr

Deposition of environmentally relevant nanoplastic models in sand during transport experiments

Keywords: Conceptualization, Methodology, Investigation, Formal analysis, Writingoriginal draft, review & editing, Visualization Characterization of Particle Hydrophobicity Transmission Electron Microscopy Analysis Conceptualization, Formal analysis, Writing -review Nanoplastic, Size, Shape, Transport, Soil, Column

Nanoplastics (NPTs) are defined as colloids that originated from the unintentional degradation of plastic debris. To understand the possible risks caused by NPTs, it is crucial to determine how they are transported and where they may finally accumulate. Unfortunately, although most sources of plastic are land-based, risk assessments concerning NPTs in the terrestrial environmental system (soils, aquifers, freshwater sediments, etc.) have been largely lacking compared to studies concerning NPTs in the marine system. Furthermore, an important limitation of environmental fate studies is that the NPT models used are questionable in terms of their environmental representativeness. This study describes the fate of different NPT models in a porous media under unfavorable (repulsive) conditions, according to their physical and chemical properties: average hydrodynamic diameters (200 to 460 nm), composition (polystyrene with additives or primary polystyrene) and shape (spherical or polymorphic). NPTs that more closely mimic environmental NPTs present an inhomogeneous shape (i.e., deviating from a sphere) and are more deposited in a sand column by an order of magnitude. This deposition was attributed in part to physical retention, as confirmed by the straining that occurred for the larger size fractions. Additionally, different Derjaguin-Landau-Verwey-Overbeek (DLVO) models -the extended DLVO (XDLVO) and a DLVO modified by surface element integration (SEI) methodsuggest that the environmentally relevant NPT models may alter its orientation to diminish repulsion from the sand surface and may find enough kinetic energy to deposit in the primary energetic minimum. These results point to the importance of choosing environmentally relevant NPT models.

Introduction:

It has been established that every environmental system contains plastic debris [START_REF] Bank | The Plastic Cycle: A Novel and Holistic Paradigm for the Anthropocene[END_REF]. The marine environment is the first system in which plastic pollution has been extensively studied, probably due to the highly visible accumulation that started in the 1970s [START_REF] Carpenter | Plastics on the Sargasso Sea Surface[END_REF][START_REF] Cozar | Plastic debris in the open ocean[END_REF][START_REF] Law | Plastic Accumulation in the North Atlantic Subtropical Gyre[END_REF][START_REF] Schwarz | Sources, transport, and accumulation of different types of plastic litter in aquatic environments: A review study[END_REF]. Since then, the extent of this contamination has been more closely described. Plastic fragments or additives are found in remote areas as they can travel long distances [START_REF] Jamieson | Bioaccumulation of persistent organic pollutants in the deepest ocean fauna[END_REF][START_REF] Peng | The ocean's ultimate trashcan: Hadal trenches as major depositories for plastic pollution[END_REF]. They have reached uninhabited forests after undergoing atmospheric transport [START_REF] Allen | Atmospheric transport and deposition of microplastics in a remote mountain catchment[END_REF] and have been incorporated into Artic sea ice after being transported by oceanic currents [START_REF] Obbard | Global warming releases microplastic legacy frozen in Arctic Sea ice[END_REF]. It is generally agreed upon that the terrestrial system is the major source of plastic waste [START_REF] Jambeck | Plastic waste inputs from land into the ocean[END_REF]. However, plastics in terrestrial and freshwater environments have received little attention [START_REF] Horton | Microplastics in freshwater and terrestrial environments: Evaluating the current understanding to identify the knowledge gaps and future research priorities[END_REF][START_REF] Rillig | Microplastic in Terrestrial Ecosystems and the Soil?[END_REF][START_REF] Wagner | Microplastics in freshwater ecosystems: what we know and what we need to know[END_REF]. Plastic accumulates in the terrestrial system due to common sources, such as faulty waste management systems, spreading of sewage sludge, use of plastic mulches and other plastic products in agriculture, and tire wear particles [START_REF] Carr | Transport and fate of microplastic particles in wastewater treatment plants[END_REF][START_REF] He | Influence of Nano-and Microplastic Particles on the Transport and Deposition Behaviors of Bacteria in Quartz Sand[END_REF][START_REF] Liu | Aging Significantly Affects Mobility and Contaminant-mobilizing Ability of Nanoplastics in Saturated Loamy Sand[END_REF][START_REF] Scheurer | Microplastics in Swiss Floodplain Soils[END_REF]. This rapidly accumulating plastic in soils and sediments can be degraded by photo-oxidation, thermo-oxidation, abrasion and biodegradation, which causes the inevitable release of small plastic particles [START_REF] Huerta Lwanga | Field evidence for transfer of plastic debris along a terrestrial food chain[END_REF][START_REF] Ng | An overview of microplastic and nanoplastic pollution in agroecosystems[END_REF]. In the terrestrial system (soils, aquifers, sediments, etc.), studies have generally focused on microplastics (1 μm to 5 mm) (Dong et al., 2019a(Dong et al., , 2019b;;[START_REF] He | Influence of Nano-and Microplastic Particles on the Transport and Deposition Behaviors of Bacteria in Quartz Sand[END_REF][START_REF] Liu | Aging Significantly Affects Mobility and Contaminant-mobilizing Ability of Nanoplastics in Saturated Loamy Sand[END_REF][START_REF] Tufenkji | Breakdown of Colloid Filtration Theory: Role of the Secondary Energy Minimum and Surface Charge Heterogeneities[END_REF], although it is clear that plastic fragmentation does not stop at the micrometer size [START_REF] Gigault | Marine plastic litter: the unanalyzed nanofraction[END_REF][START_REF] Lambert | Characterisation of nanoplastics during the degradation of polystyrene[END_REF].

The fragmentation of macro-and microplastics to nanoplastics (<1 μm) adds complexity to the global plastic waste issue. Even if no clear definition of nanoplastics (NPTs) is proposed by policymakers, a scientific consensus has been proposed to define NPTs as colloids that originate from the unintentional degradation of plastic [START_REF] Gigault | Current opinion: What is a nanoplastic?[END_REF][START_REF] Hartmann | Are We Speaking the Same Language? Recommendations for a Definition and Categorization Framework for Plastic Debris[END_REF]. Compared to microplastics, NPTs acquire colloidal properties and become unaffected by gravity. At this scale, the fate of NPTs will be controlled by surface properties, shape and diffusion rather than by the bulk properties [START_REF] Hiemenz | Principles of colloid and surface chemistry[END_REF]. The fact that the majority of plastic debris cannot be located by mass budgets at the ocean surface suggests that part of the missing fraction of plastic debris could be composed of fractions smaller than 300 µm in size (Albert A Koelmans et al., 2017;[START_REF] Van Sebille | A global inventory of small floating plastic debris[END_REF]. Since the nanoscale characteristics of plastic debris may pose potential risks, it is therefore crucial to determine where this plastic originates from and where it accumulates [START_REF] Galloway | Micro-and Nano-plastics and Human Health[END_REF]Albert A. Koelmans et al., 2017;[START_REF] Lehner | Emergence of Nanoplastic in the Environment and Possible Impact on Human Health[END_REF].

Due to limited data on the occurrence of NPTs, lab-scale models can aid in determining the transfer of NPTs from one environmental system to another. For the terrestrial system, a saturated sand column is generally used as a proxy for the study of colloid transfer in soils, sediments and aquifers [START_REF] Geitner | Measuring Nanoparticle Attachment Efficiency in Complex Systems[END_REF][START_REF] Lecoanet | Laboratory Assessment of the Mobility of Nanomaterials in Porous Media[END_REF][START_REF] Petosa | Aggregation and Deposition of Engineered Nanomaterials in Aquatic Environments: Role of Physicochemical Interactions[END_REF][START_REF] Redman | Bacterial Adhesion and Transport in Porous Media: Role of the Secondary Energy Minimum[END_REF][START_REF] Syngouna | Cotransport of clay colloids and viruses in water saturated porous media[END_REF][START_REF] Tufenkji | Correlation Equation for Predicting Single-Collector Efficiency in Physicochemical Filtration in Saturated Porous Media[END_REF]. Quevedo and Tufenkji compared the transport of spherical and size standardized polystyrene nanoparticles in quartz sand and loamy soil in the presence of monovalent or divalent salts [START_REF] Quevedo | Mobility of Functionalized Quantum Dots and a Model Polystyrene Nanoparticle in Saturated Quartz Sand and Loamy Sand[END_REF]. Recently, Hu et al. used quartz sand columns to investigate the cotransport of naphthalene with polystyrene nanoparticles [START_REF] Hu | Cotransport of naphthalene with polystyrene nanoplastics (PSNP) in saturated porous media: Effects of PSNP/naphthalene ratio and ionic strength[END_REF]. Another recent work demonstrated the influence of the soil type on the same nanoparticles [START_REF] Wu | Effect of salinity and humic acid on the aggregation and toxicity of polystyrene nanoplastics with different functional groups and charges[END_REF]. While these works highlight the importance and complexity of the transport mechanisms in soils, they only focused on NPT models that are not representative of the NPTs found in our environment [START_REF] Ter Halle | Nanoplastic in the North Atlantic Subtropical Gyre[END_REF]. Indeed, the composition, shape, size, and surface are known to play key roles in the transport of nanoparticles in the environment [START_REF] Hotze | Nanoparticle Aggregation: Challenges to Understanding Transport and Reactivity in the Environment[END_REF][START_REF] Pelley | Effect of particle size and natural organic matter on the migration of nano-and microscale latex particles in saturated porous media[END_REF][START_REF] Salerno | Transport of Rodlike Colloids through Packed Beds[END_REF]. Such parameters have only recently been investigated, with studies concerning UV-aged NPTs [START_REF] Liu | Aging Significantly Affects Mobility and Contaminant-mobilizing Ability of Nanoplastics in Saturated Loamy Sand[END_REF] and polystyrene NPTs containing different surface functionalities (Dong et al., 2019b).

Before investigating the influence of the soil properties, the objective of the present work is to investigate how the properties of NPTs affect their transport or accumulation in soils. Thus, different NPT models were used. All models were composed of polystyrene, since this is the only commercially available nanometer-sized plastic particle available. On the one hand, both commercial and noncommercial polystyrene latex spheres (PSL), PSL COOH and PSL COOH-P, respectively, are synthesized by a bottom-up process [START_REF] Pessoni | Soap-and metal-free polystyrene latex particles as a nanoplastic model[END_REF]. This makes them highly spherical and monodisperse and it allows their surface to be modified by the addition of carboxylate functional groups (COOH), which are common on weathered plastics [START_REF] Gewert | Pathways for degradation of plastic polymers floating in the marine environment[END_REF]. On the other hand, NPT-P are synthesized from a top-down method, which mimics the environmental mechanism of abrasion [START_REF] Hadri | Nanoplastic from mechanically degraded primary and secondary microplastics for environmental assessments[END_REF]. They have polymorphic (asymmetrical & irregular) shapes, a polydisperse size distribution, and they may contain some oxidized functional groups. The transport of NPT-P models through a sand column shows that their irregular shape will increase their probability of accumulating in porous media. These results stress the importance of choosing NPT models that have environmentally relevant physicochemical properties. μm purchased from (Polysciences© Warrington USA). These particles contain additives for stabilizing purposes, whose composition and concentration are unknown. To study the effect of particle composition, additive-free particles composed of primary (-P) polystyrene were also used. Soap and metal-free carboxylated polystyrene spheres of 430 nm (PSL COOH 430-P) were produced by emulsion polymerization according to the method described by [START_REF] Pessoni | Soap-and metal-free polystyrene latex particles as a nanoplastic model[END_REF]. Finally, NPT models with polymorphic (asymmetrical and irregular) shapes and with an average hydrodynamic diameter of 350 nm or 460 nm (NPT 350-P and NPT 460-P, respectively) were formulated from mechanically crushed PS pellets, according to the method previously described [START_REF] Hadri | Nanoplastic from mechanically degraded primary and secondary microplastics for environmental assessments[END_REF].

Methods

Dispersions of nanoplastic models

All suspensions were composed of NaCl at a concentration of 5.0 10 -3 M with pH fixed at 6.5.

Charge characterization

The zeta potential of the particles was assessed using a Wallis zetameter by Cordouan Technologies (Pessac, France). The particles were suspended in the same mobile phase as that used for the experimental conditions: NaCl at 5.0 10 -3 M and pH 6.5. The electrophoretic mobility of the colloidal particles was converted into a zeta potential by using Smoluchowski's formula. A zeta potential for sand of -50 mV was chosen according to [START_REF] Vinogradov | Measurement of streaming potential coupling coefficient in sandstones saturated with natural and artificial brines at high salinity[END_REF].

Size characterization

Hydrodynamic diameters (d H ) were determined by dynamic light scattering (DLS) using a

Vasco-Flex particle size analyzer (Cordouan Technologies). The d H of the injected NPT models and the initial sand column signal, as well as 8 mL aliquots of the eluate, were characterized by DLS. Each DLS measurement corresponds to an average of six measurements of 60 seconds each. Each sample was assigned a z-average hydrodynamic diameter (d zH ) using the cumulant algorithm. Additionally, the Sparse Bayesian learning (SBL) algorithm provides a multimodal analysis to determine the presence of several d H (see supplemental information, Figure S1). TEM images were obtained by a Jeol 2100 High Resolution Microscope (Figure S2). Data was analyzed with ImageJ software and the NanoDefine plugin to obtain the area equivalent diameters and the aspect ratios (length of major axis/length of minor axis), presented in Table 1 [START_REF] Verleysen | Evaluation of a TEM based Approach for Size Measurement of Particulate (Nano)materials[END_REF].

Transport in porous media

Transport and deposition were studied in a lab-scale sand-packed column using a liquid chromatography system (Äkta™ Pure by General Electric Healthcare). This porous media was

Theory

Colloid Filtration Theory (CFT)

Colloid filtration theory (CFT) is used to separate hydrodynamic processes from surface interactions that are occurring in the sand column. It is summarized by equation ( 1), which expresses the porous media's capacity to trap colloids. The single-collector removal efficiency (η) is the product of the single-collector contact efficiency (η 0 ) and the attachment efficiency of the particles (α) [START_REF] Petosa | Aggregation and Deposition of Engineered Nanomaterials in Aquatic Environments: Role of Physicochemical Interactions[END_REF].

= α (1)
The single-collector contact efficiency (η 0 ) is the capacity of colloids to be trapped by porous media due only to hydrodynamic mechanisms, without any favorable physicochemical conditions. It depends upon a number of factors, including the diameter of the colloids, the diameter of the sand grains in the porous media, the porosity of the sand column, and the flow rate. It was calculated according to the equation developed by [START_REF] Tufenkji | Correlation Equation for Predicting Single-Collector Efficiency in Physicochemical Filtration in Saturated Porous Media[END_REF]. The attachment efficiency (α) is a function of the physicochemical properties of the colloid dispersion and of the sand grain surfaces, such as ionic strength, surface charge and the Hamaker constant of the system. The attachment efficiency (α) is obtained according to equation (2):

= -( ) ln ( ) ⁄ (2)
where d 50 is the median diameter of the sand grains, ε and L are the sand column porosity and length, respectively, and (C/C 0 ) is the height of the breakthrough curve.

DLVO theory of colloidal stability

Surface energetics between particles and an infinite flat plane were modeled by the Derjaguin-Landau-Verwey-Overbeek (DLVO) theory, such as the total energy of interaction G tot as the sum of the Lifshitz-van der Waals G LW attraction and the repulsion due to the electronic double layer G el .

DLVO theory can be expanded upon and modified to take into account the effect of hydrophobicity, roughness and particle shape, separately. First an extended version of the DLVO theory (XDLVO)

was also modeled by adding the Lewis acid-base (hydrophopic) component G AB to the total energy of interaction. Secondly, the effect of particle shape and orientation was studied thanks to the surface element integration (SEI) method, since is what sets apart PSL COOH 430-P and both NPT-P particles.

The Lifshitz-van der Waals component G LW of the free energy of interaction between a colloid and a surface is given by:

= -[ -!" ] ( 3
)
where H is the Hamaker constant for the polystyrene-silica-water system, h is the separation distance between the colloid and the sand, and d p is the radius of the colloid. A Hamaker constant of 1.45 10 -20 was calculated from Israelachvili [START_REF] Israelachvili | Intermolecular and surface forces[END_REF] for the interaction between quartz, water and polystyrene.

The electronic double layer component G el is given by:

$% = &'( ) (* + + * -)[ . / . 0 . / 1 . 2 1 . !" 456( 7 ) 456( 7 ) + !"{1 -exp (-2>ℎ)} (4)
where ζ N and ζ S are the surface charges of the plastic colloids and of the sand surface, respectively.

The zeta potential was used instead of the surface potential. κ is the double layer thickness, which is determined by the following equation:

> = [ $ 1 A B C ∑EF G " G ] (5)
where e is the charge of the electron; ε is the permittivity of the medium, which is equal to 6.95 10 -10 C 2 .J -1 .m -1 , k B the Boltzmann constant, T the temperature, z i the valency of the ions i, and n i the number of ions i per unit volume.

The Lewis acid-base energy of interaction G AB of our system is:

HI = &J( ) KL HI M [( N )/P] ( 6
)
where λ is the correlation length, chosen as 1.65 nm, according to [START_REF] Valsesia | Direct quantification of nanoparticle surface hydrophobicity[END_REF], and h 0 is the minimum distance of separation between the particle and the surface, taken as 0.158 nm. The acidbase potential ΔG AB is expressed as:

L HI = -2 QRS + HI -RS HI T . (RS - HI -RS HI )
with US HI referring to the polar component of the surface free energy for W = water, S = sand surface and N = NPT model. RS HI is 51.00 mJ.m -2 and RS - HI is 15.00 mJ.m -2 according to Barhoumi et al.

( [START_REF] Barhoumi | Experimental Study of Thermodynamic Surface Characteristics and pH Sensitivity of Silicon Dioxide and Silicon Nitride[END_REF] The DLVO theory was modified by the SEI method described by [START_REF] Wu | DLVO Interactions of Carbon Nanotubes with Isotropic Planar Surfaces[END_REF] to understand the effect of NPT-P's nonspherical shape and their different orientations, described by the angle φ that is formed between the major axis of the particle and the collector surface [START_REF] Bhattacharjee | DLVO interaction energy between spheroidal particles and a flat surface[END_REF].

The interaction energy between a plane surface (collector) and a curved object (particle) were calculated by integrating this energy over the exact surface geometry of the object. Our particles were considered to be rod-like. The semi-major axis (L) and semi-minor axis (a) were determined according to the aspect ratio and the d zH presented in Table 1. [START_REF] Tufenkji | Breakdown of Colloid Filtration Theory: Role of the Secondary Energy Minimum and Surface Charge Heterogeneities[END_REF]. Finally, the absence of a tailing at the end of the BTCs shows that deposited particles did not detach from the porous media [START_REF] Bradford | Physical factors affecting the transport and fate of colloids in saturated porous media: FACTORS AFFECTING THE FATE OF COLLOIDS[END_REF].

Results & Discussion

Based on the BTCs, it appears that the particle shape strongly increases the deposition rate in the porous media under dynamic conditions. Surprisingly, NPT 350-P is significantly more deposited in porous media than PSL COOH 430-P, although they have a smaller initial hydrodynamic diameter and area equivalent diameter (Table 1). This confirms that the asymmetrical and irregular shape of the NPT model (NPT-P) significantly increases deposition in the porous media. It should be noted that these experiments were performed at equivalent mass concentrations instead of equivalent particle concentrations. This method was chosen to obtain a quality DLS signal at the outlet of the sand column since light scattering is proportional to the particle diameter. However, it creates a bias that favors elution of the smaller, more numerous particles (Table 1) [START_REF] Bradford | Concentration dependent transport of colloids in saturated porous media[END_REF].

Despite this, NPT 350-P, which has a higher particle concentration, was more deposited than PSL COOH 430-P. This highlights the fact that an irregular asymmetrical shape has a greater impact on the deposition rate than that of size and particle concentration.

To investigate the retention mechanisms of these NPT models, both colloidal filtration theory (CFT) and various Derjaguin-Landau-Verwey-Overbeek (DLVO) theories were studied. According to the CFT, the single-collector contact efficiency (η 0 ) increases as sub-micrometric particle sizes decrease [START_REF] Tufenkji | Correlation Equation for Predicting Single-Collector Efficiency in Physicochemical Filtration in Saturated Porous Media[END_REF]. Hence, smaller particles should be more easily retained, which is not the case in the present study. This suggests that nanoplastic models used here have considerable differences in attachment efficiency (α) and non-classical mechanisms of particle collision are operating. The attachment efficiency (α) obtained by CFT, can be used to compare the affinity of particles with sand across different experimental conditions. Caution must be taken in interpreting α since the CFT is meant to study monodisperse and spherical particles and collectors.

Since the sand grains are not monodisperse (d 60 /d 10 of 1.4) and nanoplastic models have different degrees of polydispersity (see Table 1), an average α range, namedα, and a range of α were presented in Table 2 with the corresponding heat-maps presented in supplementary information (Figure S3).

Anα of 0.0034 ± 0.0005 was obtained for PSL COOH 200, whereas PSL COOH 430-P had anα = 0.012 ± 0.0008 (Table 2). These results are in agreement with other studies, such as [START_REF] Bradford | Physical factors affecting the transport and fate of colloids in saturated porous media: FACTORS AFFECTING THE FATE OF COLLOIDS[END_REF], who found that the 450 nm carboxylated PSL attachment efficiency in 10 -3 M KCl at pH 6.8 varied from α = 0.0037 to α = 0.014 according to the polydispersity and average grain size of the quartz sand [START_REF] Bradford | Physical factors affecting the transport and fate of colloids in saturated porous media: FACTORS AFFECTING THE FATE OF COLLOIDS[END_REF]. [START_REF] Tufenkji | Breakdown of Colloid Filtration Theory: Role of the Secondary Energy Minimum and Surface Charge Heterogeneities[END_REF] found α values of 0.0058 and 0.014 for 320 nm PSL COOH in 20 10 -3 M KCl at pH 8. For the polymorphic NPT models, the average attachment efficiency is one order of magnitude greater than that for the spherical models, with α = 0.069 for NPT 350-P and 0.14 for NPT 460-P. This large increase in attachment efficiency is more likely to represent different hydrodynamic flow processes than a change in surface affinity. Indeed, the CFT theory is unsuitable to study the flow of nonspherical particles and largely overestimates α (and underestimates η 0 ) when the particle shape deviates from that of a sphere [START_REF] Salerno | Transport of Rodlike Colloids through Packed Beds[END_REF]. PSL COOH 430-P 67 ± 1 6.9 10 -3 ± 0.8 10 -3 1.2 10 -2 ± 0.1 10 -2 1.9 10 -2 ± 0.2 10 -2

NPT 350-P 28 4.0 10 -2 6.9 10 -2 1.2 10 -1

NPT 460-P 10 7.9 10 -2 1.4 10 -1 2.5 10 -1

Since the attachment efficiencies (α) obtained by the CFT theory should only reflect surface energetic interactions, they can be compared with the profiles of the interaction energy between the particles and the sand. Total surface energetic interactions (G tot ) are modeled according to DLVO theory by taking the sum of the Lifshitz-van der Waals attraction (G LW ) and the repulsion due to the electronic double layer (G el ). An extended version, the XDLVO theory, includes the Lewis acid-base (hydrophobic) interaction energy (G AB ). This component has been carefully quantified using the method described by [START_REF] Valsesia | Direct quantification of nanoparticle surface hydrophobicity[END_REF]. It should be noted, however, that since DLVO and XDLVO theories assume that the particles are spherical, they cannot correctly determine the interaction energy of the asymmetrical and irregularly shaped NPT-Ps [START_REF] Hotze | Nanoparticle Aggregation: Challenges to Understanding Transport and Reactivity in the Environment[END_REF]. To overcome this limitation, a SEI-modified DLVO calculates G LW and G el between an elongated particle and the sand surface. This modified DLVO does not, however, incorporate the hydrophobic component of the interaction energy. It is important to note that in these experimental conditions, advection prevails over diffusion, with a Peclet number varying from 112 to 240 for PSL COOH 200

and NPT 460-P, respectively. In these conditions, particle deposition can occur if the adhesive forces are stronger than the hydrodynamic forces. Therefore, the objective is only to investigate whether the 3 and full lines in Figure 2a). Based on the energy profiles, no colloid can overcome the energy barrier (ΔG tot max ) and irreversibly deposit in the primary minimum (ΔG tot min1 ). The secondary energetic minimum (ΔG tot min2 ) occurs between 20 and 60 nm from the sand surface. It varies from as low as -2.3 k B T for NPT 460-P, -2.2 k B T for PSL COOH 430-P, -1.6 k B T for NPT 460-P to -0.7 k B T for PSL COOH 200 (Table 3 and full lines in Figure 2c). For the three larger particles, reversible retention in the secondary energetic minimum can occur since the depth is larger than 1.5 k B T, which is the thermal energy of diffusion [START_REF] Israelachvili | Intermolecular and surface forces[END_REF]. 313 deposition should occur in the primary minima. The depth of the secondary energetic minimum (ΔG tot min2 ) is unchanged compared to that determined by the DLVO theory (Table 3 and Figure 2d).

The dashed lines in Figures 2a and2c show the minimal interaction energy between the sand surface and the NPT-P models, according to the SEI-modified DLVO theory. This occurs when they form an angle (φ) of π/3 with the collector surface [START_REF] Gomez-Flores | Interaction energies for hollow and solid cylinders: Role of aspect ratio and particle orientation[END_REF]. ΔG tot max is reduced to 45 k B T and 52 k B T for NPT 350-P and NPT 460-P, respectively (Table 3). The interaction energy is maximal at φ = 0 (Table 3 and Figure S4). This is mostly due to the fact the particle surface is closer to the collector [START_REF] Gomez-Flores | Interaction energies for hollow and solid cylinders: Role of aspect ratio and particle orientation[END_REF]. ΔG tot min2 is only slightly less deep than predicted by DLVO and XDLVO theories (Table 3 and Figure 2c).

The high retention rate of the NPT-P particles can be in part attributed to the fact that they may find orientations where there is no energy barrier to irreversible deposition (ΔG tot min1 ) [START_REF] Gomez-Flores | Interaction energies for hollow and solid cylinders: Role of aspect ratio and particle orientation[END_REF]. This can be expected to occur when both particle orientation and hydrophobicity are accounted for simultaneously. Furthermore, the magnitude of the total energies of interaction (ΔG tot ) of PSL COOH 430-P, NPT 350-P and NPT 460-Ps and the sand surfaces may be overestimated by the all DLVO models due to their surface roughness (Table 1 and Figure S2) [START_REF] Bhattacharjee | DLVO Interaction between Rough Surfaces[END_REF][START_REF] Hoek | Extended DLVO interactions between spherical particles and rough surfaces[END_REF]. NPT-P particles are characterized by a large range of orientationdependent interaction energies (Table 3 andS5), which will give them considerable torque. This added torque can provide the particles sufficient kinetic energy to overcome an energy barrier (ΔG tot max ) that is significantly reduced at some orientations Figure 2a [START_REF] Seymour | Transport and Retention of Colloids in Porous Media: Does Shape Really Matter?[END_REF].

However, the DLVO and XDLVO theories are inadequate to predict the deposition of plastic nanoparticles in porous media. On the one hand, the DLVO and XDLVO theories globally underestimate the amount of deposition that particles undergo. This is a well-known effect when studying the transport of particles in unfavorable (repulsive) conditions [START_REF] Elimelech | Effect of particle size on collision efficiency in the deposition of Brownian particles with electrostatic energy barriers[END_REF]. On the other hand, the height of the energy barrier is not proportional to the level of repulsion, nor is the depth of the secondary energetic minimum proportional to the level of deposition. For example, PSL COOH 200 is the most transported particle, although it has the lowest energy barrier (ΔG tot max ). This occurs because all interaction energies scale linearly with particle size [START_REF] Elimelech | Effect of particle size on collision efficiency in the deposition of Brownian particles with electrostatic energy barriers[END_REF]. In other words, larger particles should undergo more repulsion and be more transported, especially since all particles present comparable polar components of the surface free energy (γ AB ) and zeta potentials (ζ N ). However, contrary to these DLVO and XDLVO predictions, NPT 460-P, which has a higher energy barrier, is more deposited than NPT 350-P. Although a high energy barrier is present, particles may be retained by 1) attachment onto chemical heterogeneities on the colloid and collector surfaces that act as favorable sites [START_REF] Seymour | Transport and Retention of Colloids in Porous Media: Does Shape Really Matter?[END_REF][START_REF] Tufenkji | Breakdown of Colloid Filtration Theory: Role of the Secondary Energy Minimum and Surface Charge Heterogeneities[END_REF] 2) reversible attachment in the secondary energetic minimum [START_REF] Tufenkji | Breakdown of Colloid Filtration Theory: Role of the Secondary Energy Minimum and Surface Charge Heterogeneities[END_REF] and 3) straining (trapping of particles in pore throats that are too small to allow the passage of particles) [START_REF] Bradford | Physical factors affecting the transport and fate of colloids in saturated porous media: FACTORS AFFECTING THE FATE OF COLLOIDS[END_REF][START_REF] Mcdowell-Boyer | Particle transport through porous media[END_REF]. Furthermore, pore geometry creates a variety of speed profiles, such that particles may also be retained in zones with a lower flow velocity where the hydrodynamic forces are lowered [START_REF] Elimelech | Effect of particle size on collision efficiency in the deposition of Brownian particles with electrostatic energy barriers[END_REF].

Indeed, pore geometry plays a dominant role in colloid retention under conditions where a high energy barrier exists [START_REF] Tong | Excess Colloid Retention in Porous Media as a Function of Colloid Size, Fluid Velocity, and Grain Angularity[END_REF].

To determine the extent to which deposition can be attributed to chemical heterogeneity, batch adsorption experiments were conducted. The results, presented in S5, show that none of our plastic nanoparticles are significantly adsorbed onto the sand. While results from batch and column experiments can be difficult to compare when they are not normalized by α V =1, our batch experiments do show that chemical heterogeneity is far from significative [START_REF] Geitner | Measuring Nanoparticle Attachment Efficiency in Complex Systems[END_REF][START_REF] Treumann | An explanation for differences in the process of colloid adsorption in batch and column studies[END_REF]. This points to the importance of straining and adsorption in the secondary energetic minimum. after flowing through the sand column (-11%). Both spherical particles' polydispersities do not vary significantly. For the asymmetrical NPT models, the decrease in d zH after elution is globally greater and even more pronounced for the larger particles (-45% for NPT 460-P) compared to that of the smaller particles (-20% for NPT 350-P). The smaller size fraction of irregularly shaped particles, approximately 200 nm in diameter, is most mobile in the sand column. Furthermore, the polydispersity of the samples containing asymmetrical and irregular particles decreases. These changes in size point to a mechanism of physical retention of spherical particles larger than 380 nm and irregularly shaped particles larger than approximately 200 nm. This confirms that large particles are deposited by physical straining in the pores of the sand column, which is unaccounted for by the DLVO theories and underestimated by CFT [START_REF] Bradford | Physical factors affecting the transport and fate of colloids in saturated porous media: FACTORS AFFECTING THE FATE OF COLLOIDS[END_REF]. In similar experimental conditions, Weiss et al. [START_REF] Weiss | Effect of Bacterial Cell Shape on Transport of Bacteria in Porous Media[END_REF] noted a similar effect concerning the transport of 14 different bacterial strains. The bacterial cells recovered in the eluent were not only smaller but also rounder than the cells injected.

The deposition rate increases as the particle size increases and as its shape deviates from that of a sphere. In other words, smaller and smoother particles are more easily transported. This can be partially attributed to straining. Straining is expected to start having an effect when the ratio of the particle radius to the median collector radius (d p /d 50 ) exceeds 0.0017 [START_REF] Bradford | Physical factors affecting the transport and fate of colloids in saturated porous media: FACTORS AFFECTING THE FATE OF COLLOIDS[END_REF]. This is satisfied in the case of PSL COOH 430-P (d p /d 50 = 0.002). Concerning NPT-Ps , since most particles larger than approximately 200 nm are retained, straining occurs at a lower d p /d 50 ratio (Figure 3).

Straining of these particles is enhanced because the hydrodynamics of nonspherical particle flow in porous media are more complex than those of spherical particles. In highly repulsive conditions (thousands of k B T), elongated particles and are not more deposited than spherical particles since they can orient their major axis parallel to the flow [START_REF] Xu | Straining of nonspherical colloids in saturated porous media[END_REF]. However, in our conditions, where there is less repulsion (10s of k B T) nonspherical particles are increasingly deposited. This is attributed to larger collision frequencies with the porous media [START_REF] Salerno | Transport of Rodlike Colloids through Packed Beds[END_REF]. The CFT underestimates the single-collector contact efficiency (η 0 ) for spheroidal particles [START_REF] Salerno | Transport of Rodlike Colloids through Packed Beds[END_REF]. Finally, the NPT-P dispersions are highly polydisperse and as such will contain a higher fraction of large particles that are likely to be retained.

As mentioned above, reversible deposition in the secondary energetic minimum is possible for PSL COOH 430-P and both polymorphic particles NPT 350-P and NPT 460-P. The asymmetrical BTC shape of the PSL COOH 430-P is typical of either particle release from the secondary energetic minimum or blocking of the flowing particles by previously deposited particles [START_REF] Bradford | Physical factors affecting the transport and fate of colloids in saturated porous media: FACTORS AFFECTING THE FATE OF COLLOIDS[END_REF][START_REF] Bradford | Concentration dependent transport of colloids in saturated porous media[END_REF][START_REF] Redman | Bacterial Adhesion and Transport in Porous Media: Role of the Secondary Energy Minimum[END_REF][START_REF] Tufenkji | Breakdown of Colloid Filtration Theory: Role of the Secondary Energy Minimum and Surface Charge Heterogeneities[END_REF]. This was not observed for NPT 350-P and NPT 460-P, whose BTCs instead reach a plateau. It indicates either that the asymmetrical models are not detached and undergoing re-entrainment from the sand or that deposited particles do not block the incoming particles. Colloidal detachment and re-entrainment from the secondary energetic minimum mostly occur by rolling, which is easiest for spherical particles.

Despite the fact that nonspherical particles have larger hydrodynamic torques than spherical particles, as particles deviate from an ideal spherical shape, they present a higher moment arm requiring more energy to induce rolling [START_REF] Bergendahl | Prediction of colloid detachment in a model porous media: hydrodynamics[END_REF][START_REF] Seymour | Transport and Retention of Colloids in Porous Media: Does Shape Really Matter?[END_REF] The flat BTC may also indicate that deposited particles create more favorable sites for deposition by increasing the physical roughness of the porous media [START_REF] Seymour | Transport and Retention of Colloids in Porous Media: Does Shape Really Matter?[END_REF]. Finally, the shape of the NPT-P's BTCs suggest that reversible deposition not be as dominant a mechanism as straining and deposition in the primary energetic minimum.

Conclusion

This study showed that the physicochemical properties of the NPT models significantly affect their propensity to be transported. The asymmetrical and irregular shape of the environmentally relevant NPT model has a strikingly positive impact on the deposition rate in the sand column compared to that impact of particle size, composition or concentration. The deposition of these models in the porous media may be attributed partially to physical retention, as witnessed by the straining of the larger size fractions. Their deposition is also attributed to an orientation-dependent level of repulsion (ΔG tot max ) which allows the particles to take on the most thermodynamically preferable orientation, and possibly to overcome the energy barrier to irreversible deposition (ΔG tot max ), given sufficient hydodynamic torque. These results point to the importance of choosing environmentally relevant NPT models. Their asymmetrical and irregular shape, which is characteristic of the mechanical abrasion of plastic debris, increases their deposition rate in a sand column by an order of magnitude. It is clear that the relevance of environmental fate studies depends upon the quality of the proxies used [START_REF] Koelmans | Proxies for nanoplastic[END_REF]. More particularly, the impact of the NPT composition (polymer type, additives and sorbed contaminants), average size and size distribution as well as NPT shape are crucial parameters to study [START_REF] Wagner | Things we know and don't know about nanoplastic in the environment[END_REF].

The importance of the physical and chemical properties of the NPT models has been demonstrated here. Further studies are needed on the one hand, to define single-collector contact efficiency (η 0 ) for nonspherical particles and to quantify the amount of deposition that can be attributed to straining, irreversible attachment (in ΔG tot min1 ) and reversible attachment (in ΔG tot min2 ). On

  composed of Fontainebleau sand (type NE34) with a median diameter (d 50 ) of 210 μm and a uniformity index d 60 /d 10 of 1.4. The sand was dry-packed into a borosilicate column (XK 26, GE Healthcare, i.d. = 2.6 cm). The packed bed height varied from 11.2 to 12.7 cm. The porosity was 0.40 ± 0.01 on average. The eluent (mobile phase) was composed of NaCl at 5.0 10 -3 M. The pH value was fixed at 6.5 using NaOH and HCl. The sand was washed with at least 14 pore volumes (PVs) of deionized water and 22 PVs of the eluent. The pore volume was determined by an injection of NaCl at 5.0 10 -3 M. Complete tracer tests were performed with KBr. All reactants were analytical grade. The particles and tracers were injected during 6 PVs at a fixed rate of 1.57 10 -5 m.s -1 , followed by at least six pore volumes of the plastic-free suspension at the same ionic strength and pH value.At the outlet, the concentration of the NPT models eluted from the porous media was continuously measured by the UV-Vis spectrophotometer paired to the chromatography system at a wavelength λ = 226 nm, where polystyrene absorption is maximal. Breakthrough curves (BTC) were obtained by plotting the outflowing concentration of the NPT models normalized by the initial concentration as a function of pore volumes eluted. Duplicate experiments were performed for PSL COOH 200 and PSL COOH 430-P. Due to the limited yield of NPT-P dispersions and to the variability in size distributions of different NPT-P batches, duplicates were not performed.

Figure 1 :

 1 Figure 1: Breakthrough curves of KCl tracer and nanoplastic models injected continuously for 6 pore volumes at an initial concentration (C 0 ) of 5.0 10 -3 g.L -1 in 5.0 10 -3 M NaCl and pH 6.5 (error bars = standard deviation, n = 2).

Figure 1

 1 Figure 1 illustrates the breakthrough curves (BTC) of the different NPT models after transport

Table 2 :

 2 Summary of the percent in column effluent and attachment efficiencies (α) for the different nanoplastic models. The attachment efficiencies are based on average particle diameter (d p ) and collector diameters (d c ) and a range of dp and dc equal to the standard deviation (σ) of the size distributions. (Error bars = standard deviation of the duplicate experiments).

  physical and chemical properties of the particles modeled by the different DLVO theories, are correlated to the attachment efficiencies in the sand column and hence a good predictor of adhesive forces. Surface energetic interactions should not be considered an absolute prediction of deposition behavior. If the energy of interaction predicts trends in transport, the most transported particles should have the highest energetic barrier to deposition (ΔG tot max ) and lowest secondary energetic minimum (ΔG tot min2 ), where reversible deposition can occur.

Figure 2 :

 2 Figure 2: Interaction energy, scaled to k B T, between the nanoplastic models and the sand grain collector as a function of distance according to a) DLVO theory (full lines) and SEI -modified DLVO theory (dashed lines) and b) XDLVO theory. The identical data is replotted in c) and d) to highlight the secondary energetic minimum of panels a) and b) respectively.

Figure 3

 3 Figure 3: z-average hydrodynamic diameters (d zH ) of the plastic nanoparticles before and after flowing through the sand column (n=6, error bars= standard deviation).

Table 1

 1

presents the different NPT models used in the study. Carboxylated polystyrene latex spheres of 200 nm (PSL COOH 200) were Polybead® Carboxylate Orange Dyed Microspheres 0.20

Table 1 :

 1 Summary of the physicochemical properties of the nanoplastic models studied

Model name z-average diameter d zH (nm) Polydispersity Index (PDI) Area Equivalent Diameter (nm) Aspect Ratio Asperity Frequency (nm -1 ) and Amplitude (nm) Zeta potential (mV) Particle Concentra- tion (# L -1 )

	PSL COOH 200	200 ± 13	0.04 ± 0.03	193 ± 3	1.01 ± 0.01	0	-38.65 ± 2.23	1.14 10 12
						0
	PSL COOH 430-P	430 ± 28	0.07 ± 0.01	430 ± 19	1.01 ± 0.01	0.013	-29.07 ± 1.81	1.14 10 11
						12 ± 2
	NPT 350-P	350 ± 10	0.11	244 ± 133	1.67 ± 0.56			2.12 10 11
						0.013	-33.54 ± 2.72
	NPT 460-P	460 ± 25	0.19	329 ± 223	1.68 ± 0.58	11 ± 6		9.34 10 10

  .. The polar component of the free energy was measured according the method of Valsesia et al. (2018) for each of the three types of NPT models. US + HI values of 33.91, 37.47 and 31.82 mJ.m -2 were obtained for PSL COOH 200, PSL COOH 430-P and NPT-P, respectively.

Table 3 :

 3 Summary of primary maximum (ΔG tot max ) and secondary energetic minimum (ΔG tot min2 ) of interaction energies according to DLVO, SEI -modified DLVO and XDLVO theories, for the different nanoplastic models.

								SEI-modified DLVO Theory
			DLVO Theory									XDLVO Theory
							φ = 0 radians		φ = π/3 radians
	Nanoplastic	ΔG tot	max	ΔG tot	min2	ΔG tot	max	ΔG tot	min2	ΔG tot	max	ΔG tot	min2	ΔG tot	max	ΔG tot	min2
	Model	(k B T)	(k B T)	(k B T)	(k B T)	(k B T)	(k B T)	(k B T)	(k B T)
	PSL COOH 200	175	-0.7	-		-		-		-		51	-0.7
	PSL COOH 430-P	235	-2.2	-		-		-		-		74	-2.0
	NPT 350-P	244	-1.6	4.5 10 10	-0.4	45	-1.5	64	-1.6
	NPT 460-P	321	-2.3	5.9 10 10	-0.4	52	-1.8	83	-2.3
	According to DLVO theory, all particles undergo significant repulsion from the sand grains,
	with ΔG tot	max equal to 175 k

B T for PSL COOH 200, 235 k B T for PSL COOH 430-P, 244 k B T for NPT 350-P, and 321 k B T for NPT 460-P (Table

the other hand, more research is needed to implement these results by taking into account the environmental conditions of the porous media (ionic strength, natural organic matter, pH value, and soil composition). Such a complete approach is needed to fully assess the fate of NPTs in freshwater and terrestrial environments.

Competing Interest Statement

The authors declare that they have no competing interests.

Deposition of Environmentally Relevant Nanoplastic Models in

Highlights

• In order to determine the environmental fate of nanoplastics, models of nanoplastics with environmentally relevant physicochemical properties must be used.

• Nanoplastic models with polymorphic (irregular and asymmetrical) shapes are more likely to be trapped in porous media, than spherical models.

• Shape has a greater impact on the deposition rate than particle size or particle concentration.

Deposition of Environmentally Relevant Nanoplastic Models in

Declaration of interests

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:


David Montani

Julien Henry

Caroline O'connell

Xavier Jaïs

Vincent Cottin

David Launay

Gilbert Habib

Arnaud Bourdin

Mitja Jevnikar

Laurent Savale

Association between Rheumatoid Arthritis and Pulmonary Hypertension: Data from the French Pulmonary Hypertension Registry

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


Karno Schwinn

Nicolas Ferré

Miquel Huix-Rotllant
email: miquel.huixrotllant@univ-amu.fr

Efficient analytic second derivative of electrostatic embedding QM/MM energy: normal mode analysis of plant cryptochrome

Analytic second derivatives of electrostatic embedding (EE) quantum mechanics/molecular mechanics (QM/MM) energy are important for performing vibrational analysis and simulating vibrational spectra of quantum systems interacting with an environment represented as a classical electrostatic potential. The main bottleneck of EE-QM/MM second derivatives is the solution of coupled perturbed equations for each MM atom perturbation. Here, we exploit the Q-vector method [J. Chem. Phys., 151, 041102 (2019)] to workaround this bottleneck. We derive the full analytic second derivative of the EE-QM/MM energy, which allows to compute QM, MM and QM-MM Hessian blocks in an efficient and easy to implement manner. To show the capabilities of our method, we compute the normal modes for the full arabidopsis thaliana plant cryptochrome. We show that the flavin adenine dinucleotide vibrations (QM subsystem) strongly mix with protein modes. We compute approximate vibronic couplings for the lowest bright transition, from which we extract spectral densities and the homogeneous broadening of FAD absorption spectrum in protein using vibrationally resolved electronic spectrum simulations.

Introduction

Analytic derivatives of the energy revolutionized quantum chemistry, [START_REF] Pulay | Analytical derivatives, forces, force constants, molecular geometries, and related response properties in electronic structure theory[END_REF] giving access to gradients and Hessians (respectively first and second derivative with respect to nuclear positions) to explore potential energy surfaces, 2 in addition to derivatives with respect to electromagnetic fields, from which properties are nowadays routinely extracted.

The development of mixed quantum mechanics (QM) and molecular mechanics (MM) models by Warshel, Levitt and Karplus, allowed a more realistic simulation of quantum systems embedded in complex environments. [START_REF] Warshel | Theoretical studies of enzymic reactions: Dielectric, electrostatic and steric stabilization of the car-bonium ion in the reaction of lysozyme[END_REF][START_REF] Field | A combined quantum mechanical and molecular mechanical potential for molecular dynamics simulations[END_REF][START_REF] Senn | QM/MM Methods for Biolomolecular Systems[END_REF][START_REF] Nobelprize | The Nobel Prize in Chemistry 2013 -Advanced Information[END_REF] Analytic derivatives of QM/MM energy have essentially focused on analytic gradient, which have been implemented for the main self-consistent field (SCF) and post-SCF QM methods. [START_REF] Kongsted | The QM/MM approach for wavefunctions, energies and response functions within self-consistent field and coupled cluster theories[END_REF][START_REF] Sinicropi | Recent applications of a QM/MM scheme at the CASPT2//CASSCF/AMBER (or CHARMM) level of theory in photochemistry and photobiology[END_REF][START_REF] Ferré | Approximate electrostatic interaction operator for QM/MM calculations[END_REF][START_REF] Melaccio | Unique QM/MM potential energy surface exploration using microiterations[END_REF][START_REF] Zhang | Combining Quantum Mechanics and Molecular Mechanics. Some Recent Progresses in QM/MM Methods[END_REF][START_REF] Tu | Combining Quantum Mechanics and Molecular Mechanics[END_REF][START_REF] Barone | Combining Quantum Mechanics and Molecular Mechanics. Some Recent Progresses in QM/MM Methods[END_REF][START_REF] Carnimeo | Analytical gradients for MP2, double hybrid functionals, and TD-DFT with polarizable embedding described by fluctuating charges[END_REF] Analytic second derivatives of QM/MM energy are less commonly available. As a matter of fact, most of the QM/MM Hessian constructions are based on the so-called Partial Hessian approach (see 15 and references therein) in which only the QM degrees of freedom are considered. A first derivation of an analytic Hessian was introduced by Morokuma and coworkers in the framework of the substractive ONIOM QM/MM scheme. [START_REF] Dapprich | A new ONIOM implementation in Gaussian98. Part I. The calculation of energies, gradients, vibrational frequencies and electric field derivatives[END_REF] One year later, Cui and Karplus proposed an analytic Hessian evaluation in their additive QM/MM approach. [START_REF] Cui | Molecular properties from combined QM/MM methods. I. Analytical second derivative and vibrational calculations[END_REF][START_REF] Cui | Molecular Properties from Combined QM/MM Methods. 2. Chemical Shifts in Large Molecules[END_REF] The latter method is expensive, mainly due to the fact that the Hessian calculation involves solving a coupled perturbed (CP) SCF equation not only for the QM atoms but also for the MM atoms. Thus, the method of Cui and Karplus resorted to approximations of the CPSCF equations. More recently, Ghysels, Woodcock III and coworkers developed a coarse-graining method for the MM environment to reduce the number of CPSCF equations to solve. [START_REF] Ghysels | Efficient Calculation of QM/MM Frequencies with the Mobile Block Hessian[END_REF] A Hessian with fragment molecular orbital method has been derived by Nakata and coworkers. [START_REF] Nakata | Analytic second derivatives of the energy in the fragment molecular orbital method[END_REF] Other approximate second derivative methods have been proposed by Sakai and Morita, [START_REF] Sakai | Ab Initio Integrated Multi-Center Molecular Orbitals Method for Large Cluster Systems:âĂĽ Total Energy and Normal Vibration[END_REF] and Gadre and coworkers. [START_REF] Rahalkar | Enabling ab initio Hessian and frequency calculations of large molecules[END_REF] The most up-to-date formulation of the Hessian is certainly the one by Capelli and coworkers, especially derived for their very promising polarizable QM/MM method which can be used to calculate infra-red spectra of very large systems in the so-called Partial Hessian Vibrational Approach, hence ignoring the major part of the Hessian that only depends on the MM atoms. [START_REF] Giovannini | Calculation of IR Spectra with a Fully Polarizable QM/MM Approach Based on Fluctuating Charges and Fluctuating Dipoles[END_REF] Herein, we develop a new method for calculating the complete second derivatives of the QM/MM energy with additive QM-MM interaction terms. We base our derivation on singledeterminant SCF wavefunctions characterizing a QM charge density interacting with a classical point-like atomic environment via the Electrostatic Potential Fitting (ESPF) interaction Hamiltonian developed by one of us. [START_REF] Ferré | Approximate electrostatic interaction operator for QM/MM calculations[END_REF][START_REF] Melaccio | Unique QM/MM potential energy surface exploration using microiterations[END_REF] For that purpose, we apply the CPSCF Q-vector method for MM atoms, [START_REF] Schwinn | Analytic QM/MM atomic charge derivatives avoiding the scaling of coupled perturbed equations with the MM subsystem size[END_REF] avoiding one of the main bottlenecks of second derivatives of QM-MM energies. The method works even when the QM and MM regions are covalently linked via the use of the link atom method of Morokuma and coworkers. [START_REF] Dapprich | A new ONIOM implementation in Gaussian98. Part I. The calculation of energies, gradients, vibrational frequencies and electric field derivatives[END_REF] For testing our implementation, we apply this method to perform a normal mode analysis of arabidopsis thaliana cryptochrome, [START_REF] Müller | ATP Binding Turns Plant Cryptochrome Into an Efficient Natural Photoswitch[END_REF] and we compute the vibronic couplings, which allows us to efficiently construct spectral densities and account of homogeneous broadening of the spectrum via vibrationally resolved spectrum simulations.

Analytic Derivatives of the ESPF QM/MM energy

We use the notation R = (x, y, z) to symbolize the coordinates of QM atom, and R = (x, ỹ, z) to symbolize the coordinates of a MM atom.

The partial derivatives with respect to a coordinate will be written in the short-hand notation

E x = ∂E ∂x (1)
The indexes i, j, • • • will refer to occupied molecular orbitals (MO), a, b, • • • will refer to virtual MOs, p, q, • • • to generic MOs (both occupied and virtual), and µ, ν, • • • to atomic orbital (AO) index. Capital roman letters will be used for QM atom indexes.

ESPF QM/MM energy

In the ESPF method, the QM/MM energy can be written as

E = E QM + E M M (2)
in which E QM is the internal energy of the QM subsystem whose electron density is polarized by the MM subsystem (thus containing the energy contribution due to the electrostatic interaction with the MM subsystem), while E M M is the internal energy of the MM part complemented with classical interaction terms between QM and MM atoms. The electrostatic interaction energy in the ESPF method is evaluated as ∆E = Tr [Ph], in which P is the atomic density and h is the QM/MM interaction hamiltonian. [START_REF] Ferré | Approximate electrostatic interaction operator for QM/MM calculations[END_REF][START_REF] Melaccio | Unique QM/MM potential energy surface exploration using microiterations[END_REF] This interaction energy is also expressed as the classical electrostatic interaction between quantum-derived atomic multipoles q A , µ A , . . . and the MM electrostatic potential φ A , the electric field vector ∇φ A , . . .

∆E = A (q A φ A -µ A • ∇φ A + . . . ) = ∆E 0 + ∆E 1 + . . . (3)
in which the index A runs over the number of QM atoms N QM . For the sake of simplicity, we will restrict the multipolar expansion (Eq.

3) to the order 0 in the following, even though our method can be easily generalized to higher orders. [START_REF] Schwinn | Analytic QM/MM atomic charge derivatives avoiding the scaling of coupled perturbed equations with the MM subsystem size[END_REF] The expressions for the derivatives of the QM-only and MM-only energy terms will be omitted here, since they have been lengthy discussed in the literature. [START_REF] Ferré | Approximate electrostatic interaction operator for QM/MM calculations[END_REF][START_REF] Pople | Derivative studies in Hartree-Fock and Møller-Plesset theories[END_REF][START_REF] Kaledin | All-Atom Calculation of the Normal Modes of Bacteriorhodopsin Using a Sliding Block Iterative Diagonalization Method[END_REF] Here, we will focus on the electrostatic interaction term and its derivatives.

At the order 0 of the QM density of charges, each atomic charge q A is expressed as the sum of the nuclear charge Z A and the mean value of an atomic one-electron operator, whose matrix elements are denoted Q A,µν . The interaction energy in atomic orbital basis reads,

∆E = A,µν (-P µν Q A,µν + Z A ) φ A (4) = A q A φ A = µν P µν h µν + A Z A φ A
in which we have introduced a new one-electron operator whose elements are:,

h µν = - A Q A,µν φ A (5)
In a nutshell, the ESPF method consists in fitting the Q A,µν matrix elements to reproduce the purely electronic QM electrostatic integrals calculated on a grid surrounding the QM subsystem. 9

First derivative

The analytic first derivatives of the ESPF interaction term were first derived in Ref. 9. For completeness, we reexpress them with the current notation. The first derivative of Eq. 4 can be expressed as a vector gradient,

∇[∆E] = ∆E x ∆E x (6)
in which the first block corresponds to the derivative with respect to QM atom positions and the second block to MM atoms. The expression of the energy gradient with respect to QM coordinates is given by

∆E x = µν P x µν h µν + µν P µν h x µν + A Z A φ x A (7)
As usual, the first term is not calculated thanks to the introduction of the so-called energyweighted density matrix -PhPS x . [START_REF] Ferré | Approximate electrostatic interaction operator for QM/MM calculations[END_REF][START_REF] Melaccio | Unique QM/MM potential energy surface exploration using microiterations[END_REF][START_REF] Pople | Derivative studies in Hartree-Fock and Moller-Plesset theories[END_REF] The derivative of the ESPF operator element h µν is given in Appendix 5.1. Regarding the derivative with respect to the MM coordinates, the energy gradient can be expressed as

∆E x = A q A φ x A , (8)
noting that S x = 0 for MM perturbations.

Second derivative

The analytic second derivative of Eq. 4 can be expressed as a Hessian matrix (10) The second derivative of the Fock operator requires the second derivative of the ESPF operator, which is given in 5.1. This second derivative block also requires the density matrix derivative P y with respect to QM atom coordinates, which is obtained by solving a CPSCF equation for each QM coordinate. [START_REF] Mcweeny | Some Recent Advances in Density Matrix Theory[END_REF] The ∆E xỹ block over a QM and an MM atom coordinates is given by

∇ 2 [∆E] = ∆E xy ∆E xỹ ∆E xy ∆E xỹ (9)
∆E xỹ = A Z A φ xỹ A + µν P µν h xỹ µν + P ỹ µν h x µν (11)
The explicit form of the second derivative of the ESPF operator is given in Appendix 5.1. This energy derivative requires the explicit construction of P x. These would require the solution of a CPSCF equation for all MM atom coordinates, which is prohibitive. A more convenient form to implement, is to compute the equivalent ∆E xy derivative block, which is given by

∆E xy = A q y A φ x A + q A φ xy A ( 12
)
in which the charge derivative is given by

q y A = - µν P y µν Q A,µν + P µν Q y A,µν (13)
In this case, only derivatives of the density matrix with respect to QM coordinates are required, which is already constructed for the QM energy second derivative. [START_REF] Pople | Derivative studies in Hartree-Fock and Møller-Plesset theories[END_REF] Finally, the ∆E xỹ derivatives over MM atom coordinates is given by

∆E xỹ = A q ỹ A φ x A + q A φ xỹ A ( 14
)
requires the derivatives of the ESPF atomic charges, which would in principle require the density derivatives

q ỹ A = - µν P ỹ µν Q A,µν (15)
The construction of P ỹ would imply solving a coupled-perturbed type equation for each MM perturbation. Usually, the number of MM atoms is hundreds to thousands, becoming rapidly intractable. In Ref. 24, we showed how to avoid this by defining an auxiliary set of coupled-perturbed equations (Q-vector method) that only scales with the number of QM atoms, and from which we can efficiently reconstruct the charge derivatives,

q ỹ A = - B,ia Q A,ia φ ỹ B QB,ia . (16)
where QB,ia is the solution of the Q-vector coupled-perturbed equation, [START_REF] Schwinn | Analytic QM/MM atomic charge derivatives avoiding the scaling of coupled perturbed equations with the MM subsystem size[END_REF] which allows to rewrite the MM contribution to second derivative as

∆E xỹ = A q A φ xỹ A + B,ia Q A,ia QB,ia φ ỹ B φ x A (17)
This equation can be interpreted in terms of classical atomistic mechanics,

∆E xỹ = A q A φ xỹ A + B Q A Q AB φ ỹ B φ x A ( 18
)
in which

Q A = -T r[PQ A ],
and

Q AB = ia QB,ia Q A,ia Q A (19)
Then, Q AB can be interpreted as a tensor that, given an electric charge Q A on QM atom A interacting with a gradient of the external potential, will give the "effective" charge on QM center B that interacts with another gradient of the external field.

Link atoms

For molecular systems in which there exist covalent bonds at the interface between the MM and QM subsystems, unphysical free valences are created in the QM side. Morokuma and coworkers proposed to use a fictitious link atom in order to saturate these dangling bonds, most often using a hydrogen atom. The link atom is constrained to remain positioned between the frontier MM and QM atoms, 16

R link = R QM + g RMM -R QM (20)
in which g is a parameter that fixes the coordinates of the link atom between the covalently bonded QM and MM atoms.

From the QM perspective, the "auxiliary" link atom is considered as a QM atom, and therefore, N QM includes all real QM atoms plus the link atoms. Consequently, there will be an extra contribution for the QM and MM frontier atoms that comes from the link atom. In the case of first energy derivatives, this is given by

E x | link = E x + (1 -g) ∂E ∂x link E x| link = E x + g ∂E ∂ xlink (21)
In other words, the link atom gradient is projected onto the frontier QM and MM atoms.

Similarly, the contribution of the link atom to the MM and QM second derivatives can be expressed as

E xy | link = E xy + (1 -g) 2 ∂ 2 E ∂x link ∂y link + (1 -g) ∂ 2 E ∂x∂y link + ∂ 2 E ∂x link ∂y E xỹ | link = E xỹ + g(1 -g) ∂ 2 E ∂x link ∂y link + (1 -g) ∂ 2 E ∂x link ∂ ỹ + g ∂ 2 E ∂x∂y link E xy | link = E xy + g(1 -g) ∂ 2 E ∂x link ∂y link + g ∂ 2 E ∂x link ∂y + (1 -g) ∂ 2 E ∂ x∂y link E xỹ | link = E xỹ + g 2 ∂ 2 E ∂x link ∂y link + g ∂ 2 E ∂ x∂y link + ∂ 2 E ∂x link ∂ ỹ (22)

Computational Details

All electronic structure computations have been performed with a modified version of Gaus-sian16 including the Q-vector method and the ESPF energy, gradient and Hessian. The MM information is obtained via an interface with a modified version of Tinker 6.3.3. [START_REF] Ferré | Approximate electrostatic interaction operator for QM/MM calculations[END_REF][START_REF] Frisch | sianËIJ16 Revision A.01[END_REF][START_REF] Ponder | Software tools for molecular design[END_REF] The vibrationally resolved spectrum for FAD has been obtained by Fourier transforming the time-dependent auto-correlation function for the multi-dimensional harmonic oscillator, as implemented in Gaussian16. [START_REF] Frisch | sianËIJ16 Revision A.01[END_REF][START_REF] Baiardi | General Time Dependent Approach to Vibronic Spectroscopy Including Franck-Condon, Herzberg-Teller, and Duschinsky Effects[END_REF] The Franck-Condon approximation of the auto-correlation function has been integrated on 2 [START_REF] Nakata | Analytic second derivatives of the energy in the fragment molecular orbital method[END_REF] time steps, over an interval of about 10.5 ps. A Gaussian broadening of 235 cm -1 has been applied to the vibrational peaks.

All QM calculations have been performed at the DFT and full TDDFT levels, [START_REF] Casida | Progress in Time-Dependent Density-Functional Theory[END_REF] using the B3LYP functional [START_REF] Becke | Density-functional thermochemistry. III. The role of exact exchange[END_REF][START_REF] Lee | Development of the Colle-Salvetti correlationenergy formula into a functional of the electron density[END_REF][START_REF] Vosko | Accurate spin-dependent electron liquid correlation energies for local spin density calculations: a critical analysis[END_REF] and the TZVP/TZVPFit basis set. [START_REF] Schäfer | Fully optimized contracted Gaussian basis sets for atoms Li to Kr[END_REF][START_REF] Schäfer | Fully optimized contracted Gaussian basis sets of triple zeta valence quality for atoms Li to Kr[END_REF][START_REF] Eichkorn | Auxiliary basis sets to approximate Coulomb potentials[END_REF][START_REF] Eichkorn | Auxiliary basis sets for main row atoms and transition metals and their use to approximate Coulomb potentials[END_REF] For the MM calculations, we have employed the AMBER99 force field. [START_REF] Case | AMBER 2018[END_REF] Van der Waals parameters for FAD have been taken from Ref. 42.

Results and Discussion

As a first example of our method, we have computed the frequencies of flavin adenine dinucleotide at the QM/MM level. In the QM region, we have placed the isoalloxazine ring while the rest of the molecule is considered MM. This selection requires the introduction of a link atom to separate the QM and MM regions, which has been placed 3 bonds far away from the isoalloxazine ring to avoid overpolarization, as shown in Fig. 1. In this figure, we show the frequencies computed at the QM/MM level compared to the QM frequencies, which we take as reference. For the sake of comparison, we also show the pure MM frequencies. Overall, the quality of the QM/MM and MM frequencies is in very good accordance with the QM frequencies. This is expected since the MM parameters have been obtained from purely ab initio data. Some appreciable errors are observed for the highest frequency modes, as well as the mid-range. The QM/MM frequencies follow a similar trend than the MM frequencies, indicating that these modes have simply a small mixing and remain rather purely MM or purely QM. Indeed, the most delocalized frequencies are found only in the low frequency range (see further below). This results indicate that the choice of force field is indeed critical for obtaining a good spectrum of frequencies that is comparable to the ab initio spectrum.

To further illustrate the potential of our method, we have performed a normal mode analysis of plant cryptochrome, including vi- bronic coupling and homogeneous broadening of absorption spectrum using vibrationally resolved spectra simulations. For this, we have selected chain A of arabidopsis thaliana cryptochrome (PDB code 2J4D, see Fig. 2). Crystallographic water molecules have been kept in the model. Hydrogens have been added for standard aminoacid protonation states at physiological pH. The entire FAD coenzyme was defined as the QM subsystem. FAD is not covalently linked to the rest of the protein, thus not requiring any special care of the QM/MM boundary. The initial PDB structure was pre-optimized at the Amber99 level using microiterations. Subsequent geometry optimization, hessian calculation and normal mode analysis were carried out at B3LYP/TZVP/TZVPFit/Amber99 level for the entire protein model. One of the main advantages of the method presented here is that the CPSCF equations for the MM perturbations are not scaling with the MM subsystem size. As we showed, we can decouple this dependence by solving the Q-vector equations for N QM atoms and reconstruct the atomic charge derivatives with respect to MM perturbations by contracting with the external field derivatives. This induces an obvious saving in computational time. For the present case (plant cryptochrome, 9432 atoms at the B3LYP/TZVP/TZVPFit/Amber99 level using 32 cores of 2.6 GHz), the solution of conventional CPSCF equations for the full system would take around 86.5 hours, while the Qvector method (84 atom) just takes 17 minutes. This allows the construction of the full Hessian including the charge polarization terms at the cost of a mechanical embedding QM/MM full Hessian.

Normal mode analysis

We performed the normal mode analysis of the plant cryptochrome shown in Fig. 2 at the minimum energy structure of the QM/MM internal energy. This geometry was obtained from a quadratically convergent optimization of the QM subsystem, with a steepest descent optimization of the MM susbsystem at every step of optimization (microiterations). At the minimum structure, the FAD remains in a Ufolded conformation due to the electrostatic interactions with the protein environment. The protein model has 8121 atoms, from which 84 atoms are in the QM subsystem. This leads to a total of 24357 normal modes. In Fig. 3, we analyse the degree of delocalization as well as the QM character of these modes. First, we computed the participation ratio in order to measure the degree of delocalization. For this, we the participation ratio as

P R(ω i ) = 3N • 3N j ( 0 L † i,j • 0 L j,i ) 2 -1 (23)
in which 0 L i is the mass-weighted normal mode displacement of the ground state. For a fully delocalized mode, all atoms contribute equally to the mass-weighted displacement, and the participation ratio is equal to one, whereas if only few atoms contribute to the displacement, the participation ratio has the value of (3N ) -1 (which is close to 0 in the case of plant cryptochrome.) Following this analysis, we can observe that the low frequency modes (<200 cm -1 ) show on average the largest degree of delocalization, involving thus a large number of atoms in each normal mode frequency. These modes are usually involved in describing protein phase transitions. The rest of the normal modes involve only a few atoms.

The participation ratio analysis is not sufficient to determine the degree of mixing between QM and MM atom motions. If one would perform a QM subsystem normal mode analysis, one would obtain 246 full QM normal modes. For the full protein diagonalization, these modes are strongly mixed. In Fig. 3, we report the QM character of each normal mode, defined as

QM (ω i ) = 3N j∈QM 0 L † i,j 2 . ( 24
)
As we can observe, the low frequency modes below 200 cm -1 have a negligible QM character. This implies a strong mixing between the QM and MM atoms in this region. Since the number of MM atoms is so large, the QM character is "diluted" within the MM motion. Only two normal modes between 1200-1400 cm -1 keep a clear QM character (>90%), and another 23 normal modes have a QM character larger than 60%, comprised in the range of 400-1600 cm -1 . Therefore, only 10% of the FAD normal modes keep in essence its original character. For the rest, QM and MM atom motions are essentially mixed, even if the modes are not fully delocalized over the whole protein. Importantly, this mixing is present in the whole range of frequencies, and it is not only appearing for the fully delocalized low-frequency modes.

Vibronic couplings, Huang-Rhys factors and spectral densities

The Huang-Rhys couplings are a fundamental measure for the electron-nuclear (vibronic) coupling. The HR factor is directly related to the reorganization energy of one mode in the excited state. They can be defined as

HR i,I = 1 2 ω i ∆q 2 i,I , (25)
in which ω i is the normal mode frequency and

∆q i,I = 0 L i M 1/2 (R I -R 0 )
is the projection in the ground state normal mode direction 0 L i of the Cartesian geometry difference between excited (R I ) and ground state (R 0 ) minima, mass weighted by the mass metric M = δ ij m i , in which m i is the atomic mass of atom i.

The ∆q I i can be approximately obtained by the so-called vertical gradient approximation, in which the mass-weighted energy gradient of the excited state I (∆q I i ) is projected to the normal mode direction, leading to

∆q i,I = - (∇E I | 0 L i ) ω 2 i ( 26
)
In order to define the effect of the bath modes, the notion of spectral density J(ω) is useful. It can be defined as

J I (ω) = π i HR i,I ω 2 i δ(ω -ω i ) (27)
Indeed, the spectral density contains information about each normal mode frequency, and its coupling to the system, which suffices to characterize the internal dissipation of an excited state due to nuclear motion. Using the information of the QM/MM analytic second derivatives together with the excited-state gradient at the ground-state minimum, we can construct the spectral densities characterizing the same cryptochrome model. In Figure 4, we have plotted such spectral density for the lowest energy bright excited state. Two spectral densities are shown, one corresponding to the normal modes and gradient of the full protein, and the other corresponding to the normal modes and gradient of the QM subsystem.

The overall contribution of MM atoms to the broadened spectral density is only minor. This is expected because the transition density is mainly located on the chromophore, which is entirely located in the QM subsystem, and thus only modes with substantial QM motion should be contributing to the spectral density. Still, an indirect contribution of the MM atoms appear in the spectral density due to the mode mixing, as shown for example in the most intense peak around 1620-1625 cm -1 . While the QM spectral density has only one peak of 40 meV, the full protein spectral density shows a spectral density contribution of 4 modes. This is due to the mixing of QM and MM atom motions, alike to the intensity borrowing mechanism of electronic spectra, in which dark states couple to bright states borrowing dipole intensity and thus becoming observable. This can be observed for almost all peaks in the spectral density, which reflects once again the strong normal mode mixed character that we showed in Fig. 3. It is important to notice though that even though there is a strong mixing, the spec-tral densities of the QM and the full protein are essentially equivalent. This is because the normal mode mixing occurs only between QM and MM modes at almost the same frequency. Eventually, the effect of mixing has only a minor impact on the frequency of the modes (less than 1 cm -1 for the most intense peak). This results justify the extraction of vibronic couplings from the normal mode analysis of the QM subsystem.

Vibrationally resolved absorption spectrum

Flavoproteins like cryptochrome act as bluelight sensors in biological systems. [START_REF] Conrad | Photochemistry of flavoprotein light sensors[END_REF] These proteins have a well-resolved vibrational progression of the absorption spectrum, characteristic of the isoalloxazine chromophore. Unlike other types of proteins in which inhomogeneous broadening dominates, the homogeneous broadening has to be explicitly taken into account to reproduce the color of absorption of the flavin chromophore. [START_REF] Kabir | Electronic spectra of flavin in different redox and protonation states: a computational perspective on the effect of the electrostatic environment[END_REF][START_REF] Mondal | Impact of the redox state of flavin chromophores on the UVâĂŞvis spectra, redox and acidity constants and electron affinities[END_REF] The approximate vibronic couplings shown in Eq. 26 can be used to compute the absorption spectrum including the vibrational structure due to homogeneous broadening.

Here, we compute the vibrational resolved spectrum of FAD in plant cryptochrome, by using the QM subsystem normal modes in conjunction with the vibrationally-resolved timedependent approach of Bloino and coworkers. [START_REF] Baiardi | General Time Dependent Approach to Vibronic Spectroscopy Including Franck-Condon, Herzberg-Teller, and Duschinsky Effects[END_REF] In Fig. 5, we compare the experimental and theoretical absorption spectra for arabidopsis thaliana plant cryptochrome with FAD in its fully oxidized form. [START_REF] Müller | ATP Binding Turns Plant Cryptochrome Into an Efficient Natural Photoswitch[END_REF] The experimental absorption spectrum is characterized by two main absorption peaks of π → π * character. The lowest bright transition has a peak maximum around 450 nm and the second bright state at around 360 nm. [START_REF] Müller | ATP Binding Turns Plant Cryptochrome Into an Efficient Natural Photoswitch[END_REF] These peaks have a marked vibrational structure, especially due to in-plane C-N stretchings. [START_REF] Klaumünzer | DFT Calculation of Vibrational and Vibronic Spectra of Riboflavin in Solution[END_REF] Indeed, the first electronic transition is distributed in resolved peaks, approximately found at 470, 450 and 420 nm. Our simulated spectrum perfectly matches the relative position and intensities of the experimental peaks. From our simulations, we observe that these peaks are in fact formed from two electronic transitions, a HOMO→LUMO transition which give the basic vibrational structure, and a HOMO-1→LUMO transition which is centered at 420 nm. For the second absorption peak, we have a reasonable agreement between simulated and experimental spectra. Indeed, the relative peak position and the relative intensities between the two vibronic peaks are well reproduced by the theory. However, our theoretical simulations underestimate by around 40% the relative intensity with respect to the lowest energy absorption peak. We attribute this either to the approximations of the level of theory or some missing contributions of low-intensity high-energy peak. with the experimental spectrum (black) of arabidopsis thaliana. The grey dotted lines indicate the individual contributions from each electronic state. For the sake of comparison, both spectra have been normalized for the highest intensity peak between the range of 350 nm and 500 nm. The simulated spectrum has been red-shifted by 12 nm to match the experimental maximum. Experimental spectrum has been taken from Ref. 25

Conclusion

We have derived a fully analytic expression for the second derivative of the electrostatic embedding ESPF QM/MM energy. In this approach, the atomic charges representing the QM charge distribution are not only polarizable by the usual QM/MM electrostatic interactions, but also by the motion of the MM atoms. QM/MM electrostatic energy second derivative contains MM-block contributions too expensive to compute, and usually simply ignored or approximated. Using the ESPF method, we have shown that these terms can be equivalently computed with a set of auxiliary coupledperturbed self-consistent field equations that only scales with the QM subsystem size, [START_REF] Schwinn | Analytic QM/MM atomic charge derivatives avoiding the scaling of coupled perturbed equations with the MM subsystem size[END_REF] reducing thus the total computational cost of the full electrostatic embedding QM/MM hessian to a similar cost than mechanical embedding QM/MM hessian.

To show the capabilities of our method, we have computed the full hessian of arabidopsis thaliana plant cryptochrome. We have evidenced that most of the FAD chromophore (QM) vibrational modes are strongly mixed with the modes of the protein, and only a few high-frequency modes around 1200-1500 cm -1 keep a full QM character. In addition, we have computed the vibronic couplings and spectral densities for the bright excited state, together with the vibrationally resolved spectra of FAD in the protein. These applications open up the use of QM/MM hessians construct vibronic models for photoactive proteins and QM systems in complex environments.

Several limitations are currently imposed in our current model for a practical use. First, one of the main limitations of our QM/MM method is the lack of periodic boundary conditions to include solvent effects. In the future, we plan to implement our method in an extension of the ESPF QM/MM method to include Ewald summations. [START_REF] Holden | Analytic gradient for the QM/MM-Ewald method using charges derived from the electrostatic potential: Theory, implementation, and application to ab initio molecular dynamics simulation of the aque-ous electron[END_REF] Meantime, the only strategy is to freeze a state of the solvent, and relax the protein within a solvation shell. Since the QM equations do not scale with the MM subsystem size, the calculations considering a large number of solvent molecules should be tractable at a reasonable computational cost. Second, the current simulation of vibrationally resolved absorption spectrum requires a positive semi-definite Hessian matrix. This implies that any initial geometry must be optimized to a minimum energy structure so that all frequencies are positive. The geometry optimization step is time-consuming, hence severely limiting the usual efficiency associated with the inhomogenous broadening as it is frequently reported in the literature, i.e., by selection of snapshots from a MD simulation and performing the QM/MM analysis at each of them. This limitation could be overcome by using for example instantaneous normal modes, [START_REF] Kalstein | Vibrational dynamics of polyatomic molecules in solution: assignment, time evolution and mixing of instantaneous normal modes[END_REF] although further investigation in this direction is required.

Finally, we plan to extend our model to address computed vibrational spectroscopies. Such simulations face frequently the problem of comparing with experimental spectra, which could ultimately be attributed to the lack of anharmonic effects, both from the QM and the MM frequency spectrum. To better match the QM frequencies to experimental data, it is common to scale the Hessian via a scaling factor that depends on the method and the basis set. Morokuma and coworkers applied a similar strategy to their ONIOM QM/MM Hessian, in which three different scaling factors for the QM, MM and QM-MM Hessian blocks are applied. [START_REF] Dapprich | A new ONIOM implementation in Gaussian98. Part I. The calculation of energies, gradients, vibrational frequencies and electric field derivatives[END_REF] Such scheme would be readily applied to our analytic QM/MM Hessian too. In addition, some recent developments at the level of force fields for better describing infra-red spectra can also be trivially integrated in our Hessian scheme. [START_REF] Nutt | Theoretical Investigation of Infrared Spectra and Pocket Dynamics of Photodissociated Carbonmonoxy Myoglobin[END_REF][START_REF] Mackie | The anharmonic quartic force field infrared spectra of three polycyclic aromatic hydrocarbons: Naphthalene, anthracene, and tetracene[END_REF]

Appendix: Derivatives of the ESPF operator

The ESPF atomic charge operator elements Q A,µν are obtained by fitting the corresponding electrostatic potential integrals V k,µν ,

V k,µν = drχ µ (r)|r -R k | -1 χ ν (r) (28)
calculated over a grid of N grid points k, at position R k , around the QM atoms:

A T A,k Q A,µν = V k,µν (29)
The kernel T is a rectangular matrix of dimensions N grid •N QM . The kernel definition depends on the multipolar expansion order. For example, for atomic charges it contains the inverse distance between QM atom and the grid point

T A,k = |R A -R k | -1
. Thus, we can give an expression for the Q A,µν = χ µ | Q A |χ ν atomic charge operator elements by solving the previous set of linear equations, leading to

Q A,µν = N grid k T † T -1 T † A,k V k,µν (30)
The derivative of the ESPF hamiltonian involves the derivative of the atomic charge operator. A recursive formula for the n-derivative of the charge operator can be obtained by directly deriving n-times Eq. 29,

Q n A,µν = N grid k T † T -1 T † A,k F n k,µν (31)
in which

F n k,µν = V n k,µν - n-1 i=0 N QM B T n-i k,B Q i B,µν (32)
Here we give the explicit first and second derivatives of the h operator defined in Eq. 5. The first derivative of the operator is given by

h x µν = A Q x A,µν φ A + Q A,µν φ x A (33)
The second term involves the derivative of the classical external field, which have simple expression given the form of the force field. The derivative of the atomic charge operator can be in term expanded (in matrix notation) as

Q x = T † T -1 T † [V x -T x Q] (34)
we can easily see that the computation of firstorder derivatives of h operators requires the computation of φ x A , V x k,µν and T x in addition to the zero-order field and charges.

The second derivatives of the h operator is given by

h xy µν = A Q xy A,µν φ A + Q A,µν φ xy A (35) + A Q x A,µν φ y A + Q y A,µν φ x A
The second term involves the second derivative of the classical external field, which again has a simple expression given the form of the force field, while the third and fourth term are constructed with elements of the first derivative of the operator. Hence, one needs only to develop the second derivative of the atomic charge operator. This is given by

Q xy = T † T -1 T † [V xy -T xy Q] -T † T -1 T † [T x Q y + T y Q x ] (36)
Therefore, in order to obtain the second derivative of the ESPF operator, we need to compute φ xy A , V xy k,µν and T xy in addition to the first derivatives and the zero-order quantities.

Figure 1 :

 1 Figure1: Comparison of the QM/MM and MM frequencies (y-axis) with respect to full QM frequencies (x-axis) of flavin adenine dinucleotide. For the computation of frequencies, a geometry optimization has been performed at each level of theory (QM: B3LYP/6-31G*, MM: Amber99; QM/MM: B3LYP/6-31G*/Amber99). For the QM/MM calculation, a link atom has been placed between the isoalloxazine ring (treated at the QM level) and the rest of the molecule as indicated in the graph.

Figure 2 :

 2 Figure 2: Representation of the arabidopsis thaliana cryptochrome (PDBID: 2J4D) used in the QM/MM simulations. The FAD fragment (ball and stick) constitutes the QM subsystem, while the protein (ribbon) is treated at the MM level. Crystallographic water molecules have been excluded from the plot, but included in the simulations.

Figure 3 :

 3 Figure 3: Normal mode analysis of plant cryptochrome. (a) participation ratio showing the degree of delocalization of normal modes, and (b) analysis of the character QM character. Modes with less than 10% QM atom displacement are classified as MM (black) while modes with more than 90% are classified as QM (yellow).

Figure 4 :

 4 Figure 4: Stick spectrum corresponding to the spectral density for the first bright excited state. The continuous spectral density has been obtained from a Lorenzian broadening of 500 cm -1 . The red line correspond to the normal modes of the full protein, while the black line to the QM subsystem normal modes. The zoom-in of the most intense peak is shown in the righthand inset.

Figure 5 :

 5 Figure5: Comparison of the simulated vibrationally resolved absorption spectrum (red) with the experimental spectrum (black) of arabidopsis thaliana. The grey dotted lines indicate the individual contributions from each electronic state. For the sake of comparison, both spectra have been normalized for the highest intensity peak between the range of 350 nm and 500 nm. The simulated spectrum has been red-shifted by 12 nm to match the experimental maximum. Experimental spectrum has been taken from Ref.25

Acknowledgement The authors acknowledge financial support by the "Agence Nationale pour la Recherche" through the project BIO-MAGNET (ANR-16CE29-0008-01). "Centre de Calcul Intensif d'Aix-Marseille" is acknowledged for granting access to its high performance computing resources.

Graphical TOC Entry

Electrostatic Embedding

Full analytic Hessian Normal mode analysis Homogeneous broadening