Carles Triguero
email: carlos-triguero@chimie-paristech.fr

François-Xavier Coudert
email: fx.coudert@chimie-paristech.fr

Anne Boutin

Alain H Fuchs

Alexander V Neimark
email: aneimark@rutgers.edu

Mechanism of Breathing Transitions in Metal-Organic Frameworks

Keywords: adsorption, metal-organic frameworks, structural transitions, elasticity, Monte Carlo simulation

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

2

Metal-organic frameworks (MOF) represent a new class of porous crystals, which has attracted physicists and material scientists worldwide due to their extraordinary physico-chemical and mechanical properties and potential applications in energy and biomedicine fields. Recent discovery of structural transformations between large pore (lp) and narrow pore (np) conformations in MOFs of MIL-53 type and others brought about a new physical phenomenon that represents a challenging problem for theoreticians [START_REF] Férey | Large Breathing Effects in Three-Dimensional Porous Hybrid Matter: Facts, Analyses, Rules and Consequences[END_REF][START_REF] Horike | Soft Porous Crystals[END_REF][START_REF] Coudert | Thermodynamic Methods and Models to Study Flexible Metal-Organic Frameworks[END_REF]. This phenomenon called breathing is displayed in abrupt changes of the framework volume triggered by adsorption of guest molecules that is explored to devise advanced adsorbents, drug delivery systems, sensors, and actuators [START_REF] Liu | Reversible Structural Transition in MIL-53 with Large Temperature Hysteresis[END_REF][START_REF] Beurroies | Using Pressure to Provoke the Structural Transition of Metal-Organic Frameworks[END_REF]. In this Letter, we suggest a multiscale cooperative mechanism of breathing transitions that involves a complex interplay of adsorption and elastic interactions on the level of crystal. This mechanism is explored with a "primitive" stochastic model of adsorption in flexible bistable frameworks that takes into account the major physical factors with the minimum input parameters and is capable of reproducing experimentally observed features. As such, we address the topical questions of the hysteretic nature of breathing transitions and the possibility of co-existence of lp and np phases in one crystal that was highly debated in the literature recently [START_REF] Boutin | The Behavior of Flexible MIL-53(Al) upon CH 4 and CO 2 Adsorption[END_REF][START_REF] Millange | Structural Effects of Solvents on the Breathing of Metal-Organic Frameworks: an In Situ Diffraction Study[END_REF][START_REF] Springuel-Huet | 129 Xe NMR Study of the Framework Flexibility of the Porous Hybrid MIL-53(Al)[END_REF]. This structure may oscillate between lp and np phases, which have a remarkable difference in unit cell volume of up to 40% [START_REF] Serre | An Explanation for the Very Large Breathing Effect of a Metal-Organic Framework during CO 2 Adsorption[END_REF][START_REF] Serre | Very Large Breathing Effect in the First Nanoporous Chromium(III)-Based Solids[END_REF] (see schematics in Fig. 1). What is truly extraordinary is that the transition from the larger volume lp phase to the smaller volume np phase is not necessarily accompanied with the release of guest molecules that would be expected for the normal "exhaling". For example, the equilibrium state of the MIL-53 crystal at 220 K in the absence of guest molecules is in the lp phase, and upon Xe adsorption there first occurs the transition from unloaded lp phase to loaded np phase [START_REF] Boutin | Breathing Transitions in MIL-53(Al) Metal-Organic Framework upon Xenon Adsorption[END_REF]. This transition is associated with a sharp uptake of Xe, from a loading of ~0.2 to ~2.5 molecules per unit cell, and a decrease of the crystal volume by ~25%.

Upon further increase of the gas pressure, adsorption gradually proceeds in the np phase up to a certain point, when the second, now "normal", breathing transiting occurs, from the np phase to lp phase. The sample abruptly "inhales", increasing the loading from ~2.7 to ~6.5 molecules per cell and expands, compensating for the volume lost upon the first lp-np transition. On the desorption pass, the reverse normal lp-np and abnormal np-lp "exhaling" transitions take place with a prominent hysteresis.

This enigmatic breathing phenomenon is engendered by guest-host adsorption interactions mediated by the elasticity of three dimensional host framework, which are currently poorly understood. The specific variations of the linker conformations in the lp and np phases during breathing transitions have been studied at the molecular level by Férey and co-authors, both experimentally (in-situ X-ray diffraction) [START_REF] Serre | An Explanation for the Very Large Breathing Effect of a Metal-Organic Framework during CO 2 Adsorption[END_REF] and using molecular simulation (single-point DFT calculations and forcefield-based dynamics) [START_REF] Salles | Molecular Dynamics Simulations of Breathing MOFs: Structural Transformations of MIL-53(Cr) upon Thermal Activation and CO 2 Adsorption[END_REF]. These works provide a useful insight into the chemistry of the transformation of linker bonds associated with the framework deformation. However, a knowledge gap exists between this molecular understanding and the question of how the adsorption of guest molecules induces the physical forces responsible for macroscopic expansion and contraction of the sample. The adsorption of molecules may induce a significant stress inside the elastic framework of the order of several to tens of MPa, as was recently shown [START_REF] Neimark | Structural Transitions in MIL-53 (Cr): View from Outside and Inside[END_REF]. The adsorption stress may be either positive or negative depending on the loading and may thus induce either framework expansion or contraction. In accord with the adsorption stress hypothesis [START_REF] Neimark | Stress-Based Model for the Breathing of Metal-Organic Frameworks[END_REF], the structural transformation occurs upon achieving a certain threshold stress that the respective phase cannot withhold. Although this hypothesis describes on a semi-quantitative level the breathing transitions driven by Xe and CH 4 adsorption on MIL-53(Al) [START_REF] Neimark | Stress-Based Model for the Breathing of Metal-Organic Frameworks[END_REF][START_REF] Boutin | The Behavior of Flexible MIL-53(Al) upon CH 4 and CO 2 Adsorption[END_REF], it is macroscopic in nature and does not shed light on the mechanism and dynamics of structural transformations on the unit cell level. As shown below, adsorption interactions in the elastic framework are long range and bring about a cooperative mechanism of cell filling: adsorption in one cell triggers adsorption in the 1D channel, which this cell belongs to, and this, in turn, triggers adsorption and the respective phase transformation in the whole 2D layer of channels.

From the geometrical standpoint, the lp and np phases of MIL-53 differ mainly by the shape of cells, which can be characterized by the angle θ of its rhombus cross-section; in lp phase, with a pore diameter of 9.5 Å, θ ~ 79°, while in np phase (with pore diameter of 3.5 Å) θ ~ 40° (Fig. 1). At the same time, the variation in the linker lengths between the two phases is of the order of 0.4 Å (just only 3.8%), making the length of the rhombus sides essentially invariant upon deformation of the structure. This strong constraint on the framework geometry plays a crucial role in the dynamics of the phase transformation, as the deformation of a single cell can be quantified by just one degree of freedom represented by the angle θ. The assumption of the preserved rhombus shape with the fixed length of sides imposes a strong correlation on the deformation of neighboring cells in order to ensure the lattice integrity. First, the cells along a channel must coherently deform, as two stacked rhombuses cannot match unless they have the same angle θ. Secondly, the channels connected by sides of fixed length must possess equal θ along the shear direction. This defines a 2D layer of cells inside which all the channels have the same cross section. Mismatches in the rhombus angle can exist along the direction perpendicular to the shear plane. Thus, the only possible mechanism of framework deformation is the layer-by-layer shear and the elementary deformation consists of the shear of the layer of cells in the direction perpendicular to the channel axis (represented in Fig. 1). This conclusion can be formally derived from the Saint-Vénant principle by writing the equations for the elastic compatibility conditions linking the strain in neighboring cells [START_REF] Shenoy | Martensitic textures: Multiscale Consequences of Elastic Compatibility[END_REF]. The elastic compatibility equations yield a long-range interaction, similar to that of martensitic transformations [START_REF] Rasmussen | Three-Dimensional Elastic Compatibility and Varieties of Twins in Martensites[END_REF]. In the case of the MIL-53 structure, the long-range interaction works as a rigid constraint requiring the similitude of the cell shapes within the 2D layer. This means that the phase transformations occur in a cooperative manner and necessarily involve the entire layer of cells, all of which must be in the same phase. As a consequence, the layer of cells, rather than the single cell, can be taken as a base unit of the framework mechanical model. The 3D framework is thus presented as a stack of 2D layers of identical cells that makes possible to formulate a onedimensional minimalistic model of the framework deformation that captures the main system properties with a minimum number of input parameters.

Within our primitive 1D model, two main variables describe the state of each cell layer at given external thermodynamic conditions. First, we quantify the cell deformation by one geometrical parameter, the symmetrized strain e, which is derived from discrete elasticity theory [START_REF] Ariza | Discrete Crystal Elasticity and Discrete Dislocations in Crystals[END_REF] by renormalization of the shear deformation strain so that the stable lp and np structures correspond to e = -1 and e = +1 respectively. The second variable is the adsorption loading n, or the mean number of guest molecules per unit cell in the layer. By averaging the loading over the whole framework as a function of the gas external pressure (or chemical potential), one obtains the adsorption isotherm. This quantity is the key observable of the system, as it is measured in isothermal adsorption experiments. The experimental adsorption isotherms typically display two sharp yet continuous transitions from empty lp phase to almost fully loaded np phase and from np phase to fully loaded lp phase. The loading capacities of both phases differ significantly and represent the main quantitative parameters determining the adsorption isotherm behavior. Our primitive model is discrete; it allows for adsorption of either 0, 1, or 2 molecules in the cell, that corresponds schematically to empty cell, np cell capacity, and lp cell capacity. This simplification can be easily generalized to introduce additional loading levels, or to describe adsorption by a continuous variable. The proposed model describes the interplay between guest adsorption and host framework deformation in terms of a Hamiltonian that depends on the loading n and strain e fields and is expressed per unit cell,

H (n,e) = c 0 2 e i ! s i ( ) 2 + "F 2 s i + c 1 2 #e i ( ) 2 +! n i , s i ( ) !! a n i , s i ( ) e i ! s i ( ) $ % & ' ( ) i * . ( 1
)
The first term in the right-hand side accounts for elastic deformation of individual cells with the local elastic energy is modeled as a biparabolic potential [START_REF] Müller | Model for an Elastic-Plastic Body[END_REF], see Fig. 2; s i = sgn(e i ) denotes the phase of the cell (s i = -1 in lp phase and s i = +1 in np phase), and (e i -s i ) is the local deformation from the reference state of the respective phase. The effective elastic constant c 0 is taken equal for both phases. The second term ΔF represents the difference in the free energy between empty np and empty lp phases, which is positive, reflecting the fact that the initial "dry" state corresponds to the stable lp phase. The third term corresponds to the non-local cell-cell elastic energy proportional to the strain gradient squared with ; it penalizes the formation of interfaces between lp and np phases by the interfacial energy of 2c 1, and also levels elastic deformations in neighboring cells of the same phase. The fourth and fifth terms determine the host-guest interactions as the energy of adsorption ε(n i ,s i ) at given loading n i and deformation e i expanded around the adsorption energy ε(n i ,e i ) in non-deformed reference lp or np state, e i = s i . This expression gives rise to a quantity of prominent physical significance: the adsorption stress induced on the host framework due to its interactions with the guest molecules, which is defined as

! a n i , s i ( ) = !"! n i , e ( ) / "e e=s i ,n i , ( 2
)
in line with the thermodynamic definition of the adsorption stress [START_REF] Salles | Molecular Dynamics Simulations of Breathing MOFs: Structural Transformations of MIL-53(Cr) upon Thermal Activation and CO 2 Adsorption[END_REF][START_REF] Ravikovitch | Density Functional Theory Model of Adsorption Deformation[END_REF]21]. The adsorption stress couples the host-guest interactions with the elastic deformation and accounts for the forces exerted by the guest molecules on the host framework. It can be either negative or positive depending on the loading and thus cause either elastic contraction or expansion. As such, the number of input parameters characterizing adsorption in our model is reduced to 8: 4 adsorption energies ! n i , e i ( ) and 4 adsorption stresses! a n i , s i ( ) , n i =1 or 2, and s i = ±1.

Modeling the dynamics of the coupled adsorption and deformation in the process of incremental stepwise variation of the chemical potential of the adsorbed gas, we make a further assumption that the local elastic relaxation of the framework occurs on a much smaller time scale than the establishment of adsorption equilibrium at a given chemical potential. This assumption allows us to determine the continuous elastic strain field e as a function of the discrete fields s and n by minimizing the Hamiltonian (1) with respect to e, !H / !e n,s = 0 . Due to the quadratic nature of the elastic potential, this minimization yields a system of linear equations, which are easily solved. As such, we discretize the model and each cell layer is characterized by its state variables s i = ±1 and loading n i = 0, 1, or 2, with the Hamiltonian being a function H(n,s).

In mimicking the adsorption process, we performed grand canonical Monte Carlo (MC) simulations of the 1D array of cell layers with periodic boundary conditions. In this molecular simulation, each element in the 1D system is a layer of MIL-53 unit cells; no atomistic detail of the framework is included. The whole system is in contact with a fictitious gas reservoir. The thermodynamic variables are the size of the system L, the temperature T and the chemical potential µ of the adsorbing gas. Two kinds of MC moves are considered, which lead to the change of the cell state: adsorption or desorption of one particle (n i →n i ±1) and the phase switch between lp and np structures (s i →-s i ). The exponential in the acceptance probabilities for particle insertions and deletions moves in the standard Metropolis scheme is written as exp[-β(H({n i±1 ,n j≠i ,s j })-H({n j ,s j })±µ)]. The probability of the phase switch includes the energy landscape of the host (1) as a biparabolic function with two very narrow wells, separated by a wide region where the behavior of the system is not elastic. We approximated this "no man's land" region of energy, which is not sampled in experiments, by an energy barrier E B (green horizontal line on Fig. 2) following the scheme employed by Kang et al. [START_REF] Kang | Dynamic Monte Carlo with a Proper Energy Barrier: Surface Diffusion and Two-Dimensional Domain Ordering[END_REF] in studies of activated diffusion. As such, we introduced the energy barrier directly in the acceptance probabilities of the phase change, which becomes 1, coupling parameter c 1 = 1.

proportional to exp[-β(H({n j ,-s i ,s j≠i })-H({n j ,s j })+E B /2)] if ∆H ≤ 2E B .
c 0 = 100 L = 2000 E B = 10.5 ∆F host = 5.0 ε(1,lp) = 0 ε(1,np) = -10 ε(2,lp) = -2 ε(1,np) = 3 σ a (1,lp) = -10 σ a (1,np) = -10 σ a (2,lp) = 10 σ a (1,np) = 10
As the result of simulation, we determined the experimentally measurable quantities: the adsorption isotherm as the average loading, the sample deformation through the phase composition, and the strain field, all as functions of the gas chemical potential. The set of model parameters is given in Table 1. These parameters were chosen in order to reproduce qualitatively some of the known characteristics of Xe adsorption in MIL-53(Al) [START_REF] Boutin | Breathing Transitions in MIL-53(Al) Metal-Organic Framework upon Xenon Adsorption[END_REF]: transition pressures, magnitude of adsorption-induced elastic deformation, and free energy difference between the phases. The coupling parameter c 1 was varied from 0 to 4. The adsorption-desorption isotherm <n>, obtained by gradually increasing and then decreasing the gas chemical potential, are reported in Fig. 3, along with the respective variation of phase composition x lp and adsorption-induced stress <σ a >. The calculated stepwise isotherm with two hysteresis loops parallels the experimental isotherm [START_REF] Boutin | Breathing Transitions in MIL-53(Al) Metal-Organic Framework upon Xenon Adsorption[END_REF]. The phase composition curve is also comparable to the one calculated using an earlier model [START_REF] Neimark | Stress-Based Model for the Breathing of Metal-Organic Frameworks[END_REF]. Each MC simulation run was 30,000 steps. The length of the simulation run plays a role of the equilibration time in adsorption experiments. While we ensured the establishment of adsorption equilibrium on the cell level, the phase transformation represents a rare event and, as in experiments, occurs in a hysteretic fashion. It is worth noting that it was our goal to reproduce the hysteretic nature of breathing transitions in our simulations rather than to avoid the metastable states and to confine ourselves on the conditions of phase equilibrium. This is why we introduce the energy barrier E B in the phase change acceptance probabilities. The longer the run the larger is the probability to cross the energy barrier associated with the phase change. Without this barrier, the system would follow equilibrium configurations and the adsorption isotherms were reversible, as shown in Fig. 3. For a given simulation length, the width of the hysteresis loop sharply increases with the height of the energy barrier. Chemical potential (µ) (see black hysteresis loop in Fig. 4b). In this limit, it is easier to overcome the energy barrier E B due to thermal fluctuations. In contrast, the existence of an interfacial energy effectively increases the energy barrier of the phase change of an isolated layer in a domain. As such, the phase transformation occurs in an avalanche manner:

c1 = 0 c1 = 1 c1 = 4 (a) (b) (c) G (2) (k) ∼ A e -k ξ
once the first nucleation event consisting of lp-np transformation in one layer happens, the further lp-np transformations will most likely proceed in the neighboring layers. This prevents microstructuration, leading to a small number of large domains (Fig. 4a) and sharper adsorption and desorption hysteresis loops (blue and red loop in Fig. 4b). Furthermore, the parameter c 1 is directly linked with the system correlation length, ξ [START_REF]The correlation length ξ is obtained by fitting the correlation function as G (2)[END_REF], as shown in Fig. 4c. In the limit of large c 1 , ξ diverges and only one domain would grow since the phase coexistence is strictly prohibited.

In conclusions, the elastic compatibility conditions on the level of crystal imply that the elementary lp-np transformation in the bistable MIL-53 structure necessarily involves the shear of the entire layer of unit cells and proceeds in either an avalanche and random manner. While our dynamic model predicts the possibility of phase co-existence in the process of lp-np transformation, additional information about cell-cell correlations is required to determine to which extent the experimental systems may follow this scenario. This information might come from space-resolved experimental technique showing the evolution of the cell parameters near an lp-np interface (if it exists), or from first-principle atomistic calculations of the elastic constants of lp and np phases.

The suggested model can be upgraded to include realistic adsorption capacities of the lp and np phases for a given guest-host system and further expanded to tree-dimensional anisotropic structures.

A.V.N. and C.T. acknowledge support from the Blaise Pascal International Chair award.

Figure 1 .

 1 Figure 1. Schematic representation of the MIL-53 framework after the first event of lp-np transformation, which involves in-plane shear of a 2D layer of cells in a direction orthogonal to the channel axis.

Figure 2 :

 2 Figure 2: Schematics of the host free energy landscape for the "dry" bistable framework, as a function of symmetrized strain e. Two regions of elastic deformation around the equilibrium lp and np structures (e = ±1)

Figure 3 :

 3 Figure 3: Breathing transitions during the adsorption-desorption cycle. (a) Adsorption isotherm, (b) adsorption stress, and (c) lp phase composition as a function of the chemical potential. Adsorption (red line), desorption (blue line), reversible isotherm obtained by ignoring the energy barrier E B in simulation (gray). Model parameters from Table 1, coupling parameter c 1 = 1.

Figure 4 :

 4 Figure 4: Effect of cell-cell coupling. (a) Number of domains in the system as a function of µ for c 1 = 1 (blue), and c 1 = 4 (red), upon adsorption (solid) and desorption (dotted). (b) Adsorption hysteresis loop for the lowpressure lp-np transition, for different values of c 1 . (c) Variation of the correlation length ξ with the coupling parameter c 1 .

Table 1 :

 1 Values of the model parameters used in this work. System size is adimensional, other parameters are given in units of k B T.

(


Alexander V Neimark
email: aneimark@rutgers.edu

François-Xavier Coudert
email: fx.coudert@chimie-paristech.fr

Carles Triguero

Anne Boutin

Alain H Fuchs

Isabelle Beurroies

Renaud Denoyel

Structural Transitions in MIL-53 (Cr): View from Outside and Inside

We present a unified thermodynamic description of the breathing transitions between large pore (lp) and narrow pore (np) phases of MIL-53 (Cr) observed during the adsorption of guest molecules and the mechanical compression in the process of mercury porosimetry. By revisiting recent experimental data on mercury intrusion and in-situ XRD during CO 2 adsorption, we demonstrate that the magnitude of the adsorption stress exerted inside the pores by guest molecules, which is required for inducing the breathing transition, corresponds to the magnitude of the external pressure applied from the outside that causes the respective transformation between lp and np phases. We show that, when a stimulus is applied to breathing MOFs of MIL-53 type, these materials exhibit small reversible elastic deformations of lp and np phases of the order of 2-4%, while the breathing transition is associated with irreversible plastic deformation that leads to up to ~40% change of the sample volume and a pronounced hysteresis. These results shed light on the specifics of the structural transformations in MIL-53 (Cr) and other soft porous crystals (SPC).

The structural transitions in metal-organic frameworks (MOF) known as breathing have been attracting a lot of attention from both the material scientists, who discovered and explored potential practical applications of this intriguing phenomenon, and the theoreticians, who attempted to establish the thermodynamics foundations and simulation models to understand its physico-chemical and mechanical mechanisms. 1,2,3,4,5 The most prominent breathing transitions were observed in the process of gas adsorption in some MOFs of the MIL-53 family. The MIL-53 framework is made of parallel one-dimensional M(OH) chains (M = Al 3+ , Cr 3+ ), linked together by 1,4-benzenedicarboxylate (BDC) ligands to form linear diamond-shaped channels that are wide enough to accommodate small guest molecules. 6,7,8 This structure may oscillate (or "breathe") between two distinct conformations called a large-pore phase (lp) and a narrow-pore phase (np; see Figure S1), which have a remarkable difference in cell volume of up to 40%. The conditions, at which the breathing of MIL-53 happens, have been widely studied earlier, and it was shown in particular to be triggered by changes in temperature, 9 as well as adsorption of some gases and fluid mixtures at room temperature, but not others. 10,11,12,13,14 Recent work has been successful in rationalizing the conditions for the occurrence of breathing in MIL-53, 15 by considering the issue from the point of view of thermodynamics in the osmotic ensemble, which is the statistico-mechanical ensemble adapted for the study of adsorption in a flexible host. 16 This thermodynamic analysis was then extended to fluid mixtures and to an investigation of temperature dependence, which allowed to predict complete pressure-composition and pressuretemperature phase diagrams for various adsorbates in MIL-53. 17,18,19 Recently, the authors 20 suggested a rationale for breathing transitions driven by guest molecules adsorbed inside the pores. It was shown that the adsorbed molecules induce a significant stress inside the elastic framework that can be of the order of several to tens of MPa. This adsorption stress, ! a , may be either positive or negative depending on the loading and, thus, it may induce either framework expansion or contraction. The breathing transitions from lp to np phase and back were treated as spontaneous non-equilibrium transitions occurring upon developing a certain threshold stress !*, which the respective phase cannot withhold. Based on this adsorption stress model and the threshold stress ansatz, the authors described on a semiquantitative level the breathing transitions driven by Xe and CH 4 adsorption on MIL-53 (Al). 20,21 In a recent publication, Beurroies et al. 22 demonstrated that the structural transitions between lp and np phases can be induced by uniform mechanical compression of the sample from outside in the process known as mercury porosimetry. 23 The powder of MIL-53 (Cr) particles was immersed into liquid mercury subjected to the external pressure. Upon the increase of pressure, mercury, being a non-wetting fluid, was forced to penetrate into interparticle voids of micrometer size, but it was not intruded into the internal pores of nanometer size, thus exerting a uniform pressure on the particles. The lp-np transition in the process of mercury intrusion and the reverse np-lp transition in the process of mercury extrusion occurred at significantly different pressures, giving rise to a pronounced hysteresis loop formed by intrusionextrusion curves.

In this work, we extended the thermodynamic approach based on the osmotic ensemble 16 coupled with the adsorption stress model 20 suggested for the description of adsorption-induced breathing transitions, to describe mechanically-driven transitions observed by Beurroies et al. with MIL-53 (Cr). 22 To this end, we revisited the mercury intrusion data and compared it with CO 2 adsorption and in situ XRD data collected on MIL-53 (Cr) by Serre et al. 10 We demonstrate that the magnitude of the adsorption stress applied inside the pores by adsorbed molecules, which is required for inducing the breathing transition, corresponds to the magnitude of the external pressure applied from the outside that causes irreversible plastic transformation between lp and np phases. This conclusion allows one to predict the adsorption-induced phase transformations of MOF based on the mechanical contraction experiments.

View from outside: What can one learn from the mercury intrusion-extrusion experiments?

In the mechanical compression experiments, a powder of MIL-53 (Cr) particles is placed in the sample cell of the mercury porosimeter, which represents a piston-cylinder assembly. To ensure that the particles are in the lp state the sample was initially immersed in methanol which is known to open the MOF structure. Then it was outgassed at 150°C to remove adsorbed species, cooled and transferred as quickly as possible in the mercury porosimeter cell for analysis. In the mercury apparatus, a primary vacuum is set lower than 5Pa. Then, a given amount of mercury is added into the sample cell on the top of the powder and the mercury pressure is gradually increased by a hydraulic pump. The cell volume underneath the piston is measured as well as the equilibrium external pressure applied through the piston to the mercury and thus transmitted to the powder particles. The results of the mercury intrusion-extrusion cycles are presented in Figure 1. Note that the raw data 22 was corrected for the mercury compressibility and the pressure transducers hysteresis using blank calibration experiments on mercury intrusion-extrusion inside the empty sample cell. This correction is necessary for quantitative analyses of the experimental data at high pressures, especially in the region of hysteresis. The graph in Figure 1 gives the volume of mercury V penetrated at the given external pressure P into the volume occupied by the powder; the external pressure in the initial state is null, P = 0. The consecutive stages of the intrusion-extrusion process are shown by different colors: the blue curve depicts the primary intrusion obtained by raising P from 0 to 410 MPa, the red curve depicts the primary extrusion obtained by reverting the process at P = 410 MPa and reducing P to 0.34 MPa, the violet curve depicts the secondary intrusion obtained by reverting the process at P = 0.34 MPa and raising P again to 410 MPa, and finally the orange curve depicts the secondary extrusion obtained within the same pressure limits as the primary one. The almost ideal coincidence of the primary and secondary extrusion curves and co-linearity of the primary and secondary intrusion curves at P > 10 MPa show the reliability of the presented experiment.

The mercury porosimetry technique is a standard tool for pore size characterization of porous solids. 23 It employs the fact that being a non-wetting fluid, mercury penetrates into the pores only when the external pressure exceeds the capillary pressure determined by the Laplace-Washburn equation, P = 4 ! cos "/ D, where D is the effective pore diameter, ! is the mercury surface tension (!= 0.485 N/m), and " in the contact angle, which for most solids varies from 130° to 145° [14]. Thus, in order to intrude mercury into micrometer-wide pores, the external pressure should exceed 0.2 MPa, and, respectively, filling of micropores (D < 2 nm) may occur only at an external pressure exceeding ~600 MPa. These simple estimates help understand the different steps in the mercury intrusion-extrusion cycle. The first step on the primary intrusion curve at low pressures P < 0.2 MPa (Figure 1, point A) characterizes the powder compaction into a dense random packing without mercury penetration into the gaps between micrometer size particles of MIL-53 crystals. The second step, at 0.2 MPa < P < 10 MPa (from point A to point B), reflects mercury filling of the interparticle voids in the dense packing of particles without penetration into intraparticle pores. Indeed, the pores in MIL-53 are in the micropore range (< 2 nm) and cannot be filled at pressures smaller than 600 MPa that exceeds the maximum pressure of 410 MPa achieved in given experiments.

Particle size analysis.

From the range of pressures, 0.2-10 MPa, corresponding to the interparticle filling, one can estimate the particle size distribution using the conventional method of Mayer and Stowe. 24 The calculated differential particle size distribution is presented in Figure 2 as the mass fraction of particles of given diameter. In the basic Mayer-Stowe equation, P = ! " / D p , which relates the intrusion pressure P and the particle diameter D p , the mercury surface tension " was taken as 0.485 N/m, and the dimensionless Mayer-Stowe constant !, that typically varies in the range of 6-13, was taken as 10. [START_REF] Pabst | Characterization of particles and particle systems[END_REF] The average particle diameter is D p = 2.7 µm, and the width of the distribution is ~ 0.8 µm. This is in agreement with the results of Khan and Jhung, where scanning electron microscopy (SEM) images showed relatively homogeneous particles of size 2-3 µm, though on a different batch of MIL-53 particles. [START_REF] Khan | [END_REF] 1.2 Elastic deformation and bulk modulus estimate. Neglecting mesopores, which may exist in small amounts due to possible defects and imperfections of MIL-53 crystals, we assume that mercury does not penetrate the intraparticle pores in the whole range of external pressures in given experiment. Upon the filling of interparticle pores, mercury exerts the uniform hydrostatic pressure on the particles and the further gradual increase of the intruded volume at 10 MPa < P < 35 MPa reflects the elastic compression of particles, which are in the lp phase. This elastic compression should be reversible. From the slope of the intrusion curve in the region of linear elastic deformation, one can determine the volumetric elastic modulus (or bulk modulus) K lp of the lp phase of MIL-53 crystals, as K = V 0 (∂P/∂V) with V 0 being the reference volume of the respected phase. We have calculated the bulk modulus of the lp phase from the second mercury intrusion curve, for which the linear regime is clear (see Figure S2). The value obtained is K lp ~ 2 GPa, which is in line with the values of bulk modulus reported in the literature for other MOF structures, which typically fall in the range 1.5-12 GPa. 27,28 (It is worth noting that HKUST-1 was reported to have a much higher bulk modulus than other members of the MOF family, at ~ 30 GPa. 29 ) Similarly, the bulk modulus of the np phase can be estimated from the primary and secondary extrusion curves, which overlap reasonably in the region of elastic deformation at P > 30 MPa. Thus found value of K np ~ 10 GPa is larger than that of the lp phase. Indeed, the np phase, which is denser than the lp phase, is also less compressible. We can also estimate the elastic domain of deformation, i.e. the maximal variation of the crystal volume observed within the region of elastic deformation. The unit cell volume variation in the elastic regime for both phases is small (few percents), as indicated in Table 1.

Structural lp-np and np-lp transitions.

The region of elastic contraction terminates with a prominent step at P * ~ 55 MPa, which was treated as the transition from the lp to the np phase. This transition is highlighted on Fig. 1 with a vertical arrow. The volume of intruded mercury during this step, of 394 Å 3 /u.c., agrees within ~10% with the difference in the volumes of lp and np phases as obtained from crystallographic data, namely 440 Å 3 /u.c. The difference between these two values, of the order of 10%, may be ascribed to the continued existence of a fraction of lp phase in the experimental sample, as was observed upon Xe and CO 2 adsorption. 21,30 According to the ansatz that the lp-np transition occurs upon achieving in the lp phase a certain threshold negative stress !*, the value of the external pressure P * lp of the lp-np transition gives the estimate of the maximum negative stress !* lp , which the lp phase can withhold, namely P * lp = " (!* lp " ! 0 lp ), where ! 0 lp is the pre-stress inherent to the dry material in the absence of mechanical constraint at P = 0. Moreover, the authors hypothesized that the threshold stress !* depends on the particle size, the smaller the particle the smaller the threshold stress. 20 Following this hypothesis, we can associate the width of the lp-np transition step P * lp = 55 (± 15) MPa with the distribution of the threshold stress due to the particle size distribution.

Thus, using the particle size distribution determined from the interparticle penetration region, one can correlate the threshold pressure of the lp-np transition and the particle size.

Upon further increase of pressure the intruded volume increases linearly and reflects the elastic compression of particles, which are now in the np phase. The elastic compression is reversible as seen from the extrusion curve, which coincide with the intrusion one within experimental errors. From the slope of the intrusion-extrusion curve in this region, one can determine the volumetric elastic modulus of the np phase. The elastic expansion of the np phase in the process of mercury extrusion terminates with a stepwise np-lp transition, which occurs in the range of pressures 0 < P < 15 MPa. Similarly to the lp-np transition, according to the threshold stress ansatz, the np-lp transition occurs upon achieving a certain threshold (now positive) stress, P * np = 10 (± 10) MPa, which may depend on the particle size. The np-lp transition is, again, highlighted on Fig. 1 with a vertical arrow.

Note that the external pressure in the intrusion experiment equals the negative stress calculated from the pre-stress in the reference evacuated sample at P = 0. For the lp phase, the situation is straightforward, since this reference state is the starting point of the intrusion experiment. However, the np phase of the MIL-53 (Cr) structure is unstable at the temperature of these experiments (around 300 K), so that no such reference state exist. To remedy the situation, we can considered a hypothetical np state, which would be obtained by extrapolation of the linear elastic expansion regime down to P = 0, as the reference state for thermodynamic calculations. This extrapolation is shown in Figure 1 by a dotted line (lower panel; the reference state is marked as point C).

View from inside: adsorption-driven structural transitions.

Gas adsorption was the first physical stimulus, for which breathing transitions in MIL-53 were observed, well before it was demonstrated that temperature and mechanical pressure could also trigger the lp-np structural transition. Adsorption of many gases inside MIL-53 (Al, Cr) is characterized by step-wise isotherms indicative of breathing transitions between np and lp phases, which occur at certain gas pressures depending on the nature of the guest molecules and the temperature. The criteria for observing the breathing phenomenon in terms of MIL-53

properties and fluid/solid interactions have been defined from the thermodynamic considerations.

It was shown that breathing is a general phenomenon that occurs within a certain temperature range for all guest molecules. 19,21 Noteworthy, in the case of MIL-53 adsorbents, the two consecutive breathing transitions take place, lp-np at low gas pressure and np-lp at high gas pressure. The reverse lp-np and np-lp transitions occur in the process of desorption in a hysteretic manner. Analyses of adsorption isotherms provide the quantitative information about the breathing transitions in terms of the difference of adsorption capacities of np and lp phases and, combined with in situ XRD patterns, the difference of unit cell geometries.

In the pioneering paper on CO 2 adsorption in MIL-53 (Cr), Serre et al. reported in situ Xray diffraction experiments performed during adsorption/desorption cycles at room temperature. 10 Revisiting this data on the adsorption-driven breathing obtained on the same structure as the mercury porosimetry experiments allows us to compare quantitatively the adsorption-driven breathing to that driven by mechanical pressure. Note that this does not mean necessarily the phase coexistence within one MIL-53 (Cr) crystal;

Analysis of in

in a powder sample, some crystals can be in the lp phase and the others in the np phase.

In addition, we also see that the position of the XRD peaks of each structure shifts slightly during adsorption, which is an indication of the elastic deformation of each phase. This deformation is due to the adsorption stress created by the presence of CO 2 in the pores of the material. By using a simple geometric model of the main feature of MIL-53 flexibility, we were able to map these small displacements of the XRD peaks onto variations of unit cell volume. The details of this model (and equations) are given in Supplementary Information.

The unit cell volume changes thus calculated are reported for the lp phase in Figure 4 (the plot for the np phase is less instructive as the amplitude of variations is much smaller). The relative variation of volume upon adsorption, of around 1.5%, is comparable to that measured from the mercury intrusion experimental data, indicating a good agreement between the two completely different sets of experimental results (and validating the geometric model used to process the XRD peak positions). Moreover, it can be noted that, because the measurements used for Figure 4 correspond to the high-pressure range of the adsorption and desorption isotherms, we are in a region of positive and increasing adsorption stress: as CO 2 pressure increases, more molecules are packed into the nanopores of MIL-53, leading to an outward stress. These experimental measurement are thus in agreement with the picture provided by the stress model proposed in our earlier work. 20

Analysis of the adsorption isotherms.

We have further applied the stress model 20 to the adsorption data presented in Figure 3 left, combining the information we can extract from the adsorption isotherms with the data extracted from the mercury intrusion experiments, such as the bulk moduli of the lp and np phases and the values of threshold stresses corresponding the breathing transitions. All this information leads to a unique set of the parameters for the stress model of MIL-53 breathing, which are given in Table 1. It is noteworthy that, while we have previously used the stress model to provide for a qualitative understanding of the Xe-and CH 4driven breathing of MIL-53 (Al), introducing information from mechanical experiments has now enabled us to calculate quantitatively the actual adsorption stress in the material.

Unified description of sample deformation engendered by external pressure and adsorption inside the pores.

The equilibrium elastic deformation of a porous solid immersed in a fluid (gas or liquid) thermostated reservoir, which may contain one adsorbing species, is determined from the following fundamental equation,

where is the osmotic ensemble thermodynamic potential of the combined solid+fluid system at the given fluid pressure p ext and chemical potential µ of the adsorbate, V is the current sample volume, V 0 is the sample volume in the evacuated (or "dry") state at and , which is considered as the reference state for thermodynamic analysis, is the solid free energy in the dry state, K 0 is the bulk modulus (or volumetric elastic constant) of the solid in the dry state, and the grand thermodynamic potential of the adsorbate at given µ and V. Eq. 1 implies that the changes of the solid free energy with respect to the dry state are accounted by a harmonic elastic potential . In this simplest consideration, we ignore a possible anisotropy of the solid and characterize the elastic deformation of the framework by the volumetric strain .

(

) 2
The threshold stress ansatz implies that the regions of elastic deformations are limited by the minimum negative stress !* lp and respectively minimum negative strain #* lp in the case of lp phase and the maximum positive stress !* np and respectively maximum positive strain #* np in the case of np phase. Ignoring the necessary deviations from the linear elasticity in the vicinity of threshold points, the strain-stress dependencies for the lp and np phases are presented by two straight lines terminated at the threshold points, Figure 5. The slopes of these lines correspond to the respective bulk modulus K lp/np of the lp and np phases. The mechanically driven transitions in the process of mercury intrusion-extrusion occurring upon achieving the threshold stress are depicted by the vertical lines in Fig. 5. The resulting hysteresis loop is a mirror image of the experimental mercury intrusion-extrusion cycle in Figure 1 formed by the elastic deformation curves bounded by the vertical lines corresponding to the lp-np and np-lp transitions in the crystals of the mean size.

The equilibrium deformation is determined by minimization of the thermodynamic potential with respect to the sample volume V at given p ext and µ that brings about a general relation between the volumetric strain # and stress ! in the form of the Hooke law, ,

Here, we introduced adsorption stress 32,33 , (

as the derivative from the grand thermodynamic potential with respect to the sample volume V at given external conditions of constant µ and p ext . Note that for one component fluid p ext and µ are related through the equation of state; that is in adsorption experiments, p ext equals the gas pressure p, and in the mercury intrusion experiments p ext equals the applied pressure P.

Eq. 3 constitutes the basis for a unified treatment of the mechanical and adsorption-driven deformation of MOFs. Indeed, the sample deformation in the process of adsorption should be similar to that in the process of mechanical compression provided the equality of the stress. We may conclude that assuming that both the adsorption experiment and compression experiment begin from the same initial dry evacuated state, the same sample deformation during adsorption and compression is achieved, when

Eq. 5 allows us to predict the adsorption driven elastic deformation based of compression experiments, from which the stress-strain dependence (Fig. 5) can be derived, and vice versa.

This conclusion implies that the adsorption stress can be calculated independently, for example, by employing an empirical Langmuir adsorption equation as it was done in our previous works for Xe and CH 4 adsorption on MIL-53 (Al). 20,21 The variation of the cell volume of the MIL-53 (Cr) crystal in the process of CO 2 adsorption is a function of the adsorption stress that can be calculated based on the stress model. 20 Following our earlier treatment of Xe sorption in MIL-53 (Al), we approximated the adsorption isotherms on the lp and np crystals with the Langmuir adsorption equation, , with the unit cell capacity N 0 and the Henry constant K H determined from the fitting of the experimental data. 10 The adsorption stress was calculated according to Eq. 4, which in the case of Langmuir adsorption isotherm is directly related to the gas pressure, 20

In addition to N 0 and K H , the adsorption stress (6) depends on the variations of these quantities with respect to the cell volume, and It is noteworthy, that while is apparently positive, in general negative. This gives rise to a nonmonotonic variation of the adsorption stress in the course of adsorption, as shown in ref. 20.

In the Henry region at low pressures, is negative and decreasing, that corresponds to the sample contraction. In the saturation region at high pressures, increases and may even become positive, leading to the sample expansion. Based on the threshold stress ansatz, the values of and were determined from Eq. 5 equating the characteristic mercury pressures of lp-np and np-lp transitions observed in the porosimetry experiments (vertical lines at the threshold pressures in Figure 5) and the adsorption stress (Eq. 6) of the respective transitions in the adsorption experiments (Figure 3, left). The parameters on the CO 2 adsorption isotherms on lp and np phases are given in Table 2.

Noteworthy, the system scans the different regions on the stress-strain diagram (Figure 5)

during the adsorption-desorption cycle and during the intrusion-extrusion cycle. Indeed, in the compression experiments the stress exerted by the external pressure is always negative, while the adsorption stress may be positive due to the packing effects at high loading. At the same time, the negative adsorption stress is limited in the adsorption experiments due to the competition between the guest-host attractive interactions and the repulsive guest-guest interactions.

However, the most important conclusion is that in the process of external compression and in the process of adsorption, the host framework experiences the same strain-stress elastic behavior and the lp-np transitions occur upon achieving the same threshold stress, which can be either directly measured in the mechanical experiments or calculated from the adsorption isotherms.

Conclusions: how soft are Soft Porous Crystals?

In the recent literature, flexible MOFs which exhibit adsorption-driven elastic deformation, gating, and breathing phenomena, of which MIL-53 is one of the most characteristic examples, were named soft porous crystals (SPC). 5 Our study shows that softness of SPCs can be very different. Namely, one should distinguish between the structures like the swelling MIL-88 34 , or interpenetrated dynamic frameworks, 35 which change their volume gradually with respect to the applied stimulus of increasing magnitude, and the structures like breathing MOFs of MIL-53 type considered here. These latter structures feature regions of small reversible elastic deformations within relatively wide ranges of variation of the stimulus terminated by large irreversible plastic deformations occurring in a stepwise manner, which lead to hysteretic phase transformations.

The mechanical compression experiments with mercury porosimeters provide not only qualitative but also valuable quantitative information about the mechanical properties of lp and np phases and the breathing transitions, including elastic modules of respective phases, the extent of the elastic deformation, and threshold stresses leading to the phase transformations, as well as the particle size distribution and correlations between the threshold stress and the particle size. The quantitative estimates extracted from the mercury porosimetry data on MIL-53 (Cr) sample correlate with the results of interpretation of the independent data on CO 2 adsorption isotherms and in situ XRD patterns collected in the process of CO 2 adsorption. The combined body of data from adsorption and mechanical experiments provides a solid foundation for a reliable theoretical consideration of the breathing transitions in SPCs.

The osmotic ensemble thermodynamic analysis based on the assumption of the additive contributions of the adsorption guest-host interactions and the host-host interactions treated within the framework of linear elasticity into the guest-host free energy allows one to describe in a unified manner the mechanically and adsorption driven deformation of SPCs. The adsorption stress exerted by the guest molecules adsorbed inside the pores causes the same framework deformation as the externally applied pressure of the same magnitude. The threshold stress ansatz implies that the lp-np structural transitions take place abruptly upon achieving the limit of elastic stability of the respective phase, similarly to the transition from elastic to plastic deformation. Noteworthy, while the volume changes between the lp and np phases of MIL-53 are quite substantial (~40%), the elastic deformations are limited to 2.7% (for lp) and 3.6% (for np). The coherent quantitative analysis of the lp-np transitions in MIL-53 (Cr) based on three types of measurements, mercury porosimetry, adsorption isotherms, and in situ XRD, confirms the validity of the proposed theory.

The proposed thermodynamic treatment of adsorption deformation is not limited to ordered SPCs considered here and can be adopted to any compliant porous solids of different origin, ranging from zeolites 32 and mesoporous adsorbents 36 to coal seams in the process of carbon dioxide sequestration at geological conditions 37 .

  -situ XRD patterns. The isotherms and pressure-dependent XRD patterns from Ref. 10 are reproduced in Figure3, for one adsorption/desorption cycle. The step seen on the adsorption isotherm can be attributed to the np-lp transition, i.e. the second, high pressure breathing transition upon adsorption. The first, low-pressure lp-np transition happening around 0.3 bar31 is not seen here. The evolution of the X-ray diffraction pattern upon adsorption and desorption shows a clear phase transition between the np and lp structures, respectively characterized by the peaks at 2! around 5° and 3.9°. The transition can clearly be seen as a jump from one structure to another, with no intermediate XRD peak being observed. Moreover, it can be seen that, within a narrow range of pressure around the transition, the two phases can coexist in the experimental sample, as the two characteristic XRD peaks are observed simultaneously.

Figure 1 .

 1 Figure 1. Intrusion-extrusion experiments, featuring the volume variation of the sample (#V) as

Figure 2 .

 2 Figure 2. Particle size distribution (in black) and Gaussian fit (in red).

Figure 3 :

 3 Figure 3: CO 2 adsorption and desorption isotherms at 304 K (left) and structural transformations

Figure 4 :

 4 Figure 4: Unit cell volume variation for the lp phase upon CO 2 adsorption and desorption at

Figure 5 :

 5 Figure 5: Variation of the unit cell volume in the elastic regime, including threshold pressures

3 Table 2 .

 32 Model parameters for adsorption of CO 2 on MIL-53 (Cr) at 303 K: Langmuir parameters (K H and N 0 ) are taken from fitting the experimental data,10 and their derivatives with respect to unit cell volume were optimized to reproduce the threshold stress corresponding to the lp-np and np-lp transitions in the mercury intrusion measurements (Figure1) at the characteristic experimental pressures of the lp-np and np-lp transitions in adsorption measurements (Figure3). Pressure transitions for the first transition, not visible in data from ref. 10, were approximated to 0.2 and 0.3 bar, bracketing the experimental observations by calorimetry.31

Table 1 .

 1 Mechanical parameters for the two phases of MIL-53 (Cr): reference state volume, V 0 , as determined from crystallography; threshold stress (P*) and bulk modulus (K) extracted from mercury intrusion/extrusion; relative deformation #V/V0 during the mercury intrusion experiment,22 calculated as #P/K where #P is the pressure range of existence of each phase.

	Sampled in desorption
	Sampled in adsorption

Acknowledgment. A.V.N. acknowledges the ARO grant W911NF-09-1-0242 and the Région Île-de-France for the support via a Blaise Pascal International Research Chair, administered by the Fondation de l'École normale supérieure.

Supporting Information Available: Representation of the MIL-53 lp and np structures, fits used to determine the bulk moduli of the two phases, and detailed description of the geometrical model used for interpretation of in situ X-ray diffraction data. This material is available free of charge via the Internet at http://pubs.acs.org.


Dario Corradini

François-Xavier Coudert
email: fx.coudert@chimie-paristech.fr

Rodolphe Vuilleumier
email: rodolphe.vuilleumier@ens.fr

Carbon dioxide transport in molten calcium carbonate occurs through an oxo-Grotthuss mechanism via a pyrocarbonate anion

Carbonatite melts, mostly composed of XCO 3 (X = Li 2 , Na 2 , K 2 , Mg, Ca), have recently received a good deal of attention [START_REF] Gaillard | Carbonatite Melts and Electrical Conductivity in the Asthenosphere[END_REF][START_REF] Jones | Carbonate Melts and Carbonatites[END_REF]. In fact they have been proposed to be the most conductive phase in Earth's upper mantle. Carbonatite melts may determine the mobility and long-term storage of deep carbon in Earth, especially if they are found in the lower mantle [START_REF] Kaminsky | Nyerereite and nahcolite inclusions in diamond: evidence for lower-mantle carbonatitic melts[END_REF][START_REF] Stoppa | Nyerereite from carbonatite rocks at Vulture volcano: implications for mantle metasomatism and petrogenesis of alkali carbonate melts[END_REF][START_REF] Kaminsky | Mineralogy of the lower mantle: A review of 'super-deep' mineral inclusions in diamond[END_REF]. Estimations based on conductivity data have led to propose a 0.035-0.35% content of carbonatites by volume in the asthenosphere, in line with the estimated CO 2 content in magmas [START_REF] Gaillard | Carbonatite Melts and Electrical Conductivity in the Asthenosphere[END_REF]. Carbon dioxide from geological formations is mostly released into the atmosphere during volcanic activity, although the efficiency of CO 2 degassing remains poorly known [START_REF] Dasgupta | The deep carbon cycle and melting in Earth's interior Earth Planet[END_REF]. Despite the importance of CO 2 species for the deep carbon cycle [START_REF] Dasgupta | Melting in Earth's deep upper mantle caused by carbon dioxide[END_REF], the reactivity, speciation and solvation structure of CO 2 in carbonatite melts have so far remained unexplored.

New technologies capable of reducing the carbon footprint are highly demanded, in particular in order to limit the impact of CO 2 on the environment. One of the strategies considered in the context of carbon capture and valorisation consists in its dissolution in a molten carbonate medium and its subsequent electroreduction to CO [START_REF] Chery | CO 2 electrochemical reduction into CO or C in molten carbonates: a thermodynamic point of view[END_REF][START_REF] Chery | Thermodynamic and experimental approach of electrochemical reduction of CO 2 in molten carbonates[END_REF][START_REF] Kanai | Mass transfer in molten salt and suspended molten salt in bubble column[END_REF]. Therefore the interaction between CO 2 and carbonates is relevant for environmental science and engineering beyond geochemistry.

In the context of the electrochemical reduction of CO 2 , it has been observed that the solubility of CO 2 in the carbonate electrolyte is higher than what one could expect simply from x l CO 2 = k H P CO 2 , where x l CO 2 is the concentration of CO 2 in the liquid, k H is Henry's constant and P CO 2 is CO 2 partial pressure. In order to explain this increased solubility, a reaction of CO 2 with the carbonate anion has been hypothesised:

CO 2 + CO 2- 3 ⇌ C 2 O 2- 5 (1)
with the product of the reaction being named pyrocarbonate (or dicarbonate) [START_REF] Chery | CO 2 electrochemical reduction into CO or C in molten carbonates: a thermodynamic point of view[END_REF][START_REF] Chery | Thermodynamic and experimental approach of electrochemical reduction of CO 2 in molten carbonates[END_REF][START_REF] Kanai | Mass transfer in molten salt and suspended molten salt in bubble column[END_REF]. Taking K pyro as the equilibrium constant of this chemical reaction, one can then define an apparent Henry's constant k app H = (1 + K pyro )k H . Thus it is possible that the CO 2 /pyrocarbonate equilibrium controls the solubility of CO 2 in molten carbonates. While this equilibrium can contribute to a larger CO 2 uptake in molten carbonates, it is not clear whether this is favourable for electrochemical processes, as the CO 2 involved in pyrocarbonate may not be readily accessible depending on the rate of interconversion. It is thus crucial to obtain a better understanding of CO 2 solvation in molten carbonates at the molecular level.

A limited number of previous studies have provided some support to the possibility of the formation of pyrocarbonate. The existence of a CO 2 /pyrocarbonate equilibrium has been first hypothesised by Claes et al. on the basis of experiments performed on the Li-Na-K carbonate eutectic mixture [START_REF] Claes | Solubility and Solvation of Carbon Dioxide in the Molten Li 2 CO 3 /Na 2 CO 3 /K 2 CO 3 (43.5:31.5:25.0 mol-%) Eutectic Mixture at 973 K. I. Experimental Part[END_REF]. This has been also supported by gas phase calculations [START_REF] Peeters | Solubility and Solvation of Carbon Dioxide in the Molten Li 2 CO 3 /Na 2 CO 3 /K 2 CO 3 (43.5:31.5:25.0 mol-%) Eutectic Mixture at 973 K. II. Theoretical Part[END_REF].

It has later been shown, by density functional theory (DFT) calculations, that in the gas phase C 2 O 2-

5

and larger oligomers C n O 2- 2n+1 are structurally stable with their energetics strongly dependent on the counterion [START_REF] Burna | Density functional theory (DFT) calculations on the structures and stabilities of [C n O 2n+1 ] 2and [C n O 2n+1 ]X 2 polycarbonates containing chainlike (CO 2 ) n units (n=2-6; X=H or Li)[END_REF]. Polymeric phases of CO 2 have been also predicted by DFT [START_REF] Frapper | Search for New Allotropic Forms of Carbon Dioxide and Carbon Disulfide: A Density Functional Study of CX 2 -Based Oligomers (X= O, S)[END_REF]. Zhang et al. have proposed evidence of the presence of pyrocarbonate in Raman spectroscopy experiments performed exposing molten carbonates mixtures to a CO 2 atmosphere [START_REF] Zhang | First spectroscopic identification of pyrocarbonate for high CO 2 flux membranes containing highly interconnected three dimensional ionic channels[END_REF]. Finally, pyrocarbonate has been also invoked to explain NMR spectra in aqueous carbonate solutions [START_REF] Zeller | The Hidden Equilibrium in Aqueous Sodium Carbonate Solutions -Evidence for the Formation of the Dicarbonate Anion[END_REF].

Recently, two first principles molecular dynamics (FPMD) studies have addressed the properties of molten carbonates [START_REF] Vuilleumier | Structure, equation of state and transport properties of molten calcium carbonate (CaCO 3 ) by atomistic simulations[END_REF] and of CO 2 in silicate melts [START_REF] Vuilleumier | Carbon dioxide in silicate melts at upper mantle conditions: Insights from atomistic simulations[END_REF]. The first has provided relevant information for the study of liquid CaCO 3 , until then poorly known by experiments, such as the liquid structure, density, atomic vibration motions, diffusion coefficients and electrical conductivity. In the second one, the speciation of CO 2 in basaltic and kimberlitic melts in the CaO -MgO -Al 2 O 3 -SiO 2 system has been investigated. Importantly,

the formation of C 2 O 2- C O O C O O O 2- + C O O O 2- + C O O C O O O C O O 2- (a)

RESULTS AND DISCUSSION

Formation of pyrocarbonate anion

We study by FPMD simulations the behaviour of one CO 2 molecule in the CaCO 3 solvent.

From the analysis of the simulation trajectory, we observe the spontaneous formation of the pyrocarbonate ion C 2 O 2- 5 out of the reaction of CO 2 with the carbonate anion CO 2- 3 , as shown in Fig. 1

(a).

There is a rapid exchange between the species with multiple instances of formation and dissociation of the pyrocarbonate ion (79 events each in the simulation time span of 65 ps). In order to define the molecular species out of the FPMD simulation trajectory, it is necessary to use a distance criterion, usually determined by looking at the minimum of the pair radial distribution functions (RDFs). We have used a "hysteresis" With the mentioned criterion, we observe the presence of pyrocarbonate in the simulation box during 24.4% of the total simulation time. From this we can estimate a free energy of formation ∆A = RT ln((1 -p pyro )/p pyro ) = 16.8 ± 1.4 kJ.mol -1 (RT = 14.7 kJ.mol -1 at T = 1773 K). The minimal/maximal intervals for the continual presence of CO 2 or pyrocarbonate are respectively 0.024/4.84 ps and 0.016/1.52 ps. By analysing the histograms of the durations of CO 2 or pyrocarbonate intervals, see Supplementary Figure 1, we estimate a CO 2 average life time τ CO 2 = 0.88 ps and a pyrocarbonate average life time τ pyro = 0.28 ps.

The snapshots in Fig. 2(a), 2(b) and 2(c) show respectively the approach of CO 2- 3 to CO 2 , the formed pyrocarbonate, and its dissociation to CO 2- 3 and CO 2 . As we can see in Fig. 2(d) the identity of the C atoms belonging to CO 2 or pyrocarbonate varies along the simulation trajectory. We look then at all the pyrocarbonate formation/dissociation events and we keep track of the C atom part of the CO 2 molecule before the formation of pyrocarbonate and after its dissociation. We find that 72.2% of the times the same C is part of CO 2 before and after a pyrocarbonate ion time interval. Conversely, 27.8% of the times a different C atom composes CO 2 after the dissociation of the pyrocarbonate.

To go further, we measure the correlation between events in which C has remained the same and between events in which C has changed. Given an event in which C has remained the same before the formation and after the dissociation of pyrocarbonate, we measure the frequency of occurrence of n (after the first) consecutive events of the same kind. We measure the same for the events in which instead the C has changed before the formation and after the dissociation of pyrocarbonate. The results are shown in Fig. 2(e). Therein we also plot the theoretical curves γ(1 -γ) n for "stay" events and (1 -γ)γ n for "change" events, with γ = p stay→change = p stay (0), which assume no correlation between successive events. The probability of n consecutive "stay" events is a slowly decaying function of n.

Conversely, the probability of consecutive "change" events decays very rapidly, with a small probability of another change in C after the first one and a very small probability of having two changes after the first one. Therefore the cascading of C transfers, through a rapid succession of jump events (as in Fig. 1(b)) appears unlikely. Hence, in the carbonate melt, the transport of CO 2 (free or bound in pyrocarbonate form) occurs via independent events of formation/dissociation of the pyrocarbonate molecule. This is similar to the Grotthuss mechanism in water, which involves the Zundel cation:

H 3 O + + H 2 O ⇌ H 5 O +
2 , with oxygen playing the role of the proton. By a random walk model, we estimate the diffusion coefficient of carbon dioxide to be D Grotthuss = 8.3 10 -9 m 2 s -1 , i.e. 2.8 times faster than molecular diffusion of CO 2 (see Supplementary Information for details). It would be of interest for future work to verify if in presence of an external electric field a cascade of C transfer is likely to occur, similarly to the correlated H transfers that have been observed in water under intense electric fields [START_REF] Saitta | Ab Initio Molecular Dynamics Study of Dissociation of Water under an Electric Field[END_REF].

Pyrocarbonate geometry

Having observed that the pyrocarbonate ion exists in a significant portion of our simulation run, we then proceed to characterise its geometry by measuring distances and angles, whose distributions are shown in Fig. 3. We plot in Fig. 3(a Supplementary Table 1 compares the mean values of the distances and angles in the pyrocarbonate molecule as found in our study with those obtained in silicate melts [START_REF] Vuilleumier | Carbon dioxide in silicate melts at upper mantle conditions: Insights from atomistic simulations[END_REF] and in the gas phase, as studied by MP2 calculations by Peeters et al. [START_REF] Peeters | Solubility and Solvation of Carbon Dioxide in the Molten Li 2 CO 3 /Na 2 CO 3 /K 2 CO 3 (43.5:31.5:25.0 mol-%) Eutectic Mixture at 973 K. II. Theoretical Part[END_REF]. The extracted values in the CaCO 3 solvent are quite similar to what has been found in silicate melts.

Therefore the geometry of the pyrocarbonate ions appears to stay approximately constant in ionic melts. In contrast, the central C-O * -C angle is much more bent in the condensed phase with respect to the gas phase, by approximately 20 • . Correspondingly, a smaller C-C distance is found in the condensed phase. We attribute the bending of the pyrocarbonate ion to charge screening occurring in the condensed phase: the charge in the pyrocarbonate ion is mainly carried by the external O atoms, while the central O * is mainly neutral, see for example the O-Ca RDFs shown in Supplementary Figure 2. Therefore, charge screening by the surrounding cations allows the charged molecular ends to get closer to each other, with respect to the gas phase. The O-C-O angle and consequently the O-C-O * angle instead, remain quite similar to the gas phase case. To conclude the geometrical analysis on the pyrocarbonate ion, we calculate the O-C-O * -C dihedral angle, see Supplementary Figure 3.

We find that the relative orientation in the molecule is rather disordered, probably due to the transients of formation of the pyrocarbonate molecule out of CO 2 and the carbonate ion. However a preference for periplanar configurations can be observed.

Solvation of molecular CO 2

Not only does CO 2 form C 2 O 2- 5 but even in its molecular state it exhibits specific interactions with CO 2- 3 . In order to characterise the solvation structure of molecular CO 2 , we build two-dimensional histograms of the structure around the CO 2 molecule (see Methods for details). In Fig. 4 We find that the probability of having at least one O in that region is roughly 90%. When we consider the central ovoidal red area, the probability of having more than three O, thus necessarily more than one carbonate, is close to 80%. This shows that the CO 2 molecule is closely approached by carbonate anions on a preferred pathway, eventually leading to the formation of pyrocarbonate -similarly to what has been observed for CO 2 in imidazolium acetates [START_REF] Kelemen | An Abnormal N-Heterocyclic Carbene-Carbon Dioxide Adduct from Imidazolium Acetate Ionic Liquids: The Importance of Basicity[END_REF]. In order to get a clearer idea of the relative orientation of CO 2 and carbonate in this precursor arrangement, we perform a gas phase geometry optimisation on a cluster composed of CO 2 , CO 2- and O, a dark spot is visible for Ca. However Ca does seem to play a role in the approach of carbonate to CO 2 . In fact, the position of Ca found in the gas phase (see Fig. 4(d))

corresponds to the intense region just outside the excluded volume. Small bright spots are also visible at ρ ≃ 0 and at distances (along ζ) 3-4 times the C-O distance (≃ 1.17 Å),

where carbonate ions are already interposed to screen the Ca 2+ charge.

Pyrocarbonate solvation

In contrast to the case of molecular CO 2 , the solvation structure of pyrocarbonate is more diffuse. Given the mostly periplanar symmetry indicated by the dihedral angle, see Supplementary Figure 3, we go on in neglecting one dimension and we produce two-dimensional histograms, in analogous way to what is shown in Fig. 4 for the CO 2 molecule. In the present case, the vertical, ζ axis coincides with the C-C direction, and ρ is the distance from the this axis. The centre of the coordinate system is taken at the mid-point of the C-C vector, see Methods for further details. We also keep track of the C atom being part of the CO 2 molecule before the formation of the pyrocarbonate and we make the charts in such a way that the C formerly belonging to CO 2 is always found at the bottom of the plots. Furthermore, we notice a small bright spot at ρ ≃ 0, ζ ≃ 4.5, probably corresponding to the fact that Ca is close to the previous CO 2- 3 unit, rather than to the previous CO 2 . Finally, in Fig. 5(d) we plot the three-dimensional isoprobability surfaces for the atoms composing the pyrocarbonate, derived from three-dimensional histograms calculated as described in Methods. The plot shows how the position of the O in pyrocarbonate is orientationally disordered in the condensed phase, probably due to the transients of formation/dissociation.

When the isoprobability value is increased distinct spots appear that are consistent with periplanar configurations, 0 • < ϕ < 30 • , or 150 • < ϕ < 180 • , and with configurations with the two pairs of external O nearly perpendicular to each other 80 • < ϕ < 110 • , see also Supplementary Figure 3.

Conclusions

In conclusion, our study of carbon dioxide solvation in molten calcium carbonate by first principle molecular dynamics confirms the presence of the pyrocarbonate anion and of the

equilibrium CO 2 + CO 2- 3 ⇌ C 2 O 2- 5 ⇌ CO 2- 3 + CO 2 .
We have observed how the transport of CO 2 in molten carbonate occurs in a similar fashion to the Grotthuss mechanism in water:

we name this mechanism oxo-Grotthuss mechanism. We have characterised the geometry of pyrocarbonate and found it has a periplanar configuration, with the central C-O * -C bridge much more bent than in the gas phase, due to charge screening by the counterions. Furthermore, we have observed that the interaction between CO 2 and CO 2- 3 leading to the formation of pyrocarbonate occurs along a preferred pathway with the carbonate ion approaching CO 2 "laterally". The enhanced transport of CO 2 by O 2-exchange -three times faster than molecular diffusion -may also explain the diffusivity of CO 2 in mixed oxide-carbonate materials [START_REF] Zhang | First spectroscopic identification of pyrocarbonate for high CO 2 flux membranes containing highly interconnected three dimensional ionic channels[END_REF].

The CO 2 /CO 2- 3 equilibrium can also be described in the framework of oxo-acidity, with CO 2 being the acid and CO 2- 3 its conjugated base. In this context, Grotthuss-like transport may be rather general when an acid is dissolved in its conjugated base.

METHODS

First-principles molecular dynamics simulations

We study a system originally composed of 127 CO 2- 3 carbonate anions, 127 Ca 2+ calcium cations and one CO 2 molecule (638 atoms, 4080 valence electrons). We fix the side of cubic simulation box to L = 20.054 Å, which gives a density of ρ = 2.63 g.cm -3 . The CO 2 concentration is thus 0.206 mol L -1 , a value typical for alkali molten salts under 1 atm of CO 2 [START_REF] Chery | Overview on CO 2 valorization: challenge of molten carbonates[END_REF]. We set the temperature, controlled by the Nosé-Hoover thermostat chain [START_REF] Nosé | A molecular-dynamics method for simulations in the canonical ensemble[END_REF][START_REF] Nosé | A unified formulation of the constant temperature molecular-dynamics methods[END_REF] with time constant τ = 0.5 ps, at T = 1773 K. This thermodynamic point is located, according to the pure CaCO 3 melting curve obtained by Suito et al. [START_REF] Suito | Phase relations of CaCO 3 at high pressure and high temperature[END_REF] and Spivak et al. [START_REF] Spivak | Stability and breakdown of Ca 13 CO 3 melt associated with formation of 13 C-diamond in static high pressure experiments up to 43 GPa and 3900[END_REF], along the melting curve of disordered pure calcite (no CO 2 ).

We perform FPMD simulations based on DFT and the Born-Oppenheimer dynamics.

We use the software CP2K [START_REF] Hutter | 2K: atomistic simulations of condensed matter systems[END_REF] and in particular the Quickstep algorithm [START_REF] Vandevondele | QUICKSTEP: Fast and accurate density functional calculations using a mixed Gaussian and plane waves[END_REF]. The latter employs a hybrid Gaussian plane-wave method (GPW) [START_REF] Lippert | A hybrid gaussian and plane wave density functional scheme[END_REF]. We use the Goedecker-Teter-Hutter (GTH) norm-conserving pseudo-potentials [START_REF] Goedecker | Separable dual-space Gaussian pseudopotentials[END_REF][START_REF] Hartwigsen | Relativistic separable dual-space Gaussian pseudopotentials from H to Rn[END_REF][START_REF] Krack | Pseudopotentials for H to Kr optimized for gradient-corrected exchange-correlation functionals[END_REF] to replace the core electrons. We cut off the electronic density at 400 Ry and use NN50 smoothing for applying the exchangecorrelation potential. For C and O atoms, we use a triple-zeta valence doubly polarised (TZV2P) basis set [START_REF] Vandevondele | The influence of temperature and density functional models in ab initio molecular dynamics simulation of liquid water[END_REF]. For Ca atoms, we use a double-zeta valence plus polarisation (DZVP) basis set optimised for molecules [START_REF] Vandevondele | Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases[END_REF]. The exchange-correlation interactions are taken into account by the gradient-corrected BLYP functional [START_REF] Becke | Density-functional exchange-energy approximation with correct asymptotic behavior[END_REF][START_REF] Lee | Development of the Colle-Salvetti correlation energy formula into a functional of the electron density[END_REF]. We add dispersive interactions corrections by using the DFT-D2 scheme [START_REF] Grimme | Semiempirical GGA-type density functional constructed with a long-range dispersion correction[END_REF], with a cutoff of 40 Å.

The starting configuration has been taken from previous equilibrated FPMD simulations of the pure liquid CaCO 3 (128 units) [START_REF] Vuilleumier | Structure, equation of state and transport properties of molten calcium carbonate (CaCO 3 ) by atomistic simulations[END_REF] after having replaced one CO 2- 3 with CO 2 and having erased one Ca 2+ from the simulation box. As mentioned earlier, we run the simulation in the canonical N V T ensemble, with the simulation time step set to 0.5 fs and periodic boundary conditions. We store the trajectory at every time step for post-run-time analysis.

We accumulate the trajectory for approximately 65 ps (130000 simulation steps).

In addition, we also perform a FPMD simulation of CO 2 in the gas phase, as well as geometry optimisations of the pyrocarbonate ion in the gas phase and a cluster composed of CO 2 , CO 2- 3 and Ca 2+ (shown in Fig. 4(d)). Details for these can be found in Supplementary Information.

Analysis

We perform the analysis of the simulation trajectory by using an in-house analysis suite.

In order to define molecules from the atomic positions at each time step, we use a depth first search (DFS) closure algorithm [START_REF] Cormen | Introduction to Algorithms[END_REF]. We use a particular "hysteresis" criterion in order to define whether in a given time frame, we have molecular CO 2 or pyrocarbonate ion. If we have molecular CO 2 in the previous time frame, we account for the formation of the pyrocarbonate in the current frame only when the C -O * distance is less than 1.7 Å. On the other hand, if we have pyrocarbonate in the previous frame, we account for a dissociation event in the current frame only when the C -O * distance becomes larger than 1.95 Å. This criterion allows smoothing out the behaviour that is observed when a simple distance cutoff criterion is used and it is inspired to the criteria used for correlations in H-bonds dynamics in water [START_REF] Laage | A Molecular Jump Mechanism of Water Reorientation[END_REF].

To build the two-dimensional colour charts around the CO 2 molecules (see Fig. For the construction of the two-dimensional colour charts around the pyrocarbonate ion (see Fig. 5), we follow an analogous procedure, with the only difference that now the C-C direction is taken as ζ and the mid-point of the C-C vector as centre. Furthermore, we specifically take the beginning of the C-C vector corresponding to the C atom that was part of the CO 2 molecule just before the formation of the pyrocarbanate ion. Therefore the former C of the CO 2 molecule is always found in the negative ζ portion of the plots.

We have also built three dimensional colour charts for C, O and Ca (see for example

= r ′ i • ξ, υ = r ′ i • υ and ζ = r ′ i • ζ.
We bin the space by ∆ξ = ∆υ = ∆ζ = 0.1 Å between -6.0 Å and 6.0 Å. We normalise the intensity histograms by dividing by the number of steps, the number of particles of type α (α = C, O, Ca), and by the volume of the bin ∆ξ∆υ∆ζ. Molecular representations and iso-probability surfaces are visualised using the software VMD [START_REF] Humphrey | VMD -Visual Molecular Dynamics[END_REF].

Figure 1 .

 1 Figure 1. Oxo-Grotthus mechanism via a pyrocarbonate anion. (a) Formation of the pyrocarbonate ion C 2 O 2- 5 out of the reaction of CO 2 with the carbonate anion CO 2- 3 . (b) A cascading mechanism, with rapid sequence of pyrocarbonate formation and dissociation events.

  criterion with minimal/maximal C -O * distances, where O * is the central pyrocarbonate ion linking the two C atoms, C -O * -C (the other oxygen atoms in pyrocarbonate are simply named O in the following), see Methods for further details. The minimal/maximal distances considered are 1.7/1.95 Å since within this span the C-O pair RDF is practically zero for molten carbonates.

  ) the distribution of the C-C, C-O and C-O * distances. We observe how the C-O distribution is quite narrow, peaking at about 1.26 Å, while the C-C and C-O * distance distributions are quite wide. The latter has an asymmetric shape, which is related to the criterion chosen to define the pyrocarbonate molecule, see Methods for more details. In Fig. 3(b), we show the normalised distributions of the C-O * -C, O-C-O and O-C-O * angles for pyrocarbonate. The O-C-O and O-C-O * distributions are narrower, while the C-O * -C shows a larger dispersion. This signals that the central C-O * -C bridge is remarkably flexible. For comparison, we also show the O-C-O angle distribution for CO 2- 3 molecules, which is narrow and peaks at about 120 • . In the inset, see Fig. 3(c), we also show the logarithm of the normalised distribution of O-C-O angle of the CO 2 molecule and we compare it with the case of CO 2 in the gas phase (see Methods for further details). We remark that the fluctuations of the O-C-O angle are significantly enhanced in the condensed phase with respect to the gas phase. In fact the O-C-O angle can even bend to ≃ 120 • . Despite the low polarisability of CO 2 , we see here that the interactions with the surrounding medium, in this case composed of doubly charged cations and anions, can influence deeply the geometry of CO 2 . In particular, the bending of the O-C-O can be related to its interaction with CO 2- 3 and the formation of the pyrocarbonate. Indeed in pyrocarbonate, we observe a O-C-O angle distribution peaking at about 130 • .

  (a), Fig. 4(b) and Fig. 4(c), we show these probability density maps for C, O and Ca atoms respectively. To facilitate the visualisation and improve statistics, we neglect here the instantaneous bending of the CO 2 molecule and we build the probability histograms in a cylindrical ρ -ζ coordinate system, where the ζ direction coincides with the O-O unit vector and ρ is the distance from the ζ axis. The zero of the coordinate system is taken at the mid-point of the O-O vector, i.e. the average position of the C atom. Looking at Fig. 4(a) and Fig. 4(b), we remark the existence of a preferred position for the solvation by CO 2- 3 in the central area of the ρ -ζ representation. As a matter of fact the largest intensities for the C atoms are found in the region 2.75 Å < ρ < 3.90 Å and -1.85 Å < ζ < 1.85 Å. For O atoms instead, we notice a bright spot in the area 2.25 Å < ρ < 3.00 Å and -0.55 Å < ζ < 0.55 Å. In Supplementary Figure 4 we show the histogram of the number of O atoms found in the brightest region evidenced by Fig. 4(b).

3 and

 3 Fig 4(d). We find that one O of carbonate points directly to the C of CO 2 with the C CO 2 -O carb. -C carb. angle being ≃ 108 • . The distances found in the gas phase are compatible with the two-dimensional colour chart of Fig. 4(b) although in the condensed phase the C CO 2 -O carb. -C carb. angle is closer to ≃ 90 • . In addition to this structural analysis, we present in Supplementary Figure 5 a view of the electronic state of the molecules during the formation of the pyrocarbonate, by looking at the centers of the maximally localized Wannier functions. In Fig. 4(d) we see that a large excluded volume region is present for Ca atoms around the CO 2 molecule. Interestingly, in the region corresponding to the brightest spot for C

Fig. 5 (

 5 Fig. 5(a) shows the two-dimensional histograms for C atoms around the pyrocarbonate. We observe that the solvation of pyrocarbonate by C atoms is almost uniform with a small dent in the central region around ζ ≃ 0 and ρ ≃ 3 corresponding to a modest electrostatic attraction by the central O * atom of pyrocarbonate and the repulsion by the two C atoms.For the rest, we observe three consecutive solvation shells, with the first one being the most

  4), we first recenter the coordinates, with the centre r c defined as the mid-point of the CO 2 O-O vector. We define the ζ direction as the O-O unit vector. We then calculate the coordinate ζ of a given atom as ζ = r ′ i • ζ, where r ′ i = r i -r c . The coordinate ρ is given by ρ = √ ||r ′ i || 2 -ζ 2 or equivalently by ρ = ||r ′ i -ζ ζ||. We bin the space by ∆ζ = ∆ρ = 0.1 Å in the intervals between -L/2 and L/2 for ζ and 0 and √ 2L/2 for ρ. We normalise the two-dimensional intensity histogram dividing by the product of the number of steps, the number of particles of type α, with α = C, O, Ca, and the bin volume 2πρ∆ρ∆ζ.

Fig. 5 (

 5 Fig. 5(d)). In those plots, the axis ζ coincides with the C-C direction, the υ axis is obtained from CO * -(CO * • ζ) ζ and ξ = υ × ζ. We recenter the coordinates around the mid-point of the C-C vector r c , r ′ i = r i -r c . The new coordinates are thus ξ = r ′ i • ξ, υ = r ′ i • υ

Figure 2 .

 2 Figure 2. Formation/dissociation of pyrocarbonate and C transport. The snapshots in (a), (b), (c) represent one instance of the formation/dissociation of pyrocarbonate during the simulation run. In (a) CO 2 and CO 2- 3 approach, in (b) pyrocarbonate C 2 O 2- 5 is formed, and in (c) pyrocarbonate dissociates into CO 2- 3 and CO 2 . Carbonate molecules in the solvent are represented by cyan/red sticks and red spheres for O atoms, Ca atoms are represented as green spheres. The CO 2 molecule is shown in yellow, the CO 2- 3 molecule in blue and pyrocarbonate in purple. Panel (d) shows the indices of the C atoms in CO 2 (black) and pyrocarbonate (red) along the simulation.In panel (e) we follow the C atom in CO 2 before the formation and after the dissociation of pyrocarbonate; we plot in black the probability of that carbon atom staying n consecutive times the same (always reforming the same CO 2 molecule), and in red that of changing n consecutive times (cascading C transfers). The lines are as described in the text.

Figure 3 .

 3 Figure 3. Geometry of the pyrocarbonate anion. (a) Distribution of the C-C, C-O and C-O * distances measured during the time intervals of existence of C 2 O 2- 5 . Inset: average geometry of pyrocarbonate (O: red; C: grey), with values of the C-O* and C-O distances and C-O*-C and O-C-O angles shown. (b) Distribution of the C-O * -C, O-C-O, O-C-O * for the pyrocarbonate ion. For comparison, we also plot the distribution of the O-C-O angle in CO 2- 3 . We normalise the distributions by dividing by sin θ. (c) Natural logarithm of the distribution of the O-C-O angle in CO 2 during the intervals in which it is present in the CaCO 3 melt (green diamonds) and we compare it to its behaviour in the gas phase (black circles), see Methods. We also perform a parabolic fit of the form A(θ -180) 2 + B. The fits are shown as blue and red solid lines for the CaCO 3 and gas phase cases respectively. The values of the parameters extracted are: A = -0.0025, B = -1.25 (CaCO 3 ) and A = -0.0059, B = -0.42 (gas).

Figure 4 .Figure 5 .

 45 Figure 4. Solvation structure around the CO 2 molecule. We build two-dimensional colour charts (see Methods) around the CO 2 molecule for (a) C atoms, (b) O atoms and (c) Ca atoms. In the plots the CO 2 molecule is oriented with its O-O axis along the ζ direction, with ζ = 0 being the mid-point of the O-O distance. Therefore ζ is the projection of an atom position along the O-O axis, while ρ is its distance from the O-O axis. (d) Gas phase geometry obtained for a cluster composed of CO 2 , CO 2- 3 and Ca 2+ (see Methods). The Cartesian axes are shown at the bottom, C is represented in grey, O in red and Ca in yellow. The positions of the atoms of carbonate and Ca in the gas phase are marked as + symbols in the respective plots for the condensed phase. Distances are measured in Å.

has been observed in the basaltic melt, accounting for ≃ 3.7% of the C content in the system. Pyrocarbonate has been found to be a transient species with a very short mean life time of 175 fs, too short to assess its solvation structure and formation/dissociation dynamics.In this study, we focus on CO 2 in molten CaCO 3 , with the aim of characterising the formation/dissociation of the pyrocarbonate and the solvation structure of CO 2 and pyrocarbonate in the melt. Using FPMD simulations, we confirm the presence of pyrocarbonate which exists in rapid exchange with separate CO 2 and carbonate. The pyrocarbonate ion lives longer in CaCO 3 than in the basaltic melt reported in Ref.[START_REF] Vuilleumier | Carbon dioxide in silicate melts at upper mantle conditions: Insights from atomistic simulations[END_REF], allowing us to study the details of its solvation structure, beyond its geometry. In addition, we find that the transport of CO 2 in molten carbonate occurs in a similar fashion to the Grotthuss mechanism for proton transport in water, via independent events of formation and dissociation of the pyrocarbonate molecule.

ACKNOWLEDGEMENTS

We thank Michel Cassir, Virginie Lair, Bertrand Guillot, Fabrice Gaillard, Volker Haigis and Anne Boutin for fruitful discussions. The research reported herein was funded by PSL Research University (project COOCAR, grant ANR-10-IDEX-0001-02) and Agence Nationale pour la Recherche (project ELECTROLITH, grant ANR-2010-BLAN-621-03). This work has been performed using HPC resources from GENCI (grants 2013-082309, 2014-082309 and 2015-082309) and IDRIS (grant "Grand Challenge" 100577). We also acknowledge PRACE for awarding us access to resource Curie based in France at CEA Bruyères-le-Chatel (preparatory access allocation 2010PA2746).

AUTHOR CONTRIBUTIONS

R. V. performed the FPMD simulations, D. C. analysed the trajectories, prepared the figures and wrote the manuscript. All authors designed research, discussed the results and revised the manuscript.

ADDITIONAL INFORMATION

Correspondence and requests for materials should be addressed to F.-X.C. or R.V.

Competing financial interests

The authors declare no competing financial interests.

⇆

Table of contents graphic

Table of contents summary:

First principles molecular dynamics simulations have shown that solvation of carbon dioxide in molten calcium carbonate leads to the formation of the pyrocarbonate anion, C 2 O 2- 5 . Despite its short life time, the formation of the pyrocarbonate anion significantly impacts the transport of CO 2 , which occurs in a fashion similar to the Grotthuss mechanism in water.