Camille Banoun

Flavien Bourdreux

Emmanuel Magnier

Guillaume Dagousset
email: guillaume.dagousset@uvsq.fr

Intermolecular C-O Bond Formation with Alkoxyl Radicals: Photoredox-Catalyzed -Alkoxylation of Carbonyl Compounds

Due to the high reactivity of alkoxyl (RO•) radicals and their propensity to easily undergo -scission or Hydrogen Atom Transfer (HAT) reactions, intermolecular alkoxylations involving RO• radicals are barely described. We report herein for the first time the efficient intermolecular trapping of alkoxyl radicals by silyl enol ethers. This photoredox-mediated protocol enables the introduction of both structurally simple and more complex alkoxy groups into a wide range of ketones and amides.

Alkoxyl radicals (RO•) are key intermediates implied in many natural processes, ranging from atmospheric chemistry to biology. 1 Due to their high reactivity, they hold a special place in the field of radical chemistry, quite different from their alkyl analogs. Notably, their development as a reactive species in organic chemistry is mainly limited to their high propensity to undergo HAT or -scission reactions. 2 These reactivities usually lead to the generation of a more stable C-centered radical, which is then engaged in C-C or C-Heteroatom bond formation. A third reactivity of RO• radicals consists in their ability to directly add onto alkenes, thus involving the formation of a C-O bond. 3 However, these alkoxylation reactions are mostly performed in an intramolecular fashion, in order to be competitive with fast HAT/scission pathways. [START_REF] Beckwith | Generation of Alkoxy Radicals from N-Alkoxypyridinethiones[END_REF][START_REF] Kim | Visible-Light-Induced Cascade Radical Ring-Closure and Pyridylation for the Synthesis of Tetrahydrofurans[END_REF] Conversely, the development of alkoxylation reactions involving a highly challenging intermolecular RO• addition onto alkenes or other unsaturated systems is still at its infancy. [START_REF] Barthelemy | Alkoxyl Radicals Generated under Photoredox Catalysis: A Strategy for anti-Markovnikov Alkoxylation Reactions[END_REF] In this article we report the first intermolecular addition of alkoxyl radicals to silyl enol ethers. This methodology leads to the synthesis -alkoxylated carbonyls, which are versatile synthetic compounds used in many fields of application. 7 In spite of this wide interest, the C-H alkoxylation in  position to carbonyl compounds has only been scarcely studied so far although its development would be an important added-value over existing methods. Indeed, as no "RO + " electrophile is readily available, [START_REF] Rozen | Hypervalent Iodine Oxidation of Silyl Enol Ethers under Lewis-Acid Conditions in Methanol. A General Route to -Methoxy Ketones[END_REF] a classical strategy which would involve the Scheme 1. Strategies for C-H -alkoxylation of carbonyl compounds.

activation of the carbonyl substrate into a nucleophilic intermediate (enolate, silyl enol ether, enamine…) is unadapted. To date, only "umpolung" approaches have been developed: the ketone substrate is oxidized and further reacts with a nucleophilic alcohol ROH (Scheme 1a). 9 Such umpolung strategy has however two major drawbacks: i) excess amounts of strong oxidants (mCPBA, Pb(OAc)4, PhI(OAc)2, tert-butyl peroxide…) are required, which leads to poor functional-group compatibility and/or to substantial amounts of by-products; ii) the corresponding alcohol is always used as solvent or co-solvent, which limits the application scope to structurally simple alkoxy groups (mostly OMe and OEt).

Our new radical-based -alkoxylation strategy overcomes these issues. Indeed, given the recent successful developments using photoredox catalysis, 2 our approach is based on the mild photocatalytic generation of RO• radicals, and involves only stoichiometric amounts of the precursor of both structurally simple and more complex alkoxy moieties. For the success of our methodology, we hypothesized that silyl enol ethers I, easily generated from the corresponding carbonyl compounds, would be suitable substrates to efficiently trap RO• radicals in an intermolecular manner (Scheme 1b). The resulting -silyloxy C-centered radical II would be highly stabilized and other HAT/scission pathways would then be minimized.

Within this framework, we started our study by examining the alkoxylation of silyl enol ether 1a in presence of Ir(ppy)3 3a as the photocatalyst (1 mol%, Table 1). Based on our previous work, N-methoxypyridinium salt 2a was a General conditions: 1a (0.2 mmol), 2a (x equiv), 3 (1 mol%) and additive (2 equiv) in 2 mL of MeCN irradiated with blue LEDs at rt for 18 h. b Yields determined by 19 F NMR spectroscopy using fluorobenzene as an internal standard. c Yield into brackets refers to chromatographically pure product. d In the dark.

selected as the source of MeO• radical. 6a,10 We were pleased to see that when this mixture was irradiated with blue LEDs in MeCN, the expected -methoxylated ketone 4a was obtained in 27% yield (Table 1, entry 1). However, 4'fluoroacetophenone, resulting from direct hydrolysis of 1a, was the major product of this reaction. To prevent the formation of this by-product, a series of organic and inorganic bases were screened (entries 2-9), in which K2HPO4 afforded the best result (61% yield, entry 5). Further screening of other parameters (stoichiometry and photocatalyst, entries 10-14) led to the optimal conditions (1.8 equiv of 2a, with 1 mol% of Ru(bpy)3(PF6)2 3d, entry 14), giving product 4a in 89% isolated yield. Importantly, no trace of by-products resulting from a potential HAT side-reaction could be detected, pointing out the ability for silyl enol ethers to efficiently trap RO• radicals in an intermolecular manner. In addition, no product was observed in the absence of either photocatalyst or light (entries 15-16).

With these optimal reaction conditions in hand, the scope of this novel protocol was then evaluated (Scheme 2). First, various silyl enol ethers 1a-1j derived from aromatic ketones bearing electron-donating (Me, OMe) or electronwithdrawing groups (F, Cl, Br, CO2Me, CF3, CN) in ortho, meta or para positions reacted all smoothly to give the corresponding -methoxylated ketones 4a-4j in 57-89% yield. A gram-scale experiment was also successfully carried out and provided 1.48 g of ketone 1a without any loss of efficiency. Heteroaryl substrates 1k-m with an indole, benzofuran, or benzothiophene core were also suitable substrates for this reaction and led to the expected products 4k-m in up to 84% yield. It is also worth noting that the reaction is not limited to methyl ketones, and propiophenone or isobutylphenone derivatives 1n-o provided the corresponding secondary alkyl ethers 4n-o in high yields. In addition, cyclic silyl enol ethers 1p-1s were perfectly tolerated under these mild photoredox conditions and furnished the desired -methoxylated adducts 4p-4s in up to 83% yield. The reaction could also be efficiently applied to the synthesis of methoxylated -ketoester 4t in 67 % yield. Remarkably, the reaction was not limited to aromatic ketones, and 2-silyloxy dienes 1u-1v derived from aryl-or alkyl-substituted enones gave the desired adducts 4u-4v in good yields. This novel protocol was also fully compatible with the preparation of cyclic or acyclic aliphatic -methoxylated ketones 4w-4y.

To further demonstrate the broad applicability of this new -alkoxylation process, the late-stage alkoxylation of biorelevant molecules bearing various functional groups was explored. Thus, celestolide, phantolide, -ionone, and alanine derivatives 1z-1ac efficiently reacted to afford the desired -alkoxylated products 4z-4ac in up to 84% yield. Unlike 2-silyloxy dienes 1u-1v, 1-silyloxy dienes 1ad and 1ae generated from testosterone and (+)-nootkatone reacted exclusively in  position, without any trace of the methoxylated adduct: this selectivity was attributed to the allylic nature of the carbon-centered radical intermediate generated from the  addition. 11 Then, the scope of the reaction with respect to the alkoxy precursor was also explored. N-Alkoxypyridinium salts 2b-2h bearing a wide range of functionalities (nitrile, alkyne, phenoxy, ester, sulfone) reacted smoothly to provide the targeted alkoxylated ketones 4af-4al in good yields. The Scheme 2. Scope of the -alkoxylation reaction. a a General conditions: 1 (0.2 mmol), 2 (1.8 equiv), 3d (1 mol%) and K2HPO4 (2 equiv) in 2 mL of MeCN irradiated with blue LEDs at rt for 18 h. Given yields are for chromatographically pure isolated products. b Performed on a 10 mmol scale. c d.r. > 95:5. reaction efficiency was not dependent on steric hindrance, as bulkier alkoxy groups such as isopropoxy and 1,3difluoroisopropoxy were readily incorporated in up to 81% yield (4am-4an). Once again, it is worth to mention that no by-product resulting from -scission or HAT processes could be detected in the crude mixtures. [START_REF]Substrates highly prone to 1,5-HAT were not suitable for this -alkoxylation process[END_REF] To demonstrate the synthetic utility of this new photoredox-catalyzed protocol, a two-step, one-pot procedure was developed for the direct -alkoxylation of ketones. As revealed in Scheme 3a, treatment of ketone 5 with 2,4,6-collidine and triisopropylsilyl triflate provided first the corresponding silyl enol ether 1 without any purification or isolation. This intermediate was then exposed to -alkoxylation conditions to give the expected -alkoxylated ketones. This procedure was found to be applicable to the synthesis of ketones 4a,c,l,aj with good levels of reaction efficiency (57-72% yield).

We next sought to examine the applicability of this alkoxylation process to other carbonyl classes, notably silylketene N,O-acetals 6 derived from amides. To the best our knowledge, only one example of C-H -alkoxylation of amides has been recently reported, based on an "umpolung" strategy, which highlights the complexity of such transformation. [START_REF] Gonçalves | Unified Approach to the Chemoselective α-Functionalization of Amides with Heteroatom Nucleophiles[END_REF] To our delight, we found that the same reaction conditions could be applied to afford an array of functionalized -alkoxy amides 7a-g in good yields (Scheme 3b). Also, given the harsh oxidative conditions of umpolung strategies, aldehydes are quite challenging substrates for C-H alkoxylation. [START_REF] Lamhauge | Enantioselective α-Etherification of Branched Aldehydes via an Oxidative Umpolung Strategy[END_REF] We were therefore pleased to see that Nalkoxypyridinium 8 reacted smoothly under our mild photoredox conditions to furnish -alkoxylated aldehyde 9 in 52% yield (Scheme 3c). [START_REF]Our attempts to efficiently trap alkoxyl radicals intermolecularly by aldehyde-derived silyl enol ethers were all unsuccessful. It was attributed to the lower stability of radical intermediate B (secondary vs. tertiary)[END_REF] As the formation of -alkoxylated ketones 4 was inhibited in the presence of radical scavengers such as TEMPO ((2,2,6,6-tetramethylpiperidin-1-yl)oxyl) or in the absence of either photocatalyst or light (see Table 1), we presumed that a radical process was involved. To gain deeper insights into the reaction mechanism, a series of experiments were carried out. First, UV/vis experiments showed that an electron donor-acceptor complex was formed between 2a and K2HPO4 (see ESI for more details). 16 However, the involvement of this EDA complex in a productive reaction pathway is unlikely, as no desired product could be obtained in the absence of photocatalyst (Table 1, entry 15). Then, as outlined in Scheme 4a, luminescence quenching experiments showed that Ru-photocatalyst 3d was efficiently quenched by reagent 2a (kST = 2.5•10 -9 M -1 .s -1 ), compared to substrate 1a (kST = 6.2•10 -7 M -1 .s -1 ). This, combined with the values of the redox potentials of 2a and 3d (E2a = -0.47 V vs. oxidative quenching cycle, in which pyridinium 2a would be reduced by the excited state of Ru(bpy)3 2+ 3d (Scheme 4b). This SET reduction would generate the alkoxy radical A, which would be trapped by silyl enol ether 1 to lead to silyloxy radical B. Then, radical B is oxidized into carbocation C, which finally gives the desired -alkoxylated adduct 4 after desilylation. The oxidation step leading to intermediate C can follow two possible pathways: i) a photocatalytic pathway involving Ru III photocatalyst species; ii) or a radical chain mechanism involving a second molecule of Nalkoxypyridinium 2. To distinguish between these two pathways, the quantum yield of the reaction was measured (see ESI for more details). The maximum value  = 0.74 indicates that the radical chain process is not the major pathway in this transformation.

In summary, we have developed a new visible-light mediated process that allows for efficient intermolecular alkoxylation of various silyl enol ethers and silylketene N,Oacetals derived from a wide array of functionalized ketones and amides. Importantly, the synthetic value of this method was further demonstrated by the late-stage functionalization of more complex biorelevant molecules, and a one-pot protocol was readily applied to the direct -alkoxylation of carbonyl compounds. Extension to catalytic enantioselective -alkoxylation processes is currently under study and will be reported in due course.

ASSOCIATED CONTENT

Detailed experimental procedures, characterization data, and NMR spectra of new compounds. This material is available free of charge via the Internet at http://pubs.acs.org.

  Scheme 3. (a) One-pot -alkoxylation. (b-c) -Alkoxylation of other carbonyls. See ESI for more details.

Scheme 4 .

 4 Mechanistic investigations. a Luminescence quenching experiments. Stern Volmer plots I0/I = f([quencher]) using 1a or 2a as the quencher. b Plausible mechanism.

Table 1 . Optimization of the photoredox-catalyzed - methoxylation of silylated enol ether 1a. a

 1

	Entry Photocatalyst	x	Additive	Yield (%) b
	1	Ir(ppy)3 3a	1.2 -	27
	2	3a	1.2 Na2CO3	58
	3	3a	1.2 KOAc	17
	4	3a	1.2 Na2HPO4	<5
	5	3a	1.2 K2HPO4	61
	6	3a	1.2 K3PO4	45
	7	3a	1.2 Et3N	-
	8	3a	1.2 pyridine	15
		3a	1.2 2,6-lutidine 50
	10	3a	1.5 K2HPO4	75
	11	3a	1.8 K2HPO4	78
	12	Eosin Y 3b	1.8 K2HPO4	35
	13	Rose Bengal 3c	1.8 K2HPO4	40
	14	Ru(bpy)3(PF6)2 3d 1.8 K2HPO4	90 (89) c
	15	-	1.8 K2HPO4	-
	16 d	3d	1.8 K2HPO4	-

Funding Sources

This work was supported by a grant from the Agence Nationale de la Recherche with JCJC funding (ANR-20-CE07-0008-01).

AUTHOR INFORMATION

Corresponding Author * Guillaume Dagousset -Université Paris-Saclay, UVSQ, CNRS, UMR 8180 Institut Lavoisier de Versailles, 78035 Versailles Cedex (France); E-mail: guillaume.dagousset@uvsq.fr.


L Pizzagalli
email: laurent.pizzagalli@univ-poitiers.fr

Mechanical properties of C 60 at finite temperature from first-principles calculations

Keywords: mechanical properties, buckminsterfullerene C 60, first-principles molecular dynamics calculations

C 60 has been the focus of many studies since its discovery, but the mechanical properties of this unique molecule remain poorly known. To remedy that, we performed density functional theory molecular dynamics calculations of the compression of C 60 at finite temperature, within a large strain range. The possible recovery of the molecule after unloading was also investigated with first principles accuracy, for two different modes. We found that the behavior of C 60 follows three regimes depending on the compression strain. Up to 0.48 ± 0.02, the deformation is elastic with a small but distinct influence of the compression orientation. A relation between the electronic gap and the orientation is also revealed. The second regime starts for larger strains up to 0.75 ± 0.02, after a 30-32 nN force plateau. It is characterized by stochastic carbon bond breakings. Although bonds rupture is generally the fingerprint of plastic deformation in materials, we found that the C 60 structure can be dynamically recovered in a significant fraction of cases, following a fast unloading. This mode also yields defected structures with low energies compared to a slow unloading mode. Finally, the third regime corresponds to irreversible deformation of the molecule, for strains larger than 0.75 ± 0.02.

Introduction

The buckminsterfullerene C 60 belongs to the class of the most famous and studied molecules, partly due to its recognizable shape and its pioneering role in the emergence of other fullerenes like the carbon nanotubes. Since its identification in 1985, countless studies were conducted to discover and understand its characteristics. These efforts bore fruit as we have now at disposal a wealth of information for various of its properties in chemistry, electronics, optics, etc... [START_REF] Dresselhaus | Science of Fullerenes and Carbon Nanotubes[END_REF][START_REF]The Physics of Fullerene-Based and Fullerene-Related Materials[END_REF].

In contrast, the mechanical properties of C 60 remain poorly understood. The most prominent result concerns its extraordinary resilience i.e., the ability of C 60 to recover its initial shape after being severely compressed as a consequence of high energy impacts on surfaces [START_REF] Beck | Resilience of all-carbon molecules c 60 , c 70 , and c 84 : a surface-scattering time-of-flight investigation[END_REF][START_REF] Mowrey | Simulations of c 60 collisions with a hydrogen-terminated diamond (111) surface[END_REF][START_REF] Galli | Large scale quantum simulations: C 60 impacts on a semiconducting surface[END_REF][START_REF] Pan | Simulations of c 60 in collision with diamond surfaces[END_REF]. These pioneering studies were supplemented by a few theoretical and numerical investigations of the elastic properties of C 60 . The bulk modulus is determined to be in the range 710-860 GPa [START_REF] Woo | Bulk modulus of the c60 molecule via the tight binding method[END_REF][START_REF] Cai | Application of the tight-binding method to the elastic modulus of c60 and carbon nanotube[END_REF], surprisingly much higher than a predicted value of 2 GPa for the Young modulus [START_REF] Ghavanloo | Computational modeling of the effective young's modulus values of fullerene molecules: a combined molecular dynamics simulation and continuum shell model[END_REF]. A few researchers performed molecular dynamics calculations of the molecule compressed between graphene layers [START_REF] Brenner | Molecular dynamics simulations of the nanometer-scale mechanical properties of compressed buckminsterfullerene[END_REF][START_REF] Zhang | Simulations of the nanomechanical properties of compressed small fullerenes[END_REF]. Unfortunately these works do not allow for quantifying intrinsic mechanical properties of C 60 like the yield stress and the elastic limit, because of the presence of the deformed graphene layers. More recently, Shen reported the calculated variation of the compression load as a function of strain [START_REF] Shen | The compressive mechanical properties of cn(n=20, 60, 80, 180) and endohedral m@c60(m=na, al, fe) fullerene molecules[END_REF]. This work raises questions because surprisingly large maximum force and stiffness were predicted, which may be related to the compression method.

In brief, only conflicting and limited information is available concerning the mechanical properties of C 60 . This is in contrast to carbon nanotubes for which both elastic and plastic behaviors have been thoroughly investigated [START_REF] Yakobson | Mechanical properties of carbon nanotubes[END_REF][START_REF] Harik | Trends in Nanoscale Mechanics[END_REF]. Improving beyond the current state of knowledge is important for fundamental research, and would certainly be beneficial for a better understanding of the mechanical properties of aggregated or polymerized C 60based materials [START_REF] Hoen | Elastic properties of a van der waals solid:c60[END_REF][START_REF] Matsumuro | Mechanical properties of c 60 bulk materials synthesized at high pressure[END_REF][START_REF] Qiu | Physical and mechanical properties of c60under high pressures and high temperatures[END_REF][START_REF] Funamori | Large elastic deformation of c 60 nanowhiskers[END_REF].

Recently, a method allowing for the first principle modeling of nanoparticle compression was proposed, and its feasibility demonstrated in the case of C 60 [START_REF] Pizzagalli | Finite-temperature mechanical properties of nanostructures with first-principles accuracy[END_REF]. It is demonstrated in this work that the C 60 shell begins to yield for a compression strain close to 0.5, which contradicts the earliest studies using classical interatomic potentials where a much larger strain was predicted. However important aspects were left unexplored, such as the identification of the mechanical transition occurring at this strain or the influence of temperature. To fill the gap of knowledge we have performed a large number of simulations of the compression and unloading of a C 60 molecule in various conditions. They allowed us to determine with an unprecedented accuracy the mechanical properties of this molecule. In particular, three distincts regimes are identified depending on the compression strain. For low and high strain values the deformation is elastic or plastic, respectively. However, for medium strain values, we show that carbon bond breakings occurring during compression can be partially or even fully healed during unloading, especially by a dynamical recovery process. These findings are reported in the present paper. Calculations are first detailled in section 2, followed by a descrption of all results in section 3. The next two sections are a discussion of these results, and a conclusion. Compression and slow unloading simulations were performed following the method described in Ref. [START_REF] Pizzagalli | Finite-temperature mechanical properties of nanostructures with first-principles accuracy[END_REF]. For each run, two parallel planes were initially defined on either side and at equal distance of the C 60 molecule, along the z-axis (Fig. 1). The planes positions λ u and λ d were decreased and increased, respectively, at each iteration during a compression simulation (and the opposite during unloading). For all atoms whose z-coordinates were greater (lower) than λ u (λ d ), an additional compression force proportional to a constant K and to the square of the difference between the z-coordinate and λ u (λ d ), orientated along the z-axis and towards the center of the molecule, was applied (Fig. 1). This procedure allows for a gradual and smooth compression of the molecule, and has been largely used for classical molecular dynamics simulation [START_REF] Mordehai | Size effect in compression of single-crystal gold microparticles[END_REF][START_REF] Hale | Dislocation morphology and nucleation within compressed si nanospheres: A molecular dynamics study[END_REF][START_REF] Bel Haj | Influence of surface atomic structure on the mechanical response of aluminum nanospheres under compression[END_REF][START_REF] Kilymis | Uniaxial compression of silicon nanoparticles: An atom-istic study on the shape and size effects[END_REF][START_REF] Bian | Orientationdependent deformation mechanisms of bcc niobium nanoparticles[END_REF]. The originality of the present approach lies in the combination of these planar repulsion force fields with first principles molecular dynamics [START_REF] Pizzagalli | Finite-temperature mechanical properties of nanostructures with first-principles accuracy[END_REF].

Computational details

The electronic structure was calculated in the framework of density functional theory [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF][START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF], with a plane waves cutoff of 25 Ry for representing the wavefunctions. The Quantum Espresso suite [START_REF] Giannozzi | Advanced capabilities for materials modelling with quantum ESPRESSO[END_REF] was used. The exchangecorrelation contributions were computed using the PBE functional [START_REF] Perdew | Generalized gradient approximation made simple[END_REF]. Ultrasoft pseudopotentials [START_REF] Vanderbilt | Soft self-consistent pseudopotentials in a generalized eigenvalue formalism[END_REF] were employed to describe electron-ion interactions. The C 60 molecule was located in a (20 Å) 3 supercell, which is large enough to prevent artificial interactions between replicas [START_REF] Feuston | Electronic and vibrational properties of c 6 0 at finite temperature from ab initio molecular dynamics[END_REF]. A conjugate gradient energy minimization led to bond lengths equal to 1.401 Å and 1.451 Å, in excellent agreement with the literature [START_REF] Yannoni | Nmr determination of the bond lengths in c 60[END_REF].

The Car-Parrinello molecular dynamics scheme was used for compression and slow unloading simulations, with a timestep equal to 0.2 fs. The effective electron mass and the electron mass cutoff were set to 600 a.u. and 5 Ry, respectively [START_REF] Tassone | Acceleration schemes for ab initio molecular-dynamics simulations and electronic-structures calculations[END_REF]. Nose-Hoover thermostats were used for both electronic and ionic degrees of freedom, with an electronic kinetic energy target of 0.005 a.u. and an oscillation frequency of 90 THz for electrons, and a temperature target of 300 K and an oscillation frequency of 30 THz for the ions. Supplementary calculations at a ionic temperature of 1 K and an electronic kinetic energy target of 1.66 × 10 -5 a.u. were also carried out. The two planar repulsion force fields moved at a velocity of 0.05 Å/ps, yielding compression/unloading speeds of 0.1 Å/ps. The repulsion strength was set to 30 a.u. Finally, friction forces were applied on atoms already subject to compression forces, in order to prevent an excessive and steady rotation of C 60 during compression/unloading [START_REF] Pizzagalli | Finite-temperature mechanical properties of nanostructures with first-principles accuracy[END_REF]. Friction coefficients Γ equal to 50 and 500 a.u. were used in this work. The first value allows for slight and occasional rotations of C 60 with respect to the compression axis. On the contrary, using the second value nearly suppresses rotations. A typical compression run lasted 272000 iterations, which corresponds to a maximum deformation of about 80%. Two different C 60 orientations were investigated, with the compression axis passing either through an hexagon center (noted H in the following), or through a pentagon center (noted P).

Additional molecular dynamics calculations with no external forces were carried out with compressed molecules as initial configurations, to simulate a fast unloading. These calculations were performed using Born-Oppenheimer molecular dynamics with a timestep equal to 0.484 fs. In fact, the use of Car-Parrinello molecular dynamics was found to be inappropriate in that case because of the large temperatures generated by the fast release of mechanical energy. In these simulations, 1000 iterations with no thermal constraints were first carried out, followed by 1000 iterations with a Berendsen ionic thermostat set at 300 K.

Finally, nudged elastic band (NEB) calculations [START_REF] Jónsson | Nudged elastic band method for finding minimum energy paths of transitions[END_REF] were also performed to determine the activation energy associated with a carbon bond breaking under compression. These simulations included seven replicas, five being optimized using a quick-min algorithm. In order to preserve the strained state of C 60 in NEB calculations, topmost and bottommost carbon atoms were not allowed to move along the z-axis direction. These atoms were identified as being in contact to a planar repulsive force field during the compression stage. The NEB force convergence threshold was set to 0.05 eV Å -1 , and the spring constants to 0.04 hartrees. Two calculations were carried out in succession, using the climbing image method [START_REF] Henkelman | A climbing image nudged elastic band method for finding saddle points and minimum energy paths[END_REF] in the second one to determine the activation energy with high precision.

Several properties characterizing the atomic and electronic structures of the molecule were monitored during simulations. The potential energy is a direct outcome of DFT calculations. It is expressed relatively to the potential energy of a free C 60 molecule at 300 K (or at 1 K), which is obtained by an independent molecular dynamics calculation. These calculations were also used to determine the initial dimension d 0 of C 60 along the compression axis. The displacement during the compression/unloading is calculated as δ = d 0 -d, with d = λ u -λ d , and the engineering strain as ε = δ/d 0 . At last, the contact force during compression/unloading is given by summing the external compression forces due to one planar repulsion force field. Practically, it is calculated as the average of the two contact forces associated with upper and lower force fields. 1). Bonds between carbon atoms are drawn up to a separation of 1.81 Å. Grey bonds correspond to a local graphene-like atomic arrangement, as determined using the polyhedral template matching method and a root-mean-square deviation cutoff of 0.35 [START_REF] Mahler Larsen | Robust structural identification via polyhedral template matching[END_REF]. Unidentified local geometries with this method, typically due to broken bonds and local reconstructions, are shown with red bonds.

Results

A typical outcome of a simulation is presented in Fig. 2 (an extensive set of energy and force curves is shown in Fig. S1-S4 in Supplementary Material). During the first steps, the compression of the two caps in contact with the top and bottom repulsion force fields is observed (Fig. 2-(ab)). This C 60 flattening occurs with little lateral expansion as evidenced by the decrease of Poisson ratio up to a strain of about 0.3 ((Fig. 2). This initial stage is characterized by quadratic and linear increases of the potential energy and contact force, as expected for an elastic deformation. Note that non-linear contributions become clearly visible for strains greater than about 0.1, with the contact force seemingly reaching a plateau. At strains close to 0.3, a positive inflection can be seen for both energy and force curves. It is correlated with an increase of contact surfaces and the beginning of the lateral expansion of C 60 (Fig. 2-(b-c)), as evidenced by the augmentation of the Poisson ratio. The contact force saturates at strains greater than about 0.4, to reach a plateau defined by a value equal to 30-32 nN.

The contact force markedly decreases at strains greater than 0.47-0.49. In this stage, atomic bonds in the equatorial region of the compressed molecule are so severely stretched that one or two carbon bond break. This is usually followed by local reconstructions of the atomic structure, visible in the energy curve as sharp drops (Fig. 2-(d-e)). Such events occurred in all compression runs, for strains ranging between 0.50 and 0.70 approximately. Furthermore, for strains greater than 0.6, the slopes of both the energy and the contact force increase significantly. At large deformations the molecule resembles two flat carbon layers connected by bonds at their edges (Fig. 2-(e-f)). For strains equal to about 0.75-0.77, maximum energy and contact forces, close to 80 eV and 120 nN respectively, are reached. Ultimately the C 60 shell mechanically fails by the breaking of several bonds over a short strain interval (Fig. 2-(g)). A movie showing the whole sequence in relation with energy and force variations is available in the Supplemental Material. As mentioned above, both energy and contact force variations at low strains are in agreement with a linear elastic regime. The mechanical stiffness k can then be rigorously determined as the slope of the force-displacement curve (F = kδ), or as twice the quadratic coefficient of the energy-displacement curve (E = kδ 2 /2). Figure 3 shows the values calculated from all runs at 300 K, as a function of the displacement range used in the energy or force fit. For displacements lower than about 0.4 Å, accurate fitting is difficult because of the small size of the molecule and thermal vibrations, which leads to large and strongly decreasing stiffness values. The slopes of stiffness curves are significantly reduced but not completely canceled, by using larger displacements. This suggests that the elastic deformation is not perfectly linear. Note that similar trends and values are obtained from compression simulations carried out at 1 K (not shown).

Low strains

Figure 3 clearly indicates that the stiffness depends on the compression orientation. In fact, values ranging from 220 N m -1 to 180 N m -1 are computed for H, and from 185 N m -1 to 150 N m -1 for P. The difference is approximately 30 N m -1 for a large displacement range and from both energy and force fits. This compares extremely well with first-principles static calculations yielding a difference of 33 N m -1 between both orientations [START_REF] Wang | Mechanical and electronic properties of c 60 under structure distortion studied with density functional theory[END_REF]. Furthermore, Pawlak et al. previously reported that the C 60 stiffness depends on the position of the AFM tip apex used to press down on the molecule [START_REF] Pawlak | Atomicscale mechanical properties of orientated c 60 molecules revealed by noncontact atomic force microscopy[END_REF], in qualitative agreement with our work. Unfortunately, a direct quantitative comparison with their measurements is not possible since the loading processes are different.

Figure 3 also shows the Poisson ratio calculated for an increasing displacement range, in the elastic domain. As with stiffness, it is difficult to extract accurate values for displacements lower than 0.4 Å. Increasing the displacement range leads to a Poisson ratio close to 0.3 for both H and P i.e., with no clear dependence on the orientation. At 1 K however, different Poisson ratio values are calculated, with 0.31 for H and a value ranging from 0.21 to 0.24 for P. The reason for such a difference between 300 K and 1 K simulations is not straightforward. Note that a comparison with the literature is not possible since there are no published values to our knowledge. An interesting phenomenon is observed for runs with an initial H orientation and a low friction coefficient Γ = 50 a.u. For strains greater than about 0.2, small C 60 rotations are observed which change the C 60 orientation from H to P. This process typically takes place in two steps, with an intermediate state defined by the compression axis centered at the boundary between an hexagon and a pentagon. This reorientation process is in agreement with the lowest elastic energy and contact force for the P orientation (Fig. 3). It suggests that the reorientation with respect to the compression might be possible in the case of a physisorbed C 60 . By increasing the friction coefficient to 500 a.u., rotations are hindered and the initial H orientation can be preserved. Interestingly it is possible to detect the orientation change of C 60 by monitoring the variations of the HOMO-LUMO gap of the molecule during compression, shown in Fig. 4 for three different cases. The gap decreases for strains greater than 0.06, but while a large and seemingly linear reduction is obtained for the two H runs, only a weak decrease followed by a slow recovery is predicted for P. For the H run with a low friction coefficient, the HOMO-LUMO gap stops to decrease at a strain of about 0.2 and, after some fluctuations, sharply increases to reach the same value than for the P run. This behavior is obviously correlated with the reorientation of C 60 in this strain range. A maximum difference of almost 1 eV is predicted between both orientations. This would suggest that it might be possible to identify the compression axis with respect to the C 60 structure by measuring the HOMO-LUMO gap. Finally, for strains larger than 0.3, which corresponds to the beginning of C 60 lateral expansion, the HOMO-LUMO gap decreases for the P-oriented compression, whereas for the H case it first increases, saturates, then decreases towards zero. The C 60 gap closure under compression has already been reported in the literature [START_REF] Joachim | Electronic transparence of a single c 60 molecule[END_REF]. We further predict that the gap decrease rate and the threshold strains both depend on the compression orientation.

Large strains

Figure 5: Energy and contact force as a function of strain for different conditions: H orientation at 300 K (red lines) and at 1 K (black line), P orientation at 300 K (blue lines) and at 1 K (green line). Selected runs were all performed with a friction coefficient Γ = 500 a.u., in order to preserve the initial orientation during all the compression test. Three independent simulations were performed at 300 K to reveal the stochastic nature of the energy drops observed for strains greater than 0.48.

For strains greater than about 0.40, the contact force first experiences a plateau with a value of 30-32 nN, although the energy continues to increase (Fig. 2). This plateau is observed in all investigated cases i.e., independently of the orientation or the temperature, as can be seen in Fig. 5. It ends abruptly at a strain of 0.47-0.48 for the H orientation, and more smoothly for the P orientation at slightly higher strains of 0.49-0.50. Note that this approximate 0.02 strain difference between H and P cases is maybe not significant because of the aforementioned uncertainty in the definition of the zero strain. The presence of a similar force plateau was observed in indentation experiments on a C 60 molecules island [START_REF] Lantz | Measurement of forces during the modification of c60 islands[END_REF], albeit at a lower value of about 20 nN. However, our results disagree with previous theoretical studies. In fact, Shen reports a similar plateau but at a surprisingly high value of about 1100 nN from approximate quantum mechanical calculations [START_REF] Shen | The compressive mechanical properties of cn(n=20, 60, 80, 180) and endohedral m@c60(m=na, al, fe) fullerene molecules[END_REF]. Furthermore, classical interatomic potential simulations lead to a contact force maximum of 20 nN for similar strains, except that no plateau is obtained [START_REF] Seo | Molecular dynamics investigation on the nano-mechanical behaviour of c60 fullerene and its crystallized structure[END_REF]. The contact force decreases for strains greater than 0.47-0.50, due to the breaking of one or more C-C bonds in the equatorial region of the C 60 shell. This is demonstrated in Fig. 6 which shows the variation of the maximum C-C bond length for several cases. For a P orientation, one can observe a first inflection at a strain close to 0.3 corresponding to the beginning of the lateral expansion of C 60 (Fig. 2). A second inflection point occurs at a strains 0.48-0.50, the almost vertical increase revealing the breaking of a bond. The structures analysis shows that the C 60 shell exhibits an armchair configuration for a P orientation (Fig. 6). C-C bonds that break are perpendicular to the compression axis, allowing for a maximum stretching, and form the edges between hexagons and pentagons. For systems with an H orientation, inflection points occur first at about 0.38, and then at 0.48 (Fig. 6). They are characterized by a zigzag configuration along the compression axis (Fig. 6), and the bonds breaking first form the edges between two hexagons. Inflection points occur first at about 0.38, and next at 0.48 (Fig. 6).

Note that it is always the same kind of C-C bond which breaks first for either H or P orientations. However this initial event triggers diverse structural modifications of the C 60 molecule, often involving partial reconstruction in the vicinity of the broken bond. These reconstructions events correspond to well defined drops on energy curves, occurring immediately after the first bond breaking or at somewhat larger strains (Fig. 5). Reconstruction mechanisms involving the most atoms are typically associated with the largest energy drops. On the contrary, some light structural changes are unnoticeable from energy variations only, because the energy release is more than compensated by the elastic energy increase during the ongoing compression. It is the case for the first bond breaking event in P configurations for instance.

Various C 60 configurations are then obtained, with one or two shell openings which could be weakly or largely reconstructed (see fig. 2-(d-e), and Figs. S1-S4 in Supple-mentary Material). No clear correlation is found between the structural change or its associated energy drop, and the strain at which it occurs. Furthermore, reproducibility is only achieved for 1 K simulations. All these factors clearly demonstrate the stochastic character of the reconstruction mechanisms following the initial bond breaking event.

The last salient characteristics of the compression curves are the high force and energy maxima reached at large strains, followed by significant drops due to the total flattening of the C 60 shell. These maxima occur in a strain range 0.74-0.78 for the energy, and 0.73-0.76 for the contact force, depending on cases. On averages, the maximum energy and force before the drops are 83 eV and 126 nN, with extreme values as high as 96 eV and 150 nN, respectively. However, these numbers should be considered with caution, especially for the force, since they depend on the averaging procedure used to smooth fluctuations. In this section, the results of unloading simulations are described. Several configurations are first selected in various runs and at different compression strains ε i . For each of them, a molecular dynamics simulation is performed using the same conditions than for the compression, but now moving off the two planar repulsive planes. These simulations are referred to as slow (S) unloading in the following. Figure 7 shows an example of computed energy and contact force during a loading-unloading cycle. Compression and unloading curves are similar at first and begin to diverge at an approximate strain of 0.625, which corresponds to a bond breaking and reconstruction event. Energy and contact force curves in the unloading regime are typical of an elastic behavior, except for the sharp decrease (increase) in energy (force) at a strain of 0.506, suggesting structural reconstructions. We define the residual plastic strain ε p as the strain at which the contact force reaches zero during the unloading stage (alternatively the potential energy becomes constant). Another interesting quantity is the final energy E D at the end of the unloading stage, called defect energy in the following. Both ε p and E D are equal to zero when a perfect C 60 molecule is recovered after unloading. We tested a different unloading mode (F) where the C 60 molecule is relaxed with no planar repulsive fields, first with no temperature control then with a 300 K thermostat. The energy-strain relation is represented in Fig. 7. The two temperature regimes are clearly visible, with large strain fluctuations and a high energy in the first one, followed by a significant potential energy reduction in the second one. E D is obtained by averaging the energy in the second regime. S unloading simulations were carried out for initial configurations corresponding to strains ε i ranging from 0.50 to 0.74, at which one or several bond breaking mechanisms occurred. Calculated defect energies E D and plastic strains ε p are reported in Fig 8 . Overall, high E D and ε p values appear to be correlated with large initial strains. Our data also suggest an approximate linear relation between E D and ε p (Fig. 8). A first domain can be roughly defined for ε i 0.65, which corresponds to initial structures where the C 60 shell is opened at one place in most cases, with a moderate reconstruction (see Fig. S1-S4 in supplementary materials for examples). In that case S unloading calculations yield final configurations with E D <15 eV and ε p <0.1. The elastic energy is largely relaxed, but structural defects are not or imperfectly cured during unloading. Those are shell openings or local modifications of the normal pentagon-hexagon patterning (Fig. 9-(a)). Only 4 among 33 simulations allow for retrieving the perfect C 60 structure (E D = 0, ε p = 0). In these runs the initial structures (shown in Fig. S1-(a,c) and Fig. S2-(d)) include two or three broken bonds with no local reconstructions, which can be reformed straightforwardly during the unloading. It suggests that remaining defects after unloading are a consequence of the bond reconstruction events occurring during compression.

Unloading

For initial strains ε i 0.65, we find that E D ranges from 10 eV to at most 28 eV whereas ε p ranges from 0.03 to 0.35 (Fig 8). The pristine C 60 molecule is then never recovered, and a significant part of the mechanically stored energy is not relaxed during unloading (see two examples in Fig. 9-(b,c)). In fact the initial structures are heavily distorted with several reconstructed bonds, which can hardly be recovered during unloading.

Figure 10: Top: potential (black) and kinetic (green) energies as a function of the MD time (ps) during a F unloading simulation (1000 iterations i.e., 0.484 ps in a NVE ensemble, followed by 1000 iterations with a 300 K ionic thermostat). The potential energy reference is equal to -10103.003 eV, which is the average potential energy of a perfect C 60 at 300 K in our simulations. The initial configuration corresponds to a compression strain of 0.67. Middle: C 60 height (red) and height velocity (blue), for the same simulation. The C 60 height is determined at each iteration as the difference between the maximum atomic separation along the axis used for compression. The height velocity is calculated as the derivative of the C 60 height with respect to the MD time. Bottom: Selected snapshots from the MD simulation at selected times in the range [0 -0.116] ps (yellow frame).

We now focus on the results of F unloading simulations. An example of the variations of quantities monitored during a Born-Oppenheimer molecular dynamics calculation is shown in Fig. 10, together with snapshots of the C 60 structure at selected times. At startup, the potential energy of the deformed C 60 molecule is equal to 51.7 eV, and the kinetic energy to 3N k B T /2 = 2.33 eV with T = 300 K. A fast expansion with an increasing speed is observed during the following steps, similar to a spring release along the compression direction. The C 60 shell extremities reaches a maximal speed of 112 Å ps -1 at 0.046 ps, and a height equal to 7.3 Å at 0.108 ps. The release process appears to be already completed at 0.11 ps with both kinetic and potential energies oscillating between an average value of 27 eV, in agreement with the equipartition theorem. Although the vibrational mode Hg(1) corresponding to the longitudinal vibration along the compression axis [START_REF] Menéndez | Vibrational spectroscopy of c[END_REF] remains dominant, other phonon modes become populated during the following iterations. Assuming all kinetic energy could be equated to a thermal energy, an average temperature of 3480 K is determined (Fig. 10). After 1000 iterations i.e. 0.484 ps, a 300 K thermostat is applied and the potential energy becomes negligible, suggesting that a perfect C 60 molecule is recovered. A movie showing the C 60 evolution is available in the Supplementary Material. A structural analysis indeed confirms that all initially present topological defects disappeared in this simulation, and that healing occurred in a very short timespan between 0.077 ps and 0.11 ps (Fig. 8). A reconstruction mechanism is illustrated in Fig. 11 for the same system. The initial configuration includes an heptagon encircled by hexagons at 0.084 ps (Fig. 11-(a)). Two new bonds are concurrently created while one is broken, leading to the formation of three adjacent pentagons at 0.09 ps (Fig. 11-(b)). The central pentagon is next transformed into an hexagon by a bond-switching mechanism between 0.098 ps and 0.103 ps (Fig. 11-(c-d)). The regular alternating hexagonpentagon pattern is finally recovered (Fig. 11-(c-d)).

A perfect C 60 geometry (i.e. E D = 0 eV) is recovered in one third of the 45 performed F unloading simulations i.e., in a significantly higher proportion than in S unloading calculations. It occurred more frequently for initial configurations with strains ranging from 0.50 to 0.65, but also for three cases with larger strain values (Fig. 8). Another interesting result is the overall lower defect energies for F unloading compared to S unloading simulations (Fig. 8). In fact, for F unloading, E D is always lower than 17.2 eV (72% of the S simulations). This can be explained by less topological or coordination defects left after F unloading. Furthermore, 62% of the runs lead to E D lower than 5 eV, to be compared to only 18% for S unloading. Such small E D values are typically associated with only one topological defect, as seen in Fig. 9-(a). Our calculations reveal that the compression of C 60 can be analyzed according to three strain domains. For strain values lower than 0.47-0.50, the deformation is purely elastic i.e., a perfect molecule is recovered when the load is removed. For larger strains, but still lower than 0.73-0.76, bond breaking events occur in the equatorial region of the molecule. However partial or total structure recovery can be obtained, with a significant influence of unloading conditions. Finally, for strains greater than 0.73-0.76, the structural damages become severe enough to prevent any recovery. We favor the use of contact force rather than energy to determine these strain thresholds. In fact, a thorough analysis of many simulations reveals that the energy drops visible in Figs. 2 or 5 are mostly due to stochastic bond reconstruction events, occurring at higher strains than the first bond breaking. Therefore the contact force is a better quantity than energy to determine the onset of C 60 shell failure.

Discussion

With this criterion, the first strain threshold ranges from 0.47 to 0.50, taking into account all calculations. There is an intrinsic uncertainty associated with the determination of the strain for such a small system. Furthermore we find that critical strain values depend on the compression orientation, and also on temperature (Fig. 6). This thermal dependence indicates that the first bond breaking event is also a stochastic event. A critical issue is therefore the possible influence on our results of the high strain rate used in molecular dynamics simulations. Theoretical analysis indeed show that a lowering of the computed elas-tic limit and of the critical stress is expected when extrapolating to typical experimental strain rates [START_REF] Li | The mechanics and physics of defect nucleation[END_REF][START_REF] Zhu | Temperature and strain-rate dependence of surface dislocation nucleation[END_REF]. To tackle this issue we performed NEB calculations of the activation energy for the first bond breaking event in the case of a H orientated molecule. An initial reaction path was built between a perfect C 60 structure and another where a C-C bond was broken (Fig. 12-(a-b)), both configurations being relaxed and for several compression strains. The results of the NEB calculations are represented in Fig. 12-(c). As expected, the activation energy ∆E is low and strongly dependent on the compression strain. In fact, it is equal to 0.218 eV for a strain of 0.459, and extrapolation of the data shows that the barrier vanishes for a strain close to 0.476. Our calculations also reveal that the bond broken C 60 structure (Fig. 12-(b)) becomes unstable at a compression strain of 0.453, the C-C bond reforming during relaxation. This defines an absolute minimum value for the critical strain at which the first bond breaking event occurs. The thermal activation of this mechanism is therefore only possible in the strain range 0.453-0.476. This suggests that the influence of the strain rate used in MD calculations is very limited, with an overestimation of the elastic limit of at most 2%.

Quantitative information on the onset of C 60 shell failure is scarce in literature. Available studies rely on empirical potentials description and do not provide explicit critical strains [START_REF] Brenner | Molecular dynamics simulations of the nanometer-scale mechanical properties of compressed buckminsterfullerene[END_REF][START_REF] Zhang | Simulations of the nanomechanical properties of compressed small fullerenes[END_REF][START_REF] Seo | Molecular dynamics investigation on the nano-mechanical behaviour of c60 fullerene and its crystallized structure[END_REF]. Our own investigations revealed that different potentials all overestimate the elastic limit, with values greater than 0.60 [START_REF] Pizzagalli | Finite-temperature mechanical properties of nanostructures with first-principles accuracy[END_REF]. However, a recent work by Seo and Kim suggests that the first bond breaking event occurs for a C 60 height of about 4.5 Å i.e., a critical strain of 0.35 [START_REF] Seo | Molecular dynamics investigation on the nano-mechanical behaviour of c60 fullerene and its crystallized structure[END_REF]. The cause for such a disagreement with the present study and published works is unclear.

The second important result to discuss is the partial or full recovery of the defected C 60 , with an appreciable influence of the unloading conditions. This regime is defined by compression strains between 0.48±0.02 and 0.75±0.02, where the C 60 energy is in the range 40-80 eV. The various bond breaking events occurring in this regime are associated with energy drops of 4 eV on average (with extreme values of 1 eV and 9 eV), computed from all our simulations. This value has to be compared to the maximum potential energy reached during compression, that is 82 eV on average. Only a small fraction of this energy is released by the first bond breaking events, and the remaining part is stored as an elastic energy in the molecule. In F unloading simulations, this energy is quickly converted into kinetic energy and the fast uniaxial expansion of the molecule helps to heal the defects created during the compression, leading to full recovery in several cases. In contrast, in S unloading simulations the stored elastic energy is essentially converted to thermal energy, and lost by exchange with the thermostat. Full recovery in these conditions only happens when C 60 shell openings are not reconstructed.

S simulations can be compared to the unloading stage in a nanoindentation experiment. Usual strain rates of 10 -1 -10 -3 s -1 imply that the sample is always at thermal equilibrium. This is reproduced in the simulations by the use of a ionic thermostat. Conversely F unloading simulations are comparable to experiments where fullerenes with high kinetic energy are colliding with materials surfaces [START_REF] Beck | Resilience of all-carbon molecules c 60 , c 70 , and c 84 : a surface-scattering time-of-flight investigation[END_REF]. The C 60 molecule is observed to rebound with no fragmentation for a defined energy range. Pionering simulations showed that the C 60 molecule is first flattened, then leaves the surface with a certain speed [START_REF] Mowrey | Simulations of c 60 collisions with a hydrogen-terminated diamond (111) surface[END_REF][START_REF] Galli | Large scale quantum simulations: C 60 impacts on a semiconducting surface[END_REF]. The release process is equivalent to our F unloading simulations. However only a qualitative comparison can be made, since in experiments a non negligible part of the incoming molecule energy is dissipated by interacting with the substrate. Also, the elastic discharge of the latter is likely to impact the release velocity. Nevertheless, our results suggest that most rebounding molecules are not perfect and include structural defects, in agreement with the early study by Galli and Mauri [START_REF] Galli | Large scale quantum simulations: C 60 impacts on a semiconducting surface[END_REF]. Furthermore, they also mention that a significant proportion of these molecules could be free of any defects.

Finally, it is worth emphasizing that C 60 was described as a resilient molecule in these pioneering experiments, due to its ability to 'survive' these high energy collisions. In materials science resilience is formally defined as the ability of a material to withstand elastic deformation without deforming plastically, and it is usually expressed as an energy per volume. For instance, the resiliences of steel and carbon fibers are approximately 6 MJ m -3 and 25 MJ m -3 , while for spider silk it can be as large as 180 MJ m -3 [START_REF] Römer | The elaborate structure of spider silk[END_REF]. Our study shows that an elastic behavior is always obtained for compression strains lower than 0.48 ± 0.02, with a maximum elastic energy of 40 eV on average. Assuming that the diameter of C 60 is approximately 7 Å, this energy corresponds to an extraordinary resilience of 36 × 10 3 MJ m -3 i.e., 200 times higher than for spider silk.

Conclusions

Density functional theory molecular dynamics calculations have been carried out to determine the mechanical properties of the buckminsterfullerene C 60 at finite temperature and with first-principles accuracy. In particular, we focus on the compression of the molecule as well as on two different unloading modes. These calculations lead to several original findings, which further illustrate the unique features of C 60 . Hence three different domains can be defined depending on the compression strain. In the first one, for strain values below 0.48 ± 0.02, the deformation is exclusively elastic. It is shown that the stiffness and the electronic gap depend on the molecule orientation with respect to the compression. Second, in the domain defined by strains between 0.48 ± 0.02 and 0.75 ± 0.02, C-C bond breaking events occur in the equatorial region with a clear stochastic character. This onset of C 60 shell rupture is associated with a contact force plateau of 30-32 nN regardless of temperature or compression orientation. A surprising finding is the possible recovery of the C 60 shell structural integrity after unloading. This is especially true for fast unloading which yields a significant proportion of fully recovered molecules, or with low energy defects, compared to a slow unloading. This dynamical recovery phenomenon might be explained by the high stored elastic energy compared to the energy associated with C 60 shell breaking events. Finally, for strains greater than 0.75 ± 0.02, the C 60 shell damaging reaches a point of no return, with no possible recovery by unloading. The existence of three strain regimes is unusual, and further emphasizes the uniqueness of fullerenes compared to other nanostructures [START_REF] Amodeo | Modeling the mechanical properties of nanoparticles: a review[END_REF].

Finally, our study is based on first principles calculations which are more accurate and reliable than empirical potentials. This is especially true in situations of large deformations as encountered here. Another advantage of first principles calculations is that including different chemical species is straightforward. A natural extension of the present work is therefore to investigate the influence on mechanical properties of endohedral inclusions like noble gas atoms or small molecules in fullerenes. Works are currently in progress in these directions. One could also mention the potential application of the present method to ordered C 60 thin films. A recent study suggests that compression simulations lead to similar results for a single molecule and a thin film, at low and moderate strains [START_REF] Seo | Molecular dynamics investigation on the nano-mechanical behaviour of c60 fullerene and its crystallized structure[END_REF]. It would then be interesting to carry out first-principles MD calculations for further confirmation.

Figure 1 :

 1 Figure 1: Compression of the C 60 molecule between two planar repulsion force fields (represented by the thin turquoise slabs), located at z-coordinates λu and λ d , and moving towards each other.

Figure 2 :

 2 Figure 2: Energy (red line), contact force (blue line), and Poisson ratio (black line) as a function of compression strain for a P orientation (Γ = 500 a.u.). Curves are smoothed by performing rolling averages with a strain period of 0.00294 for the energy, of 0.0294 for the contact force, and 0.0588 for the Poisson coefficient. (a)-(g) pictures showing the configurations obtained at different strains: 0 (a), 0.20 (b), 0.45 (c), 0.55 (d), 0.60 (e), 0.735 (f), 0.795 (g). The turquoise thin slab symbolizes the bottom repulsive force field (The upper one was removed for a clearer picture compared to Fig.1). Bonds between carbon atoms are drawn up to a separation of 1.81 Å. Grey bonds correspond to a local graphene-like atomic arrangement, as determined using the polyhedral template matching method and a root-mean-square deviation cutoff of 0.35[START_REF] Mahler Larsen | Robust structural identification via polyhedral template matching[END_REF]. Unidentified local geometries with this method, typically due to broken bonds and local reconstructions, are shown with red bonds.

Figure 3 :

 3 Figure 3: Stiffness and Poisson ratio of C 60 as a function of displacement, during compression at 300 K. The yellow and green colors correspond to simulations with H and P orientations, respectively. The widths of the plots display the variations associated with several runs. Top: stiffness determined by an affine fit of the contact force (represented in the inset plot). Middle: stiffness determined by a quadratic fit of the C 60 potential energy (shown in the inset plot). Bottom: Poisson ratio.

Figure 4 :

 4 Figure 4: HOMO-LUMO gap of C 60 as a function of compression strain, for three simulations: strong friction (Γ = 500 a.u.) with either H (red curve) or P (blue curve) orientations, and weak friction (Γ = 50 a.u.) with H orientation (green curve).

Figure 6 :

 6 Figure 6: Maximum C-C bond length as a function of strain for different conditions: H orientation at 300 K (red line) and at 1 K (black line), P orientation at 300 K (blue line) and at 1 K (green line). The top pictures show the H and P configuration at a 0.47 compression strain, with 'equatorial' bonds highlighted in violet. Solid lines indicate bonds which are observed to break during compression.

Figure 7 :

 7 Figure 7: Energy (red lines) and contact force (blue lines) as a function of strain during the compression/unloading of a C 60 molecule (P, Γ = 50 a.u., 1 K). The zero energy reference corresponds to the initial perfect C 60 molecule. Dashed lines correspond to the compression stage, whereas thick full lines correspond to unloading obtained by slowly moving off the two planar repulsion force fields (S unloading), starting at a strain ε i = 0.668. The residual plastic strain εp is determined when the contact force reaches zero during unloading. Also shown in the figure is the energy variation (thin red line) resulting from the unconstrained evolution of the molecule (F unloading). For both S and F unloading modes, we define the defect energy E D as the energy difference between the final configuration and the perfect C 60 .

Figure 8 :

 8 Figure 8: Left: Defect energy E D as a function of the initial compression strain ε i for several configurations. The red (blue) dots correspond to S (F) unloading simulations, respectively. Right: Defect energy E D versus the final strain εp in S unloading simulations (green dots). The green line is the result of a linear regression fit.

Figure 9 :

 9 Figure 9: Final geometries obtained after three independent S unloading calculations. The initial configurations are obtained from a compression simulation with a P orientation (Γ = 50 a.u., 300 K), at strains ε i equal to 0.623 (a), 0.668 (b), and 0.728 (c). Bonds between carbon atoms are drawn for a maximum separation of 1.81 Å. Structural defects are marked with violet bonds.

Figure 11 :

 11 Figure 11: (a-e) Example of a reconstruction mechanism of the C 60 shell occurring during a F release simulation. Only a part of the shell is shown for clarity. The colored digits indicate how many C atoms belong to each ring (blue for heptagons, green for hexagons, orange for pentagons). The red arrows show how new bonds will form, whereas dashed black lines show where bonds disappeared.

Figure 12 :

 12 Figure 12: Initial (a) and final (b) atomic configurations (top view) used in nudged elastic band calculations. Carbon atoms with coordination 3 (2) are shown as black (red) spheres. (c) Calculated variation of the total energy (relative to the energy of the initial configuration) for the bond breaking reaction path (from (a) to (b)), for different compression strains: 0.459 (red), 0.465 (blue), 0.471 (green). The corresponding energy ∆E as a function of strain for the transition (a)-→(b) is reported in the inset graph. The extrapolation of this curve suggests that the energy barrier vanishes for ε 0.476.

Acknowledgements

Computer time for this study was provided by several computing facilities: the Spin Center at the University of Poitiers, the MCIA (Mésocentre de Calcul Intensif Aquitain), and GENCI-CINES (Grant 2020-A0090912035). We also acknowledge support by the French government program "Investissements d'Avenir" (EUR INTREE, reference ANR-18-EURE-0010).

Declaration of competing interest

The author declares that he has no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


P Dufour
email: dufour@lagep.univ-lyon1.fr

D Edouard

H Hammouri

Keywords: Destruction of COV, reverse flow reactor, optimal control

This paper aims to compare 2 optimal control laws for a catalytic reverse flow reactor (RFR), which topic has not often been tackled. The RFR aims to destruct volatile organic compounds (VOCs), which are noxious products and are the source of a lot of problems: pollution, acid rains, woods wasting, green house effect and health hazards. Due to the public regulations, the VOCs discharge in the atmosphere becomes therefore more and more limited. The RFR is characterized by a periodic cycle of treatment of the VOC flow, allowing trapping the heat of reaction inside the core of the RFR. This allows using the RFR in an autothermal mode. The control problem tackled here is to maintain the RFR is this mode, such that the RFR extinction is avoided, and such that the catalyst elements are not destroyed. Effects of the influencing stochastic variations of the inlet pollutant concentration (the input disturbance) have to be accounted for by the optimal constrained tuning of the dilution rate and the internal electric heating. Optimal closed-loop performances of a linear quadratic regulator (LQR) and a model predictive controller (MPC) are compared through simulations.

Destruction of volatile organic compounds by the optimal multivariable control of a catalytic reverse flow reactor INTRODUCTION

The catalytic RFR aims to destruct volatile organic compounds. These VOCs are noxious products and they are the source of a lot of environmental problems: pollution, acid rains, woods wasting, greenhouse effect and health hazards. The catalytic reverse flow reactor used here allows high temperatures in catalyst bed whereas the inlet and outlet gas stream temperatures are close to ambient temperature. Indeed, through periodic flow reversal, the heat released by the reaction is first trapped in the packing and is then used to heat up the feed when the flow direction is reversed.

Because of the high efficiency of heat exchange between gas and solid phases, autothermal operation is possible even for a feed with a low adiabatic temperature rise (below 15 K). These features make therefore RFR highly competitive for VOCs combustion.

The RFR is modeled by a nonlinear dynamic partial differential equation (PDE) system characterized by complex nonlinearities in the spatial domain [START_REF] Hammouri | Observer design for reverse flow reactor[END_REF]. Explicitly, transport reaction phenomena with significant diffusive and convective phenomena are typically characterized by severe nonlinearities and spatial variations, and are naturally described by partial differential equations. Examples of such processes include tubular reactors, packed bed reactors, absorption columns, drying or curing processes.

In control theory, due to the complexity of the problem, relatively few studies are devoted to the control of processes explicitly characterized by a PDE model, especially in the nonlinear case and even with only one spatial dimension. Even if various methods are proposed to control such distributed parameter systems, there is no general framework yet. In order to implement, with a computer, a low order model based controller, the original PDE model is usually simplified into an ordinary differential equation (ODE) model. Such a finite dimensional approximation is based on the finite differences method, the finite volume method, the orthogonal collocation method or the Galerkin's method. Other works utilized properties of the initial PDE system before finite dimension controller synthesis. In a previous work, the single input single output (SISO) control of the RFR was treated (Dufour et al., 2003a). The aim was to control the RFR such that the outlet gas concentration released in the atmosphere was maintained below a maximum level fixed by public regulations. This control strategy was based on a parabolic PDE model, an internal model control (IMC) structure and a MPC framework (Dufour et al., 2003b). According to various regimes, it was shown that the proposed controller was able to tune correctly the control action, i.e. the heating power at the core of the reactor. In the meantime, for another particular regime (rich feed), the temperature inside the reactor was exceeding a threshold temperature that reflects the deterioration of the catalytic elements. Moreover, the input disturbance was assumed to be relatively constant and measured, whereas unmeasured and large stochastic variations need to be accounted for in reality. The aim of this paper is to provide a multivariable control framework to solve all of these problems. This requires modifying the PDE model and the control problem to account for a new second manipulated variable: the cooling action.

In this paper, two model-based control strategies for the constrained optimal control of the RFR are compared: LQR and MPC which both use a high-gain observer. The LQR is used here in a classical framework [START_REF] Levine | The control handbook[END_REF]. MPC was developed for ordinary differential equation models and is well dedicated to solve a constrained problem [START_REF] Morari | Model predictive control: past, present and future[END_REF].The IMC-MPC approach (Dufour et al., 2003b) is used here. The paper is structured as follows: in Section 2, the RFR and its PDE model are briefly presented. Section 3 aims to remind the observer previously designed. Section 4 deals with the LQR and MPC strategies used here. Simulation results given allow comparing the performances of the LQR and the MPC for this process.

PROCESS DESCRIPTION

A medium-scale RFR (Nieken et al.) is considered (Figure 1). Cordierite monoliths of square cross sections with channels of 1*1 mm are packed in the reactor. Monolith in the core region is catalytically active and is inert in both end sections. A blower located downstream of the RFR keeps aspiration of the pollutant at a constant flow rate. In the core region, an electric heater maintains ignition temperature, while the temperature in the catalytic layer is decreased by fresh air dilution. The packed layer is adiabatic, except in the core region where heat loss is inevitable due to both the installation for air dilution and the high temperature in this region. High temperatures exist in catalyst bed whereas the inlet and outlet gas stream have ambient temperature. Indeed, through periodic flow reversal, heat released by reaction is first trapped in the packing and then used to heat up the feed.

The model considered here for control purpose is obtained from a counter current pseudo-homogeneous model [START_REF] Hammouri | Observer design for reverse flow reactor[END_REF], accounts for mass transfer limitation and periodic frequency correction: it features one nonlinear parabolic PDE, two algebraic equations, and nonlinear boundary conditions. The nonlinearities are due to the cooling action. The advantages of this model are that it is more accurate and faster to compute than a previous model used for control (Dufour et al., 2003a). This countercurrent reactor model is a heterogeneous model. In order to homogenize and simplify this model, the following model has been introduced by [START_REF] Hammouri | Observer design for reverse flow reactor[END_REF]. It is based on the method described in [START_REF] Balakotaiah | Effective models for packed bed catalytic reactors[END_REF] and it assumes that the kinetic reaction can be neglected under strong mass transfer limitation. It allows obtaining the following pseudo-homogeneous model: it features one PDE, two algebraic equations to account for mass transfer limitation, and a periodic frequency correction. Normalizing some variables, we obtain:

Heater T 1 T 2 H/2 x = 0 x = 1 x = ξ
2 2 2 1 1 2 2 ( , ) ( , ) ( , ) 1 1 1 ( ) ( ) ( ) ( , ) 2 2 ( , ) ( , ) 0 ( , ) ( , ) 0 s s s ad ax T x t T x t T x t P T t x t t P P x x x t P x t x x t P x t x θ θ θ θ θ α α τ ψ ω α ω ω ω  ∂ ∂ ∂ + - = + + + ∆  ∂ ∂ ∂   ∂  + =  ∂   ∂ - + =  ∂   (1) Where / 2 z x H =
is the normalized space variable, t is the time variable), s T is the temperature inside the solid parts, 1 w and 2 w are the VOC mass fraction in the upstream and downstream monoliths. Some model parameters are also parameterized by the period of flow reversal θ (assumed constant):

0 0 0 0 0 0 (1 ) 2 2 2 1 ( )(1 ) 2 (1 ) 2 c v pmg s ax ax v pmg s ps v pmg ha H P u c P H u c c H u c θ θ θ ρ τ λ θ ρ τ ε ρ τ ρ  = -     =  -   -  =   (2)
And the non constant coefficients depend on the VOC mass fraction 10 ( ) t ω fed into the RFR:

10 1 2 10 ( ) ( ) ( , ) ( , ) ( , ) ( ) 2 ( ) ad pmg H T t t Mc x t x t x t x t ω ω ω ψ ϕ ω ∆  ∆ =    +  =  
(3) where ( )

x ϕ
accounts for the type of monoliths: ( )

x ϕ =0 in the inert monoliths ( x θ ξ < ) and ( ) x ϕ =1 in the catalytic monoliths ( x θ ξ ≥ ).
The boundary conditions for (1) are at x=0:

1 10 1 0 ( , ) ( ) ( , ) ( , ) ( , ) s g s x t t T x t T x t T x t T P x θ ω ω α =   ∂  = - =  ∂  (4)
and at x=1:

2 0 1 0 0 0 1 2 1 2 (1 ')( ( , ) ) ( ( , ) ) ( , ) ( , ) ( , ) ( , ) 1 ( , ) ( , )
( , ) ( , )

j g g v pmg s g s s g s Q N T x t T T x t T S u C T x t T x t T x t P x T x t T x t T x t P x x t x t θ θ α ρ α αω ω  + - = - +    ∂  = -  ∂   ∂  = +  ∂  =   (5)
The initial conditions for (1) are at t=0:

1 2 0 ( , ) ( , ) ( , ) g g s T x t T x t T x t T = = = (6)
In the above equations, T 0 is the ambient temperature and the feed temperature and T g1 and T g2 are upstream and downstream gas temperatures respectively. More details about the model parameter may be found in [START_REF] Edouard | Observer based multivariable control of a catalytic reverse flow reactor: comparison between LQR and MPC approaches[END_REF]. The cooling manipulated variable α is accounted for at any x, and introduces a strongly non-linear behavior. The dilution rate (1-α ) is the percentage of fresh air in downstream flow. The heating power manipulated variable Q j (in watts) is accounted for in the boundary condition at x=1.

OPTIMAL CONTROL PROBLEM

According to the operating conditions, various problematic behaviors can take place during the operation of the RFR:

• When the feed of pollutant concentration 10 ( )

t ω (or ( ) ad T t ∆
as seen in ( 3)) is too rich, the release of heat due to the reaction produces thermal overheating that deteriorates the catalysts. The temperature inside the reactor has therefore to be maintained under the maximum temperature specified as 600 K [START_REF] Ramdani | Reverse flow reactor at short switching periods for VOC combustion[END_REF].

• If the feed of pollutant concentration 10 ( ) t ω is too lean, low heat released during the reaction leads to the extinction of the reactor. The temperature inside the reactor has therefore to be maintained over the lowest temperature specified as 450 K [START_REF] Ramdani | Reverse flow reactor at short switching periods for VOC combustion[END_REF].

• Ideal operation of such reactor is finally an operation without control. This is possible when the feed concentration is such that the reaction temperature falls within an envelope outlined by the two previous boundary temperatures. This ideal case is called autothermal operation and no control action is required. Except for the autothermal case, the controller to design has to compensate the influence of the input disturbance

( ) ad T t ∆
while preserving the stability of the system. On the other side, most of the reaction takes place at the inlet of the catalytic monoliths and is instantaneous [START_REF] Nieken | Limiting cases and approximate solutions for Fixed-bed reactors with periodic flow reversal[END_REF][START_REF] Ramdani | Reverse flow reactor at short switching periods for VOC combustion[END_REF]. Therefore, instead of the full temperature profile, only Ts at the entry 1 z of the catalyst zone is constrained:

1 450 ( , ) 600 K Ts z t K t < < ∀ (7)
Concerning the input disturbance

( ) ad T t ∆
, it varies between 0 and 30 K in real industrial use. But in order to prevent any accident, worst cases have to be evaluated. It is therefore assumed to vary randomly between 0 and 115 K (Fig. 2) and if no control is applied to the RFR, the hot-spot temperature cannot be maintained between both temperature limits (450 and 600 K). This clearly justifies the need for closed loop control in order to maintain the hot spot temperature within these 2 boundaries. Moreover, the stochastic input disturbance handled here has a much more realistic stochastic behavior than in the previous RFR control studies [START_REF] Budman | Control of a nonadiabatic packed bed reactor under periodic flow reversal[END_REF], Dufour et al. 2003). This disturbance and the temperature profile in the RFR are estimated on-line using a high gain observer based on three temperatures measurements [START_REF] Hammouri | Observer design for reverse flow reactor[END_REF]. The entire estimate state is injected in the LQR, whereas the MPC is only based on the estimated input disturbance. Simulation results allow comparing the closed-loop performances obtained with LQR and MPC.

COMPARISON OF TWO OPTIMAL CLOSED LOOP CONTROL APPROACHES

Since the hot spot temperature is constrained to be kept inside a specified region, some heating action may sometimes be required to avoid the RFR extinction. Some over times, some dilution may be required to avoid the over-heating of the catalyst elements, which would destroy them. From an optimal point of view, is constraints (7) are satisfied, no heating and no dilution are required.

LQR case

Approach used: The LQR is used here in a classical framework [START_REF] Levine | The control handbook[END_REF] where the two cost functions to minimize are given under the following: (8) with x is the process state vector, and u the manipulated variable vector. Q and R are the tuning parameters. The multivariable problem is decoupled here: when the feed of pollutant concentration is to rich, the solution of a first LQR is used for the minimum tuning of the dilution rate and no heating is possible. If the feed of pollutant concentration is to lean, the solution of another LQR is used to heat the system as less as possible, while no dilution takes place. Otherwise, no dilution and no heating are applied.

Simulation results: The output constraints (7) are satisfied at any time (Figure 3). Between 500s and 1550s, the low input disturbance (Figure 2) leads to a decrease of the temperature inside the RFR (Figure 3). LQR correctly tunes the internal heating (Figure 4) such that the temperature stays above the extinction temperature. No dilution is taking place until this time. After 1550s, rich feed (Figure 2) induces an increase of temperature inside the RFR. LQR then tunes the dilution rate such that the temperature is maintained below the maximum temperature and there is no more heating.

MPC case

Approach used: a special MPC was used here for such complex nonlinear PDE model and (Dufour et al., 2003b). The constrained optimal control problem to solve is under the form (k is the actual time, j is the future time index, Np is the prediction horizon):

(9) Compared to the LQR case, J is here specified as a real multivariable cost function since its aims to minimize, at the same time, both dilution and heating action under the input/output constraints (7).

( )

1 min max min max min ( , )
under some constraints

j k Np p u(k) j k p J f y (j) u : u u u y y y = + = + = ≤ ≤ ≤ ≤ ∑ ( ) ( ) ( ) ( ) 0 1 x' Qx u' Ru 2 k k k k k J ∞ = = +     ∑
Simulation results: The output constraints (7) are also satisfied at any time (Figure 5). Between 0s and 1300s, ad T ∆ is small (Figure 2) and extinction of the process is avoided feeding electrical power into the reactor (Figure 6). In the meantime, there is no cooling action and the maximum amount of gas is therefore treated as expected in these conditions. After 1500s, ad T ∆ becomes important and overheating of the process is avoided (see the upper bound constraint on Figure 5) due to the correct use of the cooling action. The drawback is that the controller may sometimes require both heating and cooling actions at the same time (at 3180s e.g.), which should not happen.

CONCLUSIONS ET PERSPECTIVES

In spite of a large influencing input disturbance due to the feed concentration, the temperature can be maintained by both observer based controllers inside the specified temperature envelope. Concerning the optimization performances, LQR leads here to better results than MPC since it requires less heating action while treating more gas: Q =83.4W and α =0.894 for LQR, Q =274.6W and α =0.849 for MPC (these mean values are calculated from t=0s to t=4500s). This difference is mostly due to the impact of the stochastic variations of the input disturbance over MPC. Indeed, the estimation of the disturbance ad T ∆ is directly used in the MPC, where it is assumed constant in the future. This may forces the MPC to over evaluate the need for heating and cooling. Combining this issue with the non minimum phase behavior of the process, the horizons tuning of the MPC is uneasy. In the meantime, a switching control structure is used for the LQR: therefore, heating and cooling action are explicitly decoupled and can not act at the same time as expected. The drawback is that this decoupling

Figure 1 :

 1 Figure 1: RFR description (left) and RFR counter current model (right).

  Figure 2: Adiabatic temperature rise

Figure 3 :

 3 Figure 3: Hot spot temperature inside the RFR (LQR case).

Figure 5 :Figure 4 :

 54 Figure 5: Hot spot temperature inside the RFR (MPC case).

Figure 6 :

 6 Figure 6: Cooling action by dilution (top) And heating action (bottom) (MPC case).

introduces a severe nonlinearity for stability analysis whereas the explicit multivariable MPC used here is more suitable than the LQR for general MIMO control problems. Moreover, in the MPC, output constraints are explicitly handled. Perspectives are concerned with the use of such tools for a new coupled thermal-catalytic process located in the LMSPC, using new catalyst products. Observation issue will also be tackled.