Lianke Wang

Marcello Gennari

Alexandre Barrozo

Jennifer Fize

Christian Philouze

Serhiy Demeshko

Franc Meyer

Maylis Orio

Vincent Artero

Carole Duboc

Role of the Metal Ion in Bio-Inspired Hydrogenase Models: Investigation of a Homodinuclear FeFe Complex vs Its Heterodinuclear NiFe Analogue

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

Introduction

Nature designed two distinct dinuclear catalysts to reversibly reduce protons into dihydrogen. Two metal ions act in synergy in such {MFeS2} cores (M = Ni or Fe, Scheme 1), present in the active site of the [NiFe] and [FeFe] hydrogenases (H2ases), respectively. [1][START_REF] Schilter | Hydrogenase Enzymes and Their Synthetic Models: The Role of Metal Hydrides[END_REF][START_REF] Nicolet | Desulfovibrio desulfuricans iron hydrogenase: the structure shows unusual coordination to an active site Fe binuclear center[END_REF][START_REF] Chen | Infrared Studies of the CO-Inhibited Form of the Fe-Only Hydrogenase from Clostridium pasteurianum I: Examination of Its Light Sensitivity at Cryogenic Temperatures[END_REF][START_REF] Peters | X-ray Crystal Structure of the Fe-Only Hydrogenase (CpI) from Clostridium pasteurianum to 1.8 Angstrom Resolution[END_REF][START_REF] Volbeda | Crystal structure of the nickel-iron hydrogenase from Desulfovibrio gigas[END_REF] Although both enzymes involve a similar ECEC mechanism [START_REF] Coutard | Molecular engineered nanomaterials for catalytic hydrogen evolution and oxidation[END_REF] to produce H2 (E = electron transfer, C = chemical reaction, here corresponding to a proton transfer), different types of intermediates and strategies are involved. 1,[START_REF] Can | Structure, Function, and Mechanism of the Nickel Metalloenzymes, CO Dehydrogenase, and Acetyl-CoA Synthase[END_REF] For the [NiFe] H2ase, only the Ni ion participates in the redox chemistry by cycling between oxidation states ranging from +I to +III. The Fe site remains in the +II oxidation states but contributes in the stabilization of hydride ligands, bridging both metal centers in two key intermediates in the catalytic cycle. [START_REF] Ash | Proton Transfer in the Catalytic Cycle of [NiFe] Hydrogenases: Insight from Vibrational Spectroscopy[END_REF][START_REF] Ogata | Structure and function of [NiFe] hydrogenases[END_REF][START_REF] Hidalgo | Infrared Spectroscopy During Electrocatalytic Turnover Reveals the Ni-L Active Site State During H2 Oxidation by a NiFe Hydrogenase[END_REF][START_REF] Ogata | Hydrogens detected by subatomic resolution protein crystallography in a [lsqb]NiFe[rsqb] hydrogenase[END_REF] Regarding [FeFe] H2ase, the redox chemistry also occurs at a single metal site (the distal iron center, Fed). It only cycles between two oxidation states (+I and +II), and is involved in the stabilization of a terminal hydride ligand. The second iron center (proximal Fe site, Fep) acts as a scaffold through the positioning of the doubly-bridging azadithiolate ligand that plays a key role as proton relay. Besides it also connects Fed with a [4Fe4S] cluster that acts as an electron reservoir and cycles between +1 and +2 charge states during the catalytic cycle. [START_REF] Tard | Structural and Functional Analogues of the Active Sites of the [Fe]-, [NiFe]-, and [FeFe]-Hydrogenases[END_REF][START_REF] Adamska | Identification and Characterization of the "Super-Reduced" State of the H-Cluster in FeFe Hydrogenase: A New Building Block for the Catalytic Cycle?[END_REF][START_REF] Mulder | EPR and FTIR Analysis of the Mechanism of H2 Activation by [FeFe]-Hydrogenase HydA1 from Chlamydomonas reinhardtii[END_REF] In this context, chemists applied two main approaches to mimic the active sites of these two enzymes. [START_REF] Schilter | Hydrogenase Enzymes and Their Synthetic Models: The Role of Metal Hydrides[END_REF][START_REF] Simmons | Mimicking hydrogenases: From biomimetics to artificial enzymes[END_REF][START_REF] Gloaguen | Small molecule mimics of hydrogenases: hydrides and redox[END_REF] Bio-inspired NiFe catalysts for H2 production have been mainly designed with a {NiFeS2} core consisting of a Ni site in a thiolate-rich environment and a Fe site with CO, CN -, Cp -or carbene ligands. The main objective was to localize the redox chemistry at the Ni site. [START_REF] Kaur-Ghumaan | NiFe] hydrogenases: how close do structural and functional mimics approach the active site?[END_REF] In the case of the FeFe mimics, even if the majority also displays a {Fe2S2} core, a special interest has been focused on the bridging azadithiolate ligand that bears an amine proton relay. When possible, the incorporation of an electron transfer mediator in the second coordination sphere, such as a [4Fe4S] cluster [START_REF] Tard | Synthesis of the H-cluster framework of iron-only hydrogenase[END_REF] or a ferrocenyl group, [START_REF] Camara | Combining acid-base, redox and substrate binding functionalities to give a complete model for the FeFe -hydrogenase[END_REF] has also been considered. [START_REF] Rauchfuss | Diiron Azadithiolates as Models for the [FeFe]-Hydrogenase Active Site and Paradigm for the Role of the Second Coordination Sphere[END_REF] Due to these specific requirements, a direct comparison between the hydrogen evolution reaction (HER) activity of NiFe and FeFe H2ase mimics with the same coordination environment has not been reported so far. Two series of NiFe and FeFe complexes with similar environments have been described but the complexes are not isostructural: [(dppe)Ni(pdt)Fe(CO)3]/[(dppe)(CO)Fe(-pdt)Fe(CO)3] 17, [START_REF] Barton | Nickel-Iron Dithiolato Hydrides Relevant to the [NiFe]-Hydrogenase Active Site[END_REF][START_REF] Zhu | Modulation of the electronic structure and the Ni-Fe distance in heterobimetallic models for the active site in [NiFe]hydrogenase[END_REF] and [(L' N2S2 )NiFe(CO)Cp] + / [(L' N2S2 )(Fe-NO)Fe(CO)Cp] + (L' N2S2 = bismercaptoethyl diazacycloheptane) [START_REF] Ding | Hemilabile Bridging Thiolates as Proton Shuttles in Bioinspired H2 Production Electrocatalysts[END_REF] (Scheme 1). Regarding the first series, even if both complexes can generate similar hydride intermediates, only the NiFe derivate catalyses H2 production. As for the second series, a {Fe(NO)} 7 is required to stabilize the FeFe complex in such an environment. Because of the peculiar redox character of the {Fe(NO)} n unit, it is difficult to conclude on the effect of the nature of the metal, Ni vs Fe, on the catalytic process in this specific case.

We have recently investigated the propensity of a NiFe complex, [L N2S2 Ni II Fe II Cp(CO)] + ( , with L N2S2 = 2,2'-(2,2'-bipyridine-6,6'-diyl)bis(1,1'-diphenylethanethiolate), Fe' = in the {FeCp moiety), to mimic the structural and functional properties of [NiFe] H2ase (Scheme 1). This efficient electrocatalyst for H2 production displays redox chemistry mainly centered on the Ni-N2S2 site, and two catalytic intermediates, including a hydride species, have been identified in homogeneous conditions. [START_REF] Brazzolotto | Nickel-centred proton reduction catalysis in a model of [NiFe] hydrogenase[END_REF][START_REF] Brazzolotto | Tuning Reactivity of Bioinspired NiFe -Hydrogenase Models by Ligand Design and Modeling the CO Inhibition Process[END_REF][START_REF] Tang | Biomimetics of [NiFe]-Hydrogenase: Nickel-or Iron-Centered Proton Reduction Catalysis?[END_REF][START_REF] Ahmed | Hydrogen Evolution from Aqueous Solutions Mediated by a Heterogenized NiFe -Hydrogenase Model: Low pH Enables Catalysis through an Enzyme-Relevant Mechanism[END_REF] In this work, we report the structure, properties and HER reactivity in homogeneous conditions of the parent FeFe complex, (Scheme 1), with the aim of comparing it with its NiFe counterpart and of understanding the role of the metal in the N2S2 site, Ni vs Fe.

Scheme 1.

Active sites of the [FeFe] and [NiFe] H2ases (resting state), selected examples of mimics discussed in the text.

Results and Discussion

Synthesis and solid-state characterization of . The tetrafluoroborate salt of the diiron complex [(CO)L N2S2 Fe II Fe II (CO)Cp] + (with L N2S2 = 2,2′-(2,2′-bipyridine-6,6′diyl)bis(1,1-diphenylethanethiolate and Cp = cyclopentadienyl), , was isolated as a brown solid by reacting the [CpFe(CO)(MeCN)2]BF4 precursor with the previously reported dinuclear Fe II complex [Fe II 2(L N2S2 )(HL N2S2 )]BF4 ([Fe2 SH ]BF4) [START_REF] Wang | A Non-Heme Diiron Complex for (Electro)catalytic Reduction of Dioxygen: Tuning the Selectivity through Electron Delivery[END_REF] in MeCN at room temperature (Scheme 2). The formation of is the result of the cleavage of the [Fe2 SH ] + complex making two thiolate-donor ligands available for binding the [CpFe(CO)] + moiety. Two iron-bound CO molecules are present in the structure of (see below), so that the use of an excess of [CpFe(CO)] + (4 equiv. vs [Fe2 SH ] + ) as CO source is required for its synthesis. The diffusion of Et2O into an MeCN solution of in the presence of tetrabutylammonium perchlorate provides black needle crystals suitable for single-crystal X-ray crystallographic analysis. Scheme 2. Synthesis of the diiron complexes discussed in this work. The notations Fe and Fe' refer to the {FeL N2S2 } and {FeCp} moieties, respectively. The X-ray structure of the perchlorate salt of (Figure 1, left, and Tables S1-S2) reveals a {Fe2S2} core in which two inequivalent iron(II) centers are thiolatebridged. The Fe1 ion (corresponding to the Fe unit in

) is coordinated by one terminally bound CO molecule, two sulfur and two nitrogen atoms of L N2S2 , resulting in a distorted square pyramidal geometry with a τ5 value of 0.12. The Fe2 ion (corresponding to the Fe' unit in

) is surrounded by a CO molecule, two bridging sulfur atoms of L N2S2 and a  5 -coordinated Cp ligand. No Fe … Fe interaction is present, in agreement with the observed Fe … Fe distance of 2.5611(12) Å. The latter is similar to the Fe … Fe distances found in [FeFe] hydrogenases in the resting Hox state (2.55-2.62 Å). [START_REF] Peters | X-ray Crystal Structure of the Fe-Only Hydrogenase (CpI) from Clostridium pasteurianum to 1.8 Angstrom Resolution[END_REF][START_REF] Mulder | EPR and FTIR Analysis of the Mechanism of H2 Activation by [FeFe]-Hydrogenase HydA1 from Chlamydomonas reinhardtii[END_REF] The Fe1S1S2 plane is almost perpendicular to the Fe2S1S2 plane with a dihedral angle of 87.93(6) °. The CO ligand coordinated to Fe2 is semi-bridging between the two iron centers (Fe1 The powder IR spectrum of (Figure 2, left) displays two CO stretching vibration bands at 1959 and 1918 cm -1 , assigned to the terminal and semi-bridging CO ligands, respectively (Table 1). The powder zero-field 57 Fe Mössbauer spectrum of recorded at 80 K (Figure 1, right) exhibits two quadrupole doublets with δ = 0.23 mm . s -1 ; ∆EQ = 0.69 mm . s -1 for FeA and δ = 0.35 mm . s -1 ; ∆EQ = 2.32 mm . s -1 for FeB, consistent with two distinct low spin iron(II) centers (Table 1). Based on the Mossbauer parameters of the Fe II center from the previously reported parent complex (δ = 0.39 mm . s -1 and ∆EQ = 1.82 mm . s -1 ), [START_REF] Brazzolotto | Nickel-centred proton reduction catalysis in a model of [NiFe] hydrogenase[END_REF] FeB can be assigned to the Fe II ion located in the Fe' subunit. Consequently, FeA corresponds to the Fe II ion in the Fe subunit. The solid state UV-vis spectrum of (Figure 2) displays two absorption bands at around 381 nm and 483 nm with a broad shoulder around 650 nm. The EPR spectrum of is silent, in agreement with the low spin character of both Fe II ions. a the corresponding spectrum is displayed in the SI. b the NIR vibration corresponds to an intervalence charge transfer transition (see Figure 7).

Spectroscopic characterization of in solution.

The IR spectrum of dissolved in CH2Cl2 displays CO stretching vibration bands at 1971 cm -1 and 1918 cm -1 , similar to those in powder spectrum (Figure 2, left). The solid-state structure is thus retained in CH2Cl2 solution. On the other hand, the IR spectrum of dissolved in MeCN is more complex to interpret. In a concentrated MeCN solution (1.0 mM), the two νCO bands are still present, but with an additional feature at 1822 cm -1 , corresponding to a bridging CO. Conversely, in a diluted MeCN solution (0.2 mM), the νCO band at 1822 cm -1 is almost exclusively observed. The disappearance of the bands associated with the terminal and semi-bridging CO ligands of together with the concomitant appearance of a bridging CO ligand vibration can be attributed to the formation of a new species, named , in which one of the CO has been released, and the other has been shifted towards a bridging position. Based on the fact that this reactivity is observed in MeCN but not in CH2Cl2, we believe that MeCN is involved in the decoordination of the CO ligand, either indirectly (polarity effect) or by direct coordination through exchange with one CO. Scheme 2 displays this proposed equilibrium reaction. Accordingly, in the 1 H NMR spectrum of in CD2Cl2 only one set of signals is observed, in accordance with the molecular structure observed in solid state, while two set of signals are present in CD3CN (1.0 mM solution), corresponding to and (Figure 3). The UV-vis signature of in CH2Cl2 and MeCN solutions (Figure 2, right) confirms the speciation determined by IR and 1 H NMR. The spectrum in CH2Cl2 (main bands: 382 nm, 4110 M -1 •cm -1 and 475 nm, 1800 M -1 •cm -1 ) is similar to the solid-state spectrum, confirming that the structure of is retained in CH2Cl2 solution. In MeCN, at high concentration (1.0 mM), the three absorption bands at 381 nm, 429 nm and 536 nm correspond to a mixture of both and , whereas at low concentration (0.2 mM), the spectrum shows only the absorption features of , i.e. two bands at 434 nm (2989 M -1 •cm -1 ) and 536 nm (2024 M -1 •cm -1 ). The equilibrium shown in Scheme 2 can be directed towards the formation of either or by tuning the experimental conditions. When CO is bubbled into the mixture of both species in MeCN (0.5 mM ), the band at 381 nm increases, while those at 429 and 536 nm disappear with the concomitant appearance of an additional absorption band at 472 nm (Figure 4, left). The resulting UV-vis spectrum is similar to that of in CH2Cl2, indicating that CO bubbling leads to the quantitative formation of the complex also in MeCN. When CO is removed by argon bubbling from the CO-saturated MeCN solution of (Figure 4, right), the two bands at 381 and 472 nm disappear and the two bands characteristic of newly appear. In both forward and backward processes, an isosbestic point at 415 nm is observed, demonstrating that only the two and species are involved in the equilibrium shown in Scheme 2. To further investigate the structural and spectroscopic properties of , DFT calculations have been performed. We initially validated our approach through the study of . Its optimized structure, *, displays structural parameters (Table 1) and predicted CO vibrations that are fully consistent with the experimental data (less than 8 pm of difference between predicted and experimental metal-ligand bond lengths; CO calc at 1900 and 1960 cm -1 vs CO exp at 1918 and 1959 cm -1 ). The Mössbauer parameters have been also predicted by DFT and the calculated parameters are well consistent the experimental data confirming our previous attribution: for the Fe unit, δcalc = 0.26 mm . s -1 ; ∆EQcalc = 0.82 mm . s -1 and for the Fe' unit, δcalc = 0.39 mm . s -1 ; ∆EQcalc = 2.37 mm . s -1 (Table 1).

For , several structures containing only one metal-bound CO have been optimized with or without an iron-bound MeCN molecule. Two calculated structures are consistent with the experimental CO vibration at 1822 cm -1 : (i) one with a bridging CO, CO calc = 1801 cm -1 and no coordinated MeCN ( *) and (ii) one with a bridging CO and a MeCN molecule coordinated to the Fe center, CO calc = 1838 cm -1 ( ( * ) ). Considering the acetonitrile binding energy, the calculated free energy difference (1.9 kcal.mol -1 ) between these two species is consistent with the presence of both complexes in equilibrium (for the sake of simplicity it will be referred to in the rest of the manuscript).

Electrochemical properties of .

The electrochemical properties of the diiron(II) complex have been examined by cyclic voltammetry (CV) in MeCN at 0.2 mM concentration: under these conditions only , and not , is present as shown by IR. The CV (Figure 5) exhibits two diffusion-controlled reversible one-electron reduction waves at E1/2 = -1.21 V (∆E = 91 mV) and E1/2 = -1.65 V (∆E = 90 mV) versus Fc +/0 . The first cathodic process is assigned to the reduction of into the mixedvalence complex (Scheme 3), which has been isolated and fully characterized spectroscopically (see below). A similar cathodic process was previously observed for at a close potential (-1.29 V), but with the reduction occurring on the Ni site (Ni II to Ni I ). This indicates that the nature of the metal does not notably affect the redox properties of this class of dinuclear complexes. The second one-electron reduction process is proposed to mainly occur at the bipyridine moiety of the L N2S2 ligand, as shown in the case of (for which E1/2 = -1.90 V). [START_REF] Brazzolotto | Nickel-centred proton reduction catalysis in a model of [NiFe] hydrogenase[END_REF] The bipyridine is reduced at a noticeably less negative potential in the case of the FeFe compound, presumably because of structural and electronic differences in the oneelectron reduced species (see below). Even if the CV displays a reversible redox system for the second reduction process, the proposed radical-based intermediate is too unstable to be spectroscopically detected, as for the NiFe complex.ref nat chem The CV of also presents an irreversible anodic signal at 0.10 V (SI), assigned to the Fe II to Fe III oxidation in the Fe unit. The corresponding cathodic peak at -0.62 V is located at the same reduction potential of the previously reported mononuclear [L N2S2 Fe III MeCN] + complex. [START_REF] Wang | Solvent-and Halide-Induced (Inter)conversion between Iron(II)-Disulfide and Iron(III)-Thiolate Complexes[END_REF] This implies that the electrochemically oxidized form of is not stable on the CV timescale and undergoes a fast dissociation into [L N2S2 Fe III MeCN] + and [CpFe(CO)(MeCN)2] + . Such oxidation process was not observed in the case of .

Figure 5. CVs of (0.2 mM) in the absence (black) and presence of various amounts of Et3NHBF4 (red: 1 equiv.; blue: 2 equiv.; green: 5 equiv.; magenta: 10 equiv.; orange: 15 equiv.; cyan: 20 equiv.; dark yellow: 30 equiv.) in MeCN. Glassy carbon working electrode, 100 mV•s -1 , 0.1 M n-Bu4NClO4. Inset: CV of Et3NHBF4 (0.2 mM) in MeCN in the absence (black) and presence (red) of one equivalent of .

Spectroscopic characterization of the one-electron reduced species . On the basis of the electrochemical properties of , cobaltocene (CoCp2, E1/2 = -1.33 V) was used to chemically reduce it to generate the mixed-valence . This compound could be isolated as a powder and characterized (Table 1). The νCO band at 1896 cm -1 in the IR spectrum suggests the presence of one CO ligand terminally bound to one iron center. We ruled out the coordination of an MeCN ligand because of the absence of a CN vibration in the infrared spectrum (see SI). The zero-field Mössbauer spectrum (powder, 80 K, Figure 6, left) of exhibits two main quadrupole doublets, with FeA: δ = 0.51 mm . s -1 , ∆EQ =2.19 mm.s -1 (35 %) and FeB: δ = 0.33 mm.s -1 , ∆EQ =0.44 mm.s -1 (41 %), which can be assigned to the dinuclear mixed valence complex. The attribution of these parameters to a specific Fe site is arduous and has been made with a combined DFT analysis (see below). The sample also contains two additional minor Fe-based species probably arising from the decomposition of , characterized by two doublets δ3(∆EQ3) =0.10 (2.16) mm . s -1 (17%) and δ4(∆EQ4) = 1.17 (3.29) mm . s -1 (7%). The former displays parameters close to those of the [CpFe(CO)2CH2SiMe3] complex (δ(∆EQ) =0.099 (1.75) mm . s -1 ),3 and could then be attributed to a mononuclear [FeCOCp(solvent)2] + complex. The latter parameters are typical of a high spin Fe II species and comparable to those found for [Fe II 2(L N2S2 )(HL N2S2 )]BF4 (δ(∆EQ) 0.87 (3.76) mm . s -1 indicating that this species could correspond to a mononuclear Fe complex containing the L N2S2 ligand. The powder EPR spectrum of recorded at 22 K exhibits an axial signal with gǁ = 2.262 and g⏊= 1.947, corresponding to a localized S = 1/2 spin system (Figure 6, right). In contrast, the EPR spectrum recorded on the powder at 100 K features a broader and more isotropic signal (gi = 2.064), consistent with a delocalized mixed-valence species. Based on this temperature behavior, the electronic structure of can be described as a class II delocalized mixed-valence Fe I Fe II species. The identity and valence delocalization of has been confirmed also in solution both by EPR and Near Infrared (NIR) spectroscopy. The EPR spectra of in THF/MeCN recorded at 7 K and 100 K are similar to the corresponding powder spectra (Figure 6). The NIR spectrum displays a band centered at 1015 nm (9852 cm -1 ), 1/2 = 2280 cm -1 in agreement with a mixed-valence compound (see below). To gain more insights on the structural and electronic properties of , different possible structures have been optimized by DFT. Only one computed structure, *, which displays a terminally bound CO to Fe, is consistent with the experimental IR data (CO calc = 1904 cm -1 , νCO exp = 1896 cm -1 ). All spectroscopic parameters have been then predicted theoretically for * (Table 1). TD-DFT calculations predict a NIR transition at 1017 nm (f = 0.025) that perfectly matches the experimental band. This transition is characterized by a strong metal-based character consistent with an intervalence charge transfer (IVCT) transition (Figure 7). The Mössbauer parameters predicted for * match well with the experimental data, and can thus be assigned to each iron center: for the Fe unit, δcalc(exp) = 0.27(0.33) mm . s -1 ; ∆EQcalc(exp) = 0.64(0.44) mm . s -1 and for the Fe' unit, δcalc(exp) = 0.58(0.51) mm . s -1 ; ∆EQcalc(exp) = 1.91(2.19) mm . s -1 . The Mülliken spin population is consistent with a valence delocalized species with 0.40 and 0.40 for Fe and Fe', respectively (0.03, each N and 0.14, remaining). Finally the g isotropic value obtained from DFT is 2.088, with a rather small g anisotropy (g1 = 1.991; g2 = 2.124; g3 = 2.149) compared to the experimental axial EPR signal recorded at low temperature (g = 0.31). The combination of all these experimental data and theoretical information are consistent with a class II delocalized mixed-valence species with a structure in which the terminal CO is bound the Fe site and not bridging the two metal centers as in the case of the parent complex. Transition density for 1017 nm band : IVCT character as metal centers exclusively involved color code is yellow for donor state and red for acceptor state The terminal coordination of CO to the [FeL N2S2 ] moiety in (*) is consistent with the redox properties of the / system: although this reduction is proposed to mainly occur at the bipyridine site of L N2S2 as in , its potential is significantly higher (E1/2 = -1.65 V for vs E1/2 = -1.90 V for ). The fact that is easier to reduce than (i.e. the one-electron reduced form of ) may be due to the stronger back-bonding of the terminal CO ligand in vs the bridging CO in .

EPR par ameters of [COax

M etrical par ameters

Fe1-Fe2 Fe1-N1 Fe1-N2 Fe1-S1 Fe1-S2 Fe1-C1 Fe2-S1 Fe2-S2
Electrocatalytic activity for H2 production of . The electrocatalytic behavior of for H2 production has been evaluated using Et3NHBF4 as a mild proton source in MeCN (pKa = 18.6), [START_REF] Izutsu | Acid-Base Dissociation Constants in Dipolar Aprotic Solvents[END_REF] conditions that are identical to those used during the investigation of .(Natchem) In the presence of increasing amounts of Et3NHBF4 (0.2-6 mM), a catalytic process (Ecat/2 in the range -1.89 to -1.94 V for a concentration in Et3NHBF4 in the range 1 to 6 mM Et3NHBF4, vs Eno catalyst = -2.16 V vs Fc +/0 , potentials being measured at the half-wave) develops at potential ~250 mV more cathodic than the second process generating (Figure 5). The lower value of Ecat/2 compared to E1/2 / (see the inset of Figure 5) indicates that the two-electron reduced species reacts with a proton to generate an intermediate that needs to be further reduced at a lower potential before to react with a second proton to generate H2 (see below). The plots of catalytic current intensity (icat) versus the concentration of added Et3NHBF4 provide a linear correlation indicating a diffusion controlled regime for catalysis under these conditions (SI). [START_REF] Peters | X-ray Crystal Structure of the Fe-Only Hydrogenase (CpI) from Clostridium pasteurianum to 1.8 Angstrom Resolution[END_REF] Bulk electrolysis experiments performed at ~ -1.8 V vs Fc +/0 on a Hg-pool cathode (SI) confirmed the electrocatalytic production of H2 (detected by GC) in the presence of (~15 turnovers achieved with ~60% Faradaic yield within ~20 min, ~60% conversion from 10 mM Et3NHBF4).

Mechanistic considerations. The parent

electrocatalyst was shown to follow an E(ECEC) mechanism for H2 production. [START_REF] Brazzolotto | Nickel-centred proton reduction catalysis in a model of [NiFe] hydrogenase[END_REF][START_REF] Tang | Biomimetics of [NiFe]-Hydrogenase: Nickel-or Iron-Centered Proton Reduction Catalysis?[END_REF] The complex undergoes a oneelectron reduction at the Ni center to afford [L N2S2 Ni I Fe II Cp(CO)] ( , which is further reduced to produce the radical based [L N2S2 °Ni I Fe II Cp(CO)] -complex that is then protonated on Fe affording a semi-bridging hydride [L N2S2 Ni II (H)Fe II Cp(CO)] -. [L N2S2 Ni II (H)Fe II Cp(CO)] -. The complex then needs to be further reduced before to react with a second proton to produce H2. Based on our present experimental results and on the mechanism of the parent catalyst, [START_REF] Brazzolotto | Nickel-centred proton reduction catalysis in a model of [NiFe] hydrogenase[END_REF][START_REF] Brazzolotto | Tuning Reactivity of Bioinspired NiFe -Hydrogenase Models by Ligand Design and Modeling the CO Inhibition Process[END_REF][START_REF] Tang | Biomimetics of [NiFe]-Hydrogenase: Nickel-or Iron-Centered Proton Reduction Catalysis?[END_REF] a similar E(ECEC) mechanism is proposed for electrocatalytic H2 production with (Scheme 3).

Scheme 3.

Proposed catalytic pathways for H2 evolution with catalyst (potentials are referred to Fc +/0 ).

The first reduction generates the spectroscopically characterized mixed-valence complex (see above). This species undergoes a second one-electron reduction to form . The fact that the catalytic wave occurs at a reduction potential lower than E1/2 = -1.65 V versus Fc +/0 (corresponding to the / redox system) evidences that: (i) is protonated to generate either a Fe-hydride or thiol-containing intermediate, and (ii) this protonated species needs to be reduced before reacting with a second proton to produce H2.

is thus regenerated for another electrocatalytic cycle.

Interestingly, for the FeFe electrocatalyst the E(ECEC) cycle can be experimentally validated, while for the parent NiFe complex it was mainly predicted by DFT calculations. In the case of the predicted potential of the third reduction process (that occurs after the first protonation step) is higher than the second one (E 0 calc = -1.42 V and E 0 calc= -1.77 V vs Fc +/0 , respectively) [START_REF] Tang | Biomimetics of [NiFe]-Hydrogenase: Nickel-or Iron-Centered Proton Reduction Catalysis?[END_REF] and thus was not detectable. [START_REF] Brazzolotto | Nickel-centred proton reduction catalysis in a model of [NiFe] hydrogenase[END_REF] In contrast, for the third reduction process can be directly observed since it occurs at a lower reduction potential (Ecat/2 = -1.94 V vs E1/2 = -1.65 V).

Discussion

In the present work, we investigated a Fe II Fe II complex as the homodinuclear counterpart of the efficient heterodinuclear NiFe HER electrocatalyst [L N2S2 Ni II Fe II Cp(CO)] + , in which the Ni II ion has now been replaced by Fe II . The objective of this study is to evaluate the influence of the nature of the metal ion located in the N2S2 site on the HER reactivity. Structural differences between the two M II Fe II complexes, both in the solid state and in solution, have been observed. In the solid state, the FeFe complex could not be isolated in an isostructural form of since it contains two CO molecules: one bound to Fe' ({FeCp}) in a semi-bridging manner and one terminally bound to the Fe center ({FeL N2S2 }). In low-concentrated MeCN solutions, one CO ligand is removed and the remaining CO ligand is bridging, whereas in it is terminally-bound to the Fe' moiety. Stabilizing a four-coordinate Fe II ion in the Fe site is difficult. In the specific case of the L N2S2 ligand, even if mononuclear tetracoordinated [ML N2S2 ] complexes have been previously isolated with M = Co II , [START_REF] Gennari | Structural, spectroscopic and redox properties of a mononuclear Co-II thiolate complex -the reactivity toward S-alkylation: an experimental and theoretical study[END_REF] Cu II , [START_REF] Gennari | A Copper Thiolate Centre for Electron Transfer: Mononuclear versus Dinuclear Complexes[END_REF] Zn II 35 and Ni II , [START_REF] Gennari | Reversible Apical Coordination of Imidazole between the Ni(III) and Ni(II) Oxidation States of a Dithiolate Complex: A Process Related to the Ni Superoxide Dismutase[END_REF] we did not succeed in isolating the corresponding [FeL N2S2 ] complex. Instead, the mononuclear [FeL N2S2 (Cl)] - and the dinuclear [Fe2(L N2S2 )(L N2S2 H)] + complexes were obtained, in which the Fe II ions are five-coordinated. [START_REF] Wang | A Non-Heme Diiron Complex for (Electro)catalytic Reduction of Dioxygen: Tuning the Selectivity through Electron Delivery[END_REF][START_REF] Wang | Solvent-and Halide-Induced (Inter)conversion between Iron(II)-Disulfide and Iron(III)-Thiolate Complexes[END_REF] Likewise, in [(N2S2)(Fe-NO)Fe(CO)Cp] + (Scheme 1) Darensbourg et al. needed a NO co-ligand to complete the N2S2 coordination of one Fe center. [START_REF] Ding | Hemilabile Bridging Thiolates as Proton Shuttles in Bioinspired H2 Production Electrocatalysts[END_REF] For [(dppe)Ni(-pdt)Fe(CO)3], an additional CO was required to isolate the homobimetallic congener [(dppe)(CO)Fe(-pdt)Fe(CO)3], as in the present case. [START_REF] Gloaguen | Small molecule mimics of hydrogenases: hydrides and redox[END_REF] Since our aim is to evaluate the role of the metal in the catalytic process, the most relevant species to be compared are the one-electron reduced complexes, , which represent the entrance in the proposed catalytic cycles. Even if only one CO is present in both structures, it is terminally bound to the {FeL N2S2 } moiety in , while it bridges the Ni and Fe ions in . The electronic structure is also notably different since displays a mixed-valence character of type II, whereas in the spin-density is almost exclusively localized on the Ni site. Concerning the electrocatalytic HER process, it occurs in both cases through an E(ECEC) mechanism: for , cyclic voltammetry experiments carried out in the presence of Et3NH + (see Figure 5) demonstrate experimentally this pathway, while for it was mainly supported by DFT. The overpotential requirement of the electrocatalyst ( = 730 mV) is slightly higher than that of ( = 690 mV). Indeed, even if the second E-step occurs at a less negative potential for , the catalytic wave develops at a lower potential for the FeFe compound (see above). Finally, the FeFe and NiFe systems display comparable catalytic stabilities in MeCN solution (TON = 15, TON = 16). All together, these data suggest that even if the electronic and structural parameters of the key intermediates, i.e. the one-electron reduced species slightly differs, the nature of the metal ion located in the N2S2 ligand does not significantly impact the catalytic HER efficiency and an E(ECEC) pathway is proposed in both cases. The reactivity is most likely dominated by the coordinating and electronic properties of the L N2S2 ligand itself, especially owing to its non-innocent redox character. In fact, for both NiFe and FeFe electrocatalysts, the second reduction process is proposed to occur on the bipyridine unit of L N2S2 . In addition, DFT calculations carried out for the NiFe system predict that the third reduction process also occurs on the bipyridine unit. Since, a similar mechanism is suggested for the FeFe system, it clearly appears that the redox activity of the L N2S2 ligand, especially its bipyridine unit, partly drives the entire catalytic process. Interestingly, the present study describes a rare example of an [FeFe] H2ase model that mimics the iron oxidation states relevant to enzyme reactivity. The mimics the Hox state (scheme 1) and an ECEC pathway generates H2 from both and the Hox state. By contrast, most of the reported [FeFe] H2ase models involve more reduced complexes, Fe I Fe I or Fe 0 Fe I , stabilized through the coordination of phosphines and/or multiple CO ligands, as the resting state. [1][START_REF] Schilter | Hydrogenase Enzymes and Their Synthetic Models: The Role of Metal Hydrides[END_REF] The intermediate formed during the electrocatalytic cycle parallels the generation of the super-reduced state, proposed to be involved during the catalytic cycle of [FeFe] hydrogenase.40 In this case, the bipyridine moiety acts as an electron reservoir for the FeFe catalyst, like the proximal [4Fe4S] cluster in [FeFe] H2ase.41, 42

Regarding the protonation steps, especially the one that is proposed to occur after the generation of the two electron reduced species, , it can lead to a Fe-hydride species or to a thiol-containing intermediate, . Based on several DFT investigations, it has been proposed that such type of (hetero)dinuclear complexes with bridging thiolates display the propensity to protonate a thiolate during the catalytic process followed by the decoordination of the corresponding thiol, because of the hemilability properties of these metallodithiolate ligands.(darensbourg et hall, acs catal) 24, 26-27, 37 . In contrast, in the present case of , we propose the generation of a terminal hydride on the Fe' site based on the fact that (i) a metal-hydride is generated during the first protonation of the parent electrocatalyst, and (ii) the Fe' site is proposed to have a vacant site based on the DFT-predicted structure of .

Conclusion

In summary, we succeed to compare the activity between a [FeFe] H2ase mimic and its previously reported [NiFe] analogue, two dinuclear HER catalysts in which only the nature of the metal (Ni vs Fe) in the {ML N2S2 } metallodithiolate ligand is modified.

Although the structure of the catalytic intermediates, especially the one electron reduced species, notably differs, the activity of the two electrocatalysts is similar. This might be explained by the non-innocent redox activity of the L N2S2 ligand that is proposed to dominate the catalytic process. Further studies are needed at probing the full catalytic cycle, deciphering between terminal or bridging hydride ligands as catalytic intermediates and investigating if hemilability of the metallodithiolate ligand is at play in the H2 evolution process mediated by .

Figure 1 .

 1 Figure 1. (left) X-ray crystal structure of ClO4•MeCN with partial thermal ellipsoids drawn at 30% probability. For clarity, the co-crystallized solvent and the perchlorate counterion have been omitted and the aromatic groups are drawn as lines. (right) Zero-field 57 Fe Mössbauer spectrum (experimental + simulation) of solid recorded at 80 K.

Figure 2 .

 2 Figure 2. (left) IR spectra of in the solid state (black), in CH2Cl2 (0.4 mM, red), and in MeCN (1.0 mM, blue; 0.2 mM, pink) solutions; (right) UV-vis spectra of in the solid state (black), in CH2Cl2 (0.4 mM, 2 mm path length, red), and in MeCN (1.0 mM, 2 mm path length, blue; 0.2 mM, 10 mm path length, pink) solutions.

Figure 3 .

 3 Figure 3. 1 H NMR (500 MHz) spectra of dissolved in CD2Cl2 (0.4 mM, above, all the shown peaks are assigned to ) or CD3CN (1.0 mM, below, peak attribution: blue triangles = ; red circles = ).

Figure 4 .

 4 Figure 4. (left) UV-vis spectra of dissolved in MeCN (0.5 mM, 5 mm path length) before (blue line), during (grey lines) and after (red line) bubbling CO. The COsaturated solution of is successively bubbled with argon and the corresponding UV-vis monitoring is shown on the right (before Ar-bubbling: red line, during: grey lines, after: magenta line).

Figure 6 .

 6 Figure 6. Zero-field57 Fe Mössbauer spectrum (left) at 80 K of in solid state. X-band EPR spectra (right) of (black solid, in solid state at 22 K; black dash, in frozen MeCN solution at 7 K; red dash, in frozen THF solution at 100 K; red solid, in solid at 100 K.

-Fe 1 . 5 -

 15 Fe

5 -Figure 7 .

 57 Figure 7. (left) Difference electron densities sketch of the NIR TDDFT-predicted transition for * at 1017 nm (yellow = negative, red = positive), (right) localized SOMO of *.

{FeFe'} + {FeFe'} + 2CO CO {FeFe'} CO

								Ph	Ph		+
								OC
										S
								N	Fe II	S	Fe II
								N		C O	only form
		Ph	Ph		Ph	Ph	+		Ph Ph	present in CH 2 Cl 2		OC Ph	Ph
	1/2	N N	Fe II S S S Ph Ph H	Fe II S Ph Ph	N N	2 [CpFe(CO)(MeCN) 2 ] + MeCN, RT	Ph -CO	Ph + CO equilibrium in MeCN	+	CoCp 2 MeCN, RT	N	N	S Ph Ph Fe 1.5+ S	Fe 1.5+
			[Fe 2	SH ] +			(MeCN)		S
								N	Fe II	S	Fe II
								N		C O
									Ph Ph
							•••C44 and Fe2-C44 distances of 2.345(6) Å and
	1.770(6) Å, respectively), like in the enzyme. Indeed, the crystallographic structure of
	the [FeFe] hydrogenase from D. desulfuricans (presence of two non-equivalent
	molecules) displays Fe-CO distances of 2.40 Å (2.59 Å) and 1.69 Å (1.69 Å) in the Hred
	state. 30

Table 1 .

 1 Experimental and (TD)-DFT predicted spectroscopic parameters of the different characterized complexes and DFT-optimized structures ( * ).

		Experimental values	Calculated values
	IR (CO in cm -1 )
	( * )	1959, 1918 (solid)	1960, 1900
	( * )	1822 (MeCN)	1801
	( * )	1822 (MeCN)	1838
	( * )	1896 (solid) a	1904
	NIR (in nm)
	( * )	1015	1017 b
	Mossbauer
	(δc (∆EQ) in mm . s -1 )

( * ) Fe 0.23(0.69) Fe' 0.35(2.32) Fe 0.26 (0.82) Fe' 0.39 (2.37) ( * ) Fe 0.33(0.44) Fe' 0.51(2.19) Fe 0.27 (0.64) Fe' 0.58 (1.91) EPR (g) ( * ) @ 100 K, giso = 2.064 @ 22 K, gǁ = 2.262 g⊥= 1.947 2.088 (1.991, 2.124, 2.149)

Acknowledgements. The authors gratefully acknowledge research support of this work by the China Scholarship Council (LW), the French National Agency for Research in the framework of the "Investissements d'avenir" program (ANR-15-IDEX-02), the Labex ARCANE, the CBH-EUR-GS (ANR-17-EURE-0003) and the ANR-DFG (ANR-16-CE92_0012_01), and the Deutsche Forschungsgemeinschaft (DFG Me1313/14-1, NiFeMim).


Alexander Schlaich
email: alexander.schlaich@univ-grenoble-alpes.fr

Benoit

Benoit Coasne
email: benoit.coasne@univ-grenoble-alpes.fr

Dispersion truncation affects the phase behavior of bulk and confined fluids: coexistence, adsorption and criticality

published or not. The documents may come     Dispersion truncation affects the phase behavior of bulk and confined fluids: Coexistence, adsorption, and criticality

I. INTRODUCTION

Nanoporous materials are widely used in industry as adsorbents for filtration, separation and catalysis, where understanding the thermodynamics and transport of the confined fluids is crucial.

Typically, fluids confined within nanometric pores (the size of a few molecular diameters) exhibit properties that are significantly different from their bulk counterpart. In particular, surface forces and the reduced dimensions of the system affect phase transitions due to capillary condensation, freezing, etc. For obtaining microscopic insights on the confinement effects, computer simulations have proven as helpful tool and shall be discussed in the present article in the context of predicting the thermodynamics of the confined systems.

For a typical pair potential, the computational cost increases roughly as N 2 with N being the number of particles in a simulation. [START_REF] Allen | Computer Simulation of Liquids[END_REF] A tremendous reduction can be achieved by assuming that the largest contribution to the potential and forces is due to the closest neighbors only. Then, a spherical cutoff is typically applied. This means, that the pair potential φ (r) is set to zero for r ≥ r cut , where r cut is the cutoff distance. In this article, we focus on the Lennard-Jones (LJ) potential, which is probably the most widely studied intermolecular interaction potential in computer simulations and performs realistically for small spherical and nonpolar molecules. It has been employed since more than 50 years in order to study phase equilibria as well as interfacial and transport properties. [START_REF] Rahman | [END_REF][3][4][5] The LJ potential is given by

φ (r) = 4ε σ r 12 - σ r 6 , (1)
where r is the distance between two interacting particles, ε is the interaction strength and σ determines the equilibrium separation. Equation (1) typically is truncated at a value r cut ≤ L/2, where L is the (shortest) length of the simulation box. In a bulk liquid, and for sufficiently large values of r cut , the pair correlation function can be assumed to be unity, g(r) = 1 for r > r cut . Integrating Eq. (1) from r cut to ∞ leads to the so-called long-range tail corrections (LRC) for the energy U and pressure P. [START_REF] Allen | Computer Simulation of Liquids[END_REF] Throughout this paper. we use reduced units in which the density is measured in terms of the particle diameter σ as ρ = (N/V )σ 3 , with N the number of particles in the volume V . All energies are measured in terms of thermal energy k B T and, correspondingly, the reduced temperature and pressure follow as T = k B T /ε and P = Pε/σ 3 . The LRC for the potential energy and pressure then read as

Ũlrc = 8 9 π ρ2 1 r9 cut - 3 r3 cut , (2)
Plrc = 32 9 π ρ2 1 r9 cut - 3 2r 3 cut . (3)
Extensions of Eqs. ( 2) and (3) to non-isotropic systems can be made only for simplifying geometries, like a planar interface. 6,7 Thus, in practice, when anisotropic heterogeneous systems are simulated, the truncated and shifted LJ potential is employed,

φ ts (r) =      φ (r) -φ (r cut ) if r ≤ r cut 0 if r > r cut . (4)
The shift -φ (r cut ) is necessary to ensure consistent thermodynamic behavior 8 and a constant chemical potential of the particles in the simulations. 9,[START_REF] Schlaich | Computational Trends in Solvation and Transport in Liquids[END_REF] In this work, using Monte Carlo (MC) and molecular dynamics (MD) simulations, we study in Section II the consequences of using Eq. ( 4) instead of Eq. ( 1) for bulk thermodynamic properties such as the critical temperature, liquid/vapor saturation pressure and surface tension, with particular impact on the choice of r cut and confirm that the mean field relations given in Eqs. (2) and (3) yield accurate predictions for the shift in the bulk critical temperature T c , the vapor/liquid saturation pressure P sat and the surface tension γ due to the cutoff. We then turn to the impact of r cut on adsorption isotherms. Among porous solids with pores at the nanometer scale, silicabased materials are one of the most common and nitrogen adsorption experiments at low temperature are routinely used for the characterization of porous solids. [11][12][START_REF] Rouquerol | Adsorption by Powders and Porous Solids[END_REF][START_REF] Kruk | [END_REF][15][16][START_REF] Fan | Microporous and Mesoporous Materials SI[END_REF] We perform corresponding Grand Canonical Monte Carlo (GCMC) simulations in Section III. In the latter, particles are exchanged with a virtual reservoir of ideal gas particles at pressure P = f /ϕ, with the fugacity

f = 1/(β Λ 3 ) exp(β µ
) that relates the chemical potential µ to P via the fugacity coefficient ϕ, β = 1/k B T is the reduced temperature and Λ the de Broglie wavelength. Obviously, ϕ = 1 for the ideal gas, but in general ϕ = ϕ(P, T ), which we derive based on the LJ equation of state (EOS) in Section II D. Finally, in Section IV we discuss the temperature T cc below which capillary condensation is observed in dependency of pore width and cutoff. We show that for a given r cut , T cc (r cut )

is related to T c (r cut ), allowing to extrapolate to r cut → ∞ in which case the full LJ potential Eq. ( 1) is recovered. In Section V we summarize our results and draw conclusions which are also relevant to the experimental analysis of adsorption isotherms.

II. BULK FLUID PROPERTIES

A. Coexistence line & critical temperature

We perform Gibbs ensemble Monte Carlo (GEMC) simulations 18,19 using MCCCS Towhee. 20 In summary, GEMC simulations mimic two coexisting simulation boxes of volume V I +V II = V consisting of N I + N II = N particles, respectively. While V I , V II , N I and N II fluctuate, N and V are constant in the (NVT) ensemble employed. Three distinct types of trial moves are allowed for a Monte Carlo step: particle displacement, volume rearrangement, and particle exchange. 21 Our GEMC simulations are performed for N = 1000 up to 1500 LJ particles, where V was adjusted such that about 20% of the particles are in the vapor phase.

Figure 1 (A) shows the coexistence densities obtained from GEMC simulations for cutoff radii r cut = 3σ , 4σ and 5σ (blue circles, red squares and green triangles, respectively). The coexistence curves depend sensitively on r cut if the truncated and shifted potential given in Eq. ( 4) is used, as has been noted at several places before in the past and current literature (see e.g. Refs. 9,[22][23][24][25][26][27] ). The hollow black symbols included in Fig. 1 (A) show the corresponding results of the GEMC simulations if the LJ potential Eq. ( 1) together with the tail corrections Eqs. ( 2) and ( 3) is employed. The corresponding data points in Fig. 1 (A) including the LRC converge excellently for all studied values of r cut . We employ the Modified Benedict-Webb-Rubin equation of state (MBWR) by Johnson et al. 22,28 for the equation of state (EOS), P(ρ). Although a multitude of models has emerged over time for the LJ EOS, see e.g. Ref. 29 for a current overview, we decided to employ the MBWR model from Johnson et al. because of its wide use in literature. Deviations for the saturating vapor pressures P sat and vapor density ρ(P), which are the important quantities for our further analysis, are small between different models.

The critical temperature and density are given by the saddle point,

∂ P ∂ ρ T T c ,ρ c =0 (5)
and

∂ 2 P ∂ ρ 2 T T c ,ρ c =0, (6)
which we obtain numerically from the EOS, c.f. the solid black line in Fig. 1 (A), in perfect agreement with our GEMC simulation data. For the truncated-shifted potential, we include the tail correction given by Eq. ( 3) to the numerical solution of Eqs. ( 5) and (6), resulting in the solid colored lines shown in Fig. 1 (A), which agree remarkably with the simulation data for r cut = 4σ (circles), 4σ (squares) and 5σ (triangles) for the truncated and shifted potential given in Eq. ( 4) (filled symbols). Hollow symbols denote the corresponding simulation data if Eq. ( 1) with the LRC is employed. Solid lines show the MWBR equation of state from Johnson et al. 22 (black line) and the mean field corrections for finite cutoff (colored lines). Colored hollow symbols and dashed lines show the rectilinear density as given by the l.h.s. of Eq. ( 8). (B) Critical temperature for the LJ fluid determined from our simulations using the truncated and shifted potential (purple diamonds) and using the LRC (orange diamonds). Empty triangles denote simulation results from Panagiotopoulos. 23 Horizontal black line shows the value Tc = 1.313 of the MWBR EOS, dashed orange line T = 1.298 resulting from the hyperscaling relation of the surface tension discussed below. Dash dotted line represents the MWBR mean field prediction. Right axis shows the temperature for the LJ parameters of nitrogen for comparison. Inset: Fits of Eq. (7) to the simulation data using the Ising critical exponent β c = 0.32 for the cut and shifted potential (solid lines, filled symbols) and the LRC (dashed lines, hollow symbols). Color encoding is the same as in (A). and 5σ . However, for r cut = 3σ , the data slightly deviate due to the breakdown of the mean field assumption g(r) = 1 for r > r cut when such a small cutoff is used.

From the simulation data shown in Fig. 1 (A) we determine the critical temperature Tc from the scaling law

ρ l sat -ρ v sat = B |T -T c | β c , (7)
where ρ l sat , ρ v sat are the liquid and vapor densities at coexistence, β c = 0.32 is the Ising critical exponent and B is a fitting parameter. Fits of Eq. (7) to the simulation data are shown in the inset of Fig. 1 (B) for the truncated and shifted potential (filled symbols and solid lines) and using the LRC (empty symbols and dashed lines). As has been noted extensively in the literature, finite size effects become relevant especially as the temperature approaches T c . 24,27,[30][31][32][33][34] The latter studies report values in the range Tc = 1.294 -1.326 with significant scatter if more data from different groups are considered. 29 The accurate determination of the critical parameters is not at the heart of the current work, rather we want to discuss the influence of r cut on capillary criticality below. Fitting of Eq. ( 7) to the simulation data gives values in the range Tc = 1.29 -1.30 for the simulations including the LRC, depending on which values are included in the fit (orange diamonds in Fig. 1 (B)), in fair agreement with the literature. The orange dashed line shown in Fig. 1 (B) indicates Tc = 1.298 determined via the critical scaling of the surface tension to be discussed below. We solve Eqs. ( 5) and ( 6) numerically resulting in Tc = 1.313 and ρc = 0.310 in perfect agreement with the values set by Johnson et al. when fitting the MWBR equations, 22 shown as solid horizontal line in Fig. 1 (B). Simulations of this work using the truncated and shifted potential are shown in Fig. 1 (B) as purple diamonds together with simulation data from Panagiotopoulos (white triangles). 23 The dash-dotted line in Fig. 1 (B) shows the corresponding mean field correction to the MWBR model for the truncated and shifted LJ potential. Excellent agreement with the simulation data is observed over the full cutoff range.

The critical density can be estimated using the rectilinear diameter

ρ l sat + ρ v sat 2 = ρ c + A (T -T c ) , (8)
where A is another fitting parameter. The rectilinear densities obtained from the GEMC simulations in this work are shown as colored hollow symbols and the predictions of the MWBR equation as dashed lines in Fig. 1 (A). The critical density is independent of r cut in our simulations within the statistical uncertainty and in line with the MBWR mean field correction for a sufficiently large cutoff, see Supplemental Figure S1. The mean field corrected EOS, however, underestimates ρ c for r cut 3.5σ (see Supplemental Figure S1), whereas T c is well predicted also for small r cut in Fig. 1 (B). In other words, T c , corresponds to the energy scale to bring a particle from the liquid to the vapor phase (or vice versa) and is strongly affected by the missing tail energy, Eq. ( 2); ρ c on the other hand is dominated by short range correlations and thus turns out to be nearly independent of r cut . The solid black line denotes the Antoine law, Eq. ( 10), with parameters taken from Kofke. 31 Solid squares denote results from direct MD simulations of liquid/vapor coexistence for r cut /σ = 3, 4, 5 and 6, respectively (blue to purple symbols). Circles denote results from GEMC simulations, crosses include the long-range correction Eq. ( 3). (B) Change of the saturation pressure for finite cutoff r cut , normalized to P sat using the parameters obtained by Kofke 31 for simulations including the long-range correction. Solid line denotes the shift in the saturation pressure due to the missing long-range contribution in the mean field approximation, Eq. ( 11). The inset shows the same data as in (a) but now also normalized to P sat

B. Liquid/vapor saturation pressure

As seen in the previous section, the choice of the cutoff r cut has drastic consequences on the thermodynamic properties when using a truncated and shifted LJ potential. Considering pore filling and capillary condensation, the precise knowledge of the liquid/vapor saturation pressure P sat is of great importance. In this work, we evaluate the pressure from the configurations of the GEMC simulations from the virial microscopic expression via the trace of the pressure tensor 1

P = Tr T ρ + 1 3 Ṽ ∑ i< j Fi j • ri j , (9)
which we analyze separately in the vapor and the liquid volumes; equality of both pressures is also a measure for proper equilibration of the simulations. The simulation data are shown in Fig. 2 (A) and reveal a monotonous increase of Psat with increasing temperature and decreasing cutoff. We additionally perform Molecular Dynamics (MD) simulations of the liquid-vapor interface in the NV T ensemble with N = 3750 LJ particles (shown as squares in Fig. 2 (A)) in a cubic box of length L = 16. To exclude finite size effects, we checked for selected data points that simulations with N = 15000 and L = 32 yield identical results. All MD simulations were carried out using LAMMPS 35 with a time-step t = t ε/(mσ 2 ) = 0.9 • 10 -4 and periodic boundary conditions.

The data obtained from GEMC and MD simulations agree excellently, thus strengthening that our results are robust w.r.t. finite size effects. Compared to GEMC simulations, MD yields the advantage of straight-forward parallelization and fast phase-space sampling for a LJ fluid. The saturation pressure of the LJ fluid with LRC is well described by an Antoine law,

log P sat = A - B C + T ( 10
)
where A = 3.31885 B = 7.31828 and C = 0.039433 are established fit parameters obtained by Kofke, 31 c.f. black line in Fig. 2 (A). The simulation data including the LRC (crosses in Fig. 2 (A))

agree well with Eq. ( 10) with significant deviations only for r cut = 3σ at low temperatures (see inset of Fig. 2 (B)). Colored lines in Fig. 2 (A) and inset of (B) show the corresponding mean field prediction following from Eq. ( 3),

Psat = P sat 1 - 32 9 π ρsat ρ sat 2 1 r9 cut - 3 2r 3 cut , (11)
where we have used ρsat / ρ sat = 1 in good agreement with the simulation data shown in Supplemental Figure S1.

The average normalized saturation pressure P sat /P sat for the different values r cut is shown in Fig. 2 (B) together with the prediction from Eq. ( 11). The temperature-dependence shown in the inset of Fig. 2 (B) reveals that deviations appear at low temperatures and with decreasing cutoff, corresponding to the breakdown of the assumption g(r) = 1 for r > r cut . As discussed in detail in the Supplemental Material, the relative error ( Psat -P sat )/ P sat is typically about 5% and, thus, by all practical means negligible for the interpretation of adsorption isotherms to be discussed below.

Significant deviations of up to 20% appear, however, for r cut = 3σ at low temperatures, see the blue data in the inset of Fig. 2 (B). Importantly, the overall agreement does not improve if an extended Antoine law for P sat (T ) is used, see Supplemental Figure S2, strengthening the breakdown of the mean field approximation for r cut 3σ .

C. Surface tension

The liquid-vapor surface tension γ further characterizes the thermodynamic properties of the fluid and becomes relevant for the analysis of phase coexistence in confinement to be discussed below. For the Lennard-Jones potential, γ has been studied extensively since more than four decades and it has soon been realized that tail corrections similar to Eq. (3) contribute significantly. 5,6,26,[36][37][38][39][40] However, most studies report the surface tension explicitly taking into account some tail correction, for which usually a flat interface with a tanh-shaped density profile is assumed and the corresponding long range correction is added a posteriori. Note that this approach explicitly relies on the assumption that the densities ρ l and ρ v in the liquid and vapor phase do not change significantly with r cut , in strong contrast with simulation data and the mean field corrected EOS shown in Fig. 1 (A). In a general approach, the tail energy equivalent to Eq. ( 2)

due to the missing interactions behind a distance r cut must be taken into account on the fly in simulations, and the explicit terms are available only in simple geometry. 7 We thus consistently use the truncated and shifted potential Eq. ( 4) in this work, allowing for interfaces of arbitrary geometry.

From our MD simulations we extract the surface tension from the normal and tangential components of the virial pressure Eq. ( 9), 1

γ = 2 (P N -P T ) , (12)
where the factor 2 accounts for the fact that due to the use of periodic boundary conditions we have two surfaces in our system and γ = γσ 2 /ε. Simulation results γ( T , r cut ) are summarized in Fig. 3. Near the critical points one expects the surface tension to vanish as 42

γ(T ) = γ 0 (1 -T /T c ) 2ν . (13)
A fit to all simulation data and using T c (r cut ) as determined in the previous section shown in the inset of Fig. 3 (A) indeed yields the critical exponent of the correlation length ν = 0.63, in perfect agreement with the three-dimensional lattice gas model. [START_REF] Le Guillou | [END_REF][44][45] The convergence of the simulation data onto Eq. ( 13), shown as solid lines in Fig. 3 (A), again justifies the use of the mean field approximation in our transformation of the EOS to the truncated-shifted potential with cutoff. We explicitly checked that the obtained surface tension does not depend on the system size by doubling the lateral area to 2 L = 36 for r cut = 4σ , shown as red crosses in Fig. 3 (A).

We also include in Fig. 3 the data by Mecke et al. 38 for r cut = 5σ and L ≈ 12 as green crosses.

Their data point at T = 1.15 is slightly off, suggesting that, when approaching the critical point Tc (r cut = 5σ ) = 1.28, finite size effects might get important, which is consistent with the deviation of our data for r cut = 3σ at high temperature (blue data in Fig. 3). Contrary, the strong deviation from the scaling behavior Eq. ( 13) at low temperature for this small cutoff is attributed to a failure of the mean field model as discussed above. 38 Black symbols denote reference data from the National Institute of Standards and Technology using Transition Matrix Monte Carlo including long range LJ interactions via the lattice sum method and with rigorous finite size scaling. [START_REF]NIST Standard Reference Simulation Website, NIST Standard Reference Database[END_REF] For a classical van der Waals fluid the interfacial density profile on a mean field level has the form ρ(z) = ρ c -∆ρ tanh(z/(2ξ )), where ∆ρ = ρ lρ v and ξ is the interfacial thickness. Using a tanh profile ρ(z) and the Kirkwood-Buff formula for the surface tension, 46

Values of the interfacial width ξ are obtained from the simulation data and shown in Supplemental Fig. S3 (A). For r cut ≥ 4σ the width ξ moderately depends on r cut for T < 1, consistent with the tail corrected surface tension γ + γlrc for our simulation data in Fig. 3 (B). Contrary, r cut = 3σ , ξ differs by about 10% from the larger cutoffs and the corresponding surface tension in Fig. 3 (B) does not converge to the data points for larger cutoff. We fix the critical exponent ν = 0.63 for the tail corrected data with r cut ≥ 5σ to obtain the critical temperature T c = 1.298 in the inset of Fig. 3 (B), in excellent agreement with the GEMC simulation results using the critical exponent of the order parameter in Fig. 1 (B). Whereas the precise value of T c for a LJ fluid is still under debate, the Ising critical exponents have been determined with high precision and we suggest that, using a priori corrections 7 with sufficient cutoffs, the estimated surface tension is reliable and accurate and can thus actually be used to determine T c with high precision. Note that we employ the MWBR EOS to determine T c (r cut ), in good agreement with the simulation data in Fig. 1 (B) and satisfactory convergence of the simulation data in Fig. 3 (B). Thus the parametrization of the MWBR EOS using T c = 1.313 yields reliable estimates T c (r cut ), although we here suggest a critical temperature of T c = 1.298. Using T c one obtains a slightly higher exponent ν = 0.65, see Supplemental Fig. S3 (B), i.e., considering the uncertainty in T c the analysis of the surface tension scaling does not give reliable estimates for the critical exponent.

D. Fugacity coefficient of a LJ fluid

In grand-canonical Monte Carlo (GCMC) simulations, the chemical potential µ of a fictive ideal gas reservoir is prescribed and particle exchange and translocation moves are performed according to a Metropolis criterion. The reduced chemical potential μ = µ/k B T is linked to the fugacity via

μ = ln f + ln Λ3 , (15)
with the de Broglie length Λ = Λσ = h 2 B T ) and the fugacity is rescaled according to f = f ε/σ 3 . In general, the pressure P for an interacting system differs from f and one thus defines the fugacity coefficient ϕ = f /P, where for ϕ = 1 the ideal gas law is recovered. Note that the same correction must be applied in a volumetric adsorption measurement, where the adsorbed amount is often calculated from an ideal gas law at measured pressure P. 47 Typically, the density and pressure are assumed to be sufficiently low such that ϕ ≈ 1, yet the correction factor for nitrogen adsorption at 77 K and close to P sat is about 5%, in line with the MWBR EOS and our simulation results discussed below. 48 The use of the fugacity in common adsorption isotherms does not have a major effect as increasing temperature at fixed pressure makes the ideal gas approxi- 1) together with the LRC and (B) the truncated-shifted potential Eq. ( 4) obtained from GCMC simulations. Solid lines show the fugacity coefficient, determined using Eq. ( 16) and the MWBR equations.

Dashed lines denote the corresponding result from the virial expansion up to second order. Temperature in brackets show physical units for the nitrogen parameters, c.f. Table I. mation better. However, it is common practice to normalize the isotherms to the corresponding saturation pressure. For increasing temperature, P sat increases (cf. Fig. 2 (A)) as well as the vapor saturation density ρ v (cf. Fig. 1 (A)), i.e., ϕ will depart from the ideal gas behavior at P sat for increasing temperature.

The fugacity coefficient follows from the (inverse) equation of state ρ(P), which we determine by numerical inversion of the MWBR equations, as

lnϕ = P 0 1 ρk B T - 1 P dP = B 2 P + O(P 2 ), (16)
where in the last step we performed a virial expansion up to second order and B 2 is the second virial coefficient, see the Supplemental Material for details. For a LJ liquid, B 2 can be evaluated analytically in terms of modified Bessel functions, 49 and we show in Supplemental Figure S4 that the second order virial EOS indeed describes ρ(P) accurately for the LJ gas at moderate temperatures T 1. For practical purposes, the virial expansion yields the advantage that B 2 can be measured straightforwardly both in experiments as well as in molecular simulations (see Supplemental Figure S4 for a comparison). Yet, at temperatures close to T c the difference amounts to more than 10%, which is clearly not negligible.

Figure 4 shows the fugacity coefficients for various temperatures and a selected cutoff r cut = 4σ obtained from GCMC simulations using LAMMPS (A) with the LRC and (B) using the cut and shifted potential, respectively. A cubic simulation box of length L = 26.6 was chosen, corresponding to up to 8200 particles. We checked that our results are not subject to finite size effects by running selected simulations at box length 2 L. Solid lines in Fig. 4 are obtained from the MBWR equation of state via Eq. ( 16), while dashed lines result from the virial expansion up to second order with the analytical virial coefficient. In all cases, the simulation data are perfectly described by the MWBR parameters, whereas with increasing temperature higher order terms in the virial expansion have to be included. One expects the virial expansion to improve with increasing temperature, but we plot in Fig. 4 the fugacity coefficient vs. the vapor/liquid saturation pressure, which according to Eq. ( 10) increases exponentially with temperature, and thus leads to stronger deviations from the ideal gas behavior with increasing temperature.

III. ADSORPTION ISOTHERMS

Having assessed the bulk properties of a LJ liquid, we now study the adsorption of a LJ fluid in an amorphous slit pore. In detail, we consider N 2 at 77 K, like is commonly used for the characterization of porous solids. [11][12][START_REF] Rouquerol | Adsorption by Powders and Porous Solids[END_REF][START_REF] Kruk | [END_REF][15][16][START_REF] Fan | Microporous and Mesoporous Materials SI[END_REF] A computational model of amorphous silica was obtained by quenching from cristobalite as described in Ref. 50 , see Fig. 5

(A). A finite slit pore of width

H is obtained by expanding the cubic unit cell of length a = 38.5 Å in the z-direction such that

L z = a + (H + 2σ N 2 /O
), where we account for the finite size of nitrogen such that for H 0 a N 2 molecule fits into the slab, see Fig. 5 (B). The unit cell was replicated three times in the x-direction and a reservoir of length 100 Å was added on both sides, see Fig. 5 (C). Periodic boundary conditions in all directions are used and for r cut = 5σ the system also was replicated in y-direction in order to exclude interactions with periodic images. Molecular nitrogen is modeled as a single LJ sphere with parameters σ = 3.75 Å and ε/k B = 95.2 K. 51 The Lennard-Jones parameters for the oxygen atoms of the silica substrate are σ OO = 2.7 Å and ε OO /k B T = 230 K, and following previous simulations of nitrogen adsorption in silica pores, we neglect interactions with silicon and hydrogen atoms of the substrate because of the low polarizability. 52 The cross-parameters for the interaction between the nitrogen molecules and the oxygen atoms of the pores were determined using the Lorentz-Berthelot mixing rules as σ N 2 /O = 3.225 Å and ε/k B = 147.97 K (Table I). We have previously shown that this LJ model can be used successfully for capturing the global experimental behavior upon adsorption and desorption in silica nanopores. 53 All simulations were performed for at least 10 6 MC steps per nitrogen molecule of which 10% were discarded as equilibration. To minimize boundary effects of the finite-sized pores, the adsorbed amount Γ has been determined only in the central pore region of length a, c.f. Fig. 5. In Figure 6 we report the adsorption/desorption isotherms for three cutoff values r cut = 3, 4 and 5σ and three pore widths H = 1.2, 2.2 and 3.2 nm. In all figures, the adsorbed amount Γ has been normalized to the number of molecules in the central filled pore at P sat (T ), Γ 0 , and plotted vs. the reduced pressure P/P sat (T ). Considering only the data for T = 0.81 (blue data in Fig. 6), which corresponds to T = 77 K in our simulations of nitrogen adsorption, the obtained isotherms conform to the typical behavior observed in experiments of adsorption/condensation in nanopores: the adsorbed amount increases continuously in the multilayer adsorption regime until an irreversible jump occurs because of capillary condensation and evaporation of the fluid within the pore. Once Table I. Lennard-Jones interaction parameters. Nitrogen molecules b and the oxygen atoms of the silica substrate b are modeled as LJ spheres, while cross terms stem from the Lorentz-Berthelot mixing rule. the pore is filled, the adsorbed amount changes only slightly upon increasing the pressure due to the small compressibility of liquid nitrogen. With increasing pore width (top to bottom) and increasing cutoff (left to right), the simulation data reveal a pronounced hysteresis loop in Fig. 6.

σ N 2 /N 2 3.75 Å ε N 2 /N 2 k B • 95.2 K σ N 2 /O 3.225 Å ε N 2 /O k B • 147.
Conceptionally, the presence of condensation and hysteresis can be sketched directly by combining the classical Kelvin and Laplace equations, [START_REF] Charlaix | Handbook of Nanophysics: Principles and Methods[END_REF] which shifts the vapor-liquid coexistence pressure according to

P v (µ) -P l (µ) = k B T ρ ln P sat P v = 2γ lv H . (17)
Here, P v and P l are the pressure in the vapor and the liquid phase, separated in confinement by a curved interface, and γ lv is the liquid-vapor surface tension of the adsorbate. Based on Eq. ( 17),

Derjaguin and Churaev have proposed a thermodynamic model describing both the formation of an adsorbed film at the pore surface and capillary condensation-evaporation. [START_REF] Derjaguin | [END_REF] Considering the simulation data for 77 K (blue symbols in Fig. 6), the condensation/evaporation pressures decrease inversely with the pore width H, as expected from Eq. ( 17) and the Derjaguin-Churaev model. 56 Noteworthy, the width of the hysteresis loop increases with H, but also with the cutoff, which we attribute to the increased stability of the adsorbed film for increasing r cut .

Whether the adsorption isotherms are determined experimentally or theoretically from molecular simulation, the data points must be fitted with analytical equations for the calculation of thermodynamic properties by numerical integration or differentiation. Instead of relying on expressions based on the underlying physical mechanisms like the Langmuir equation, the BET model or the Derjaguin-Churaev expressions, we employ here a purely phenomenological fit which captures the multilayer adsorption and capillary condensation according to

P(Γ) = w 1 1 + exp -Γ-Γ 1 ∆Γ 1 + w 2 1 + exp -Γ-Γ 2 ∆Γ 2 . ( 18
)
Note that Eq. ( 18) phenomenologically interpolates between a BET isotherm and a saturation at Γ 1,2 Γ like in a Langmuir model. Also note that Eq. ( 18) yields a Langmuir-like linear adsorption regime for Γ → 0, but with an offset P(Γ → 0) > 0 that describes the simulation data well for the strongly wetting system considered here. The latter is due to the fact that the silicanitrogen interaction is strong and thus a linear adsorption would be observed only at extremely small pressures.

The fits of Eq. ( 18) to the simulation data are shown as solid lines in Fig. 6 and reveal a drastic influence of r cut on the condensation/evaporation pressures as well as on the width of the hysteresis loop. Generally, condensation appears at lower pressures for large cutoff: Increasing r cut lowers the potential energy according to Eq. ( 2), and thus increases the relative energy difference between the adsorbed film and the filled pore. Increasing the temperature in all cases results in a transition from a discontinuous pore filling at temperature T < T cc to a continuous filling at larger temperatures T > T cc , where T cc is the capillary critical temperature and shall be discussed in the next section. Note that for regular pores such as slits or cylinders, the equilibrium condensation pressure is the one observed on the evaporation branch for wetting pores, whereas condensation occurs in the metastable coexistence region between an adsorbed film and the filled pore due to thermal fluctuations. The condensation pressures increase with the loss in energy due to the cutoff and for our finite number of MC steps condensation on the adsorption branch is not observed for r cut = 3σ at 77 K for the large pore H = 3.2 nm in Fig. 6 (G), as expected from kinetic modeling of the phase coexistence. 57 IV. CAPILLARY (PSEUDO) CRITICALITY Capillary condensation (CC) occurs due to cooperative interactions between adsorbed layers, and thus depends on temperature and pore width. Whether capillary condensation in nanopores constitutes a first order transition or not is still debated. 58,59 Pore filling occurs at a temperature T cc , that is below the bulk critical temperature T c . Above a critical temperate T cc (H) < T < T c , capillary condensation is suppressed, i.e. pore filling becomes continuous and reversible, which has been confirmed experimentally 60,61 , predicted by density functional theory 62 and is consistent with our simulation data in Fig. 6. Importantly, T cc is in general different from the temperature T ch at which the hysteresis disappears, i.e., capillary condensation still occurs but is reversible as fluctuations are large enough to prevent the system from being trapped in metastable states.

To determine if CC occurs, we evaluate the dimensionless isothermal compressibility

κT = - Γ Γ 0 ∂ (Γ 0 /Γ) ∂ (P/P 0 ) T , (19)
from a fit of Eq. ( 18) to the simulation data (lines in Fig. 4). In a true first-order transition, κ -1 T vanishes, in contrast to Eq. ( 18). However, also in experiments, the divergence of the compressibility is usually not observed, which is often attributed to pore-size heterogeneity. 61 We evaluate the full width at half maximum (FWHM) of the compressibility peak to determine whether CC takes place, i.e., whether the isotherm is discontinuous, or not. The FWHM changes by about an order of magnitude between continuous and discontinuous filling (Supplemental Figure S6) and thus provides a more rigorous definition as the numerical value of κ, see Supplemental Material for details. The resulting CC phase diagram is shown in Fig. 7 (A) for all pore sizes, temperatures and cutoff values considered. We also indicate by the orange symbols in Fig. 7 (A) the presence of a hysteresis loop if the relative difference between the condensation and evaporation pressures (i.e., the positions of the maxima in κ T ) is larger than 2%. As expected, capillary condensation occurs at temperatures lower than the corresponding bulk critical temperature T c (r cut ) determined from the MWBR EOS, shown as dashed lines in Fig. 7 (A). For the smallest pore, filling is discontinuous irrespective of T , but the hysteresis vanishes with increasing temperature.

The shift ∆T c = T c -T cc in the critical temperature is from finite size scaling arguments expected to scale as ∆T c /T c ∼ (ξ /H) 1/ν , where ξ ∼ σ is the correlation length of the density fluctuations and ν is a critical exponent. 63 Using density functional theory (DFT), Evans et al. 62

derived a relation T c -T cc T c ≈ σ H (20)
which agrees with the scaling arguments by Nakanishi and Fisher for ν → 1. 63 The DFT result given in Eq. ( 20) is shown in Fig. 7 (A) as solid line with the corresponding values T c (r cut ) determined from the MWBR EOS, revealing excellent agreement except for the smallest pore. This disagreement for the small pore is not unexpected as in the DFT calculations the fluid molecules were treated as hard spheres and in such strong confinement correlation effects will become important.

The relevance of the choice of the cutoff becomes clear when the CC data are plotted vs. σ /H as proposed by Eq. ( 20) in Fig. We also include in Fig. 7 (B) the temperature T ch at which the capillary hysteresis disappears as empty symbols. The agreement of both T ch and T cc with Eq. ( 20) is remarkable and reveals a shift due to the cutoff ranging from about 3 K for r cut = 5σ , 5 K for r cut = 4σ to a value as large as 15 K for r cut = 3σ . Remarkably, all simulation data converge on the ∆T c ∼ H -1 law, Eq. ( 20), as we highlight in the inset of Fig. 7 (B).

V. CONCLUSIONS

The computational benefit of using a finite cutoff together with the truncated and shifted LJ potential leads to significant changes in the phase diagram, and importantly the critical temperature, even in bulk. Using mean field corrections, we have shown that the capillary critical temperature is well described by the T cc ∼ T c (1σ /H) law found in DFT calculations and experiments if the consistent value T c (r cut ) is used. Care, however, has to been taken in the quantitative comparison of adsorption isotherms: On the one hand, the chemical potential, and thus the fugacity f , determines the equilibrium condensation and are straightforwardly accessible both in simulations and volumetric adsorption experiments. On the other hand, it is common practice to present the adsorbed amount vs. the pressure P normalized to the vapor/liquid saturation pressure P sat , revealing two possible pitfalls, both in experiments and simulation. First, the saturation pressure depends on temperature, which for a simple fluid like nitrogen is well accounted for by Antoine's equation, and in simulations the pressure correction due to r cut has to be taken into account. Second, and more importantly, f and P are related via the fugacity coefficient ϕ. For nitrogen at 77 K typically ϕ ≈ 1 is assumed, leading to systematic deviations close to P sat of about 5%. For higher temperatures as used for the analysis of pore criticality, ϕ becomes non-negligible and is about 0.7 close to P sat . It is therefore of utmost importance to precisely specify if fugacities or pressures are reported.

Summarizing, we have given a comprehensive guide on how to properly take into account truncation effects in molecular simulations whenever the isotropic mean field corrections cannot be applied on the fly, including critical temperature, saturation pressure, surface tension and fugacity coefficient. This allowed us to quantify cutoff effects on modeling adsorption isotherms of fluids in porous media. Performing simulations at different temperatures, we were able to rigorously determine the capillary critical temperature T cc based on the divergence of the compressibility for a first order transition, which reveals a universal scaling of T cc with the inverse pore width independent of the detailed interaction parameters.

Figure 1 .

 1 Figure 1. Liquid/vapor coexistence curve and critical temperature. (A) Coexistence data for r cut = 3σ

Figure 2 .

 2 Figure 2. Liquid/vapor saturation pressure. (A) Saturation pressure Psat as a function of the temperature.

Figure 3 .

 3 Figure 3. Surface tension of the LJ liquid/vapor interface. (A) MD Simulation results (filled symbols) using the truncated-shifted potential for various cutoff values. Simulations are performed in a box of size L × L × 2L with L/σ = 16. Red crosses show additional simulations for L/σ = 32 and r cut /σ = 4. Green crosses denote simulations by Mecke et al. for r cut /σ = 4 in a box L × L × 3L with L/σ ≈ 12. 38 Inset shows the hyperscaling relation Eq. (13) with the critical exponent ν = 0.63, solid line shows the slope 2ν = 1.26 as guide to the eye. Rescaling Eq. (13) with the cutoff-dependent critical temperature the colored solid lines are obtained. Solid Black line extrapolates T c = 1.313 to the full LJ potential from the EOS. (B) Simulated surface tension including the long range correction Eq. (14). The critical temperature T c = 1.298 was fitted using ν = 0.63 (inset), the corresponding curve according to Eq. (20) is shown as dashed black lines in (A) and (B). Dash-dotted lines denote fitting results of Mecke et al. with ν = 0.63 and T c = 1.325 for comparison. 38 Black symbols denote reference data from the National Institute of Standards and Technology

6 γlrc

 6 Blokhuis et al. have derived the long range contribution of the LJ potential assuming that ξ is independent of r cut , r -3 (3s 3s) coth sr 2ξ .

Figure 4 .

 4 Figure 4. Fugacity coefficient of a LJ liquid. Simulation data for r cut = 4σ using (A) the LJ potential Eq. (1) together with the LRC and (B) the truncated-shifted potential Eq. (4) obtained from GCMC sim-

Figure 5 .

 5 Figure 5. GCMC simulations. (A) Molecular model of amorphous silica as obtained by quenching from cristobalite. 50 In our simulations, only oxygen atoms (orange spheres) interact with molecular nitrogen, shown as blue spheres in (C). (B) A slit pore of width H is obtained by extending the periodic box in the z-direction according to L z = a + (H + 2σ N 2 /O ). (C) Typical molecular configuration obtained by means of Grand Canonical Monte Carlo simulations at the onset of capillary condensation for N 2 (blue spheres) adsorbed at T = 77 K and P = 0.7P sat in a finite silica slit pore of width H = 1.2 nm with an external reservoir of length 10 nm on each side.

r

  cut = 3ff A T =Tc(rcut) 0:70 (77 K) 0:75 (82 K) 0:79 (87 K) r cut = 4ff B T =Tc(rcut) 0:65 (77 K) 0:70 (82 K) 0:74 (88 K) r cut = 5ff C T =Tc(rcut) 0:63 (77 K) 0:68 (83 K) 0:76 (93 K)

Figure 6 .

 6 Figure 6. Truncation effect on adsorption isotherms. Adsorption isotherms of nitrogen in silica nanopores are shown as empty symbols, desorption isotherms as filled symbols for H = 1.2 nm (A)-(C), 2.2, nm (D)-(F) and 3.2 nm (G)-(I). Columns show data for different values r cut . Solid lines denote fits of Eq. (18) to the desorption data, dashed lines to adsorption. Colors denote the temperature normalized to the bulk critical point T c (r cut ), corresponding values for the nitrogen parameters are given in brackets for clarity.

Figure 7 .

 7 Figure 7. Influence of temperature and cutoff on pore filling. (A) Capillary condensation phase diagram indicating whether pore filling is continuous (open circles) or not (filled circles). Continuity of pore filling has been determined via the FHWM of the compressibility at the condensation pressure for r cut = 3, 4 and 5σ in (i)-(iii), respectively. Symbols indicate if a hysteresis loop is observed (x) or not (-). Dashed lines show the corresponding bulk critical temperature T c (r cut ). Solid lines denote the shift of T cc according to Eq. (20). (B) T cc as function of the inverse pore width. Solid line denotes the capillary shift of the critical temperature according to Eq. (20) using the LJ potential with long-range correction, broken lines denote the corresponding values T cc (r cut ) for the truncated and shifted potential. Filled symbols show T cc , empty symbols the disappearance temperature of the hysteresis loop, T ch . Right scale shows the temperature in physical units for the nitrogen parameters. The inset shows the same data but for the relative change ∆T c /T c .

  7 (B). Let us first discuss the drastic shift ∆T c due to the confinement if the LJ potential with LRC is employed, black line in Fig.7 (B). Using the nitrogen simulation parameters (right temperature scale), the critical temperature reduces from its bulk value T c = 125 K to T cc = 102 K for the pore of width H = 1.2 nm. The simulation data for finite cutoffs using the truncated and shifted potential defined in Eq. (4) are shown as solid symbols, together with corrected bulk temperature T c (r cut ) obtained from the MWBR EOS (broken lines in Fig.7 (B)).

F 6 -

 6 ts (r) = φ (r cut )δ (rr cut ), r ≤ r cut 0, r > r cut .

ACKNOWLEDGMENTS

We acknowledge funding from the ANR project TAMTAM (ANR-15-CE08-0008-01).