Rui Shang

Stephan N Steinmann

Bo-Qing Xu

Philippe Sautet

Mononuclear Fe in N-doped carbon: computational elucidation of active sites for electrochemical oxygen reduction and oxygen evolution reactions

Non-precious metal catalysts are well investigated in electrocatalysis. Fe and N codoped carbon (Fe-N-C) catalysts have drawn a great attention due to their low cost and good performance in oxygen reduction reaction (ORR) and recently in oxygen evolution reaction (OER). Based on the recent advances of a variety of physical characterization techniques, more information about the chemical environment of the catalytic active sites has been acquired. However, due to the complexity of the catalytic material and process, the real active structures are still controversial. In this work, several active sites are proposed for Fe-N-C catalysts (L-FeNx, x = 2 or 4, L = nothing, O or OH) and their performance in electrochemical ORR and the OER are investigated.

The computations are based on density functional theory (DFT), including Van der Waals interaction and solvation effect with an implicit electrolyte model. Calculations

indicate that the catalytic activity of the Fe centers depends strongly on the N coordination number and on the presence of extra ligands like OH group. In particular, HO-FeN2, but not FeN4, appears as the most active site. Scaling relations are obtained by connecting the free energy of potential-determining steps with the adsorption free energy of intermediates. Furthermore, three promising active sites suggested from scaling relations are studied by the more elaborate surface charging approach, which includes the influence of the applied potential and the electrolyte. The results show that the specific treatment of the influence of the applied potential has a minor influence at low potential, which is the case for ORR, but a major influence at higher potential, as

Introduction

The electrochemical oxygen reduction reaction (ORR) is a key process that controls the performance of proton exchange membrane fuel cells (PEMFC). PEMFCs are promising devices to efficiently transform hydrogen and oxygen into electricity. 1 However, ORR usually requires precious Pt metal as the catalyst and suffers from significant loss of energy, 1 which decreases the economical attractiveness of fuel cells and limits their usage.

In order to address these problems, non-precious metal catalysts have been investigated as potential alternatives. Co-phthalocyanine was first shown to have catalytic activity for ORR. 2 Inspired by that, phthalocyanine and porphyrin with various macro-cycle structures and cheap transition metals, such as Fe and Co, were studied. [3][4][5][6] Then, these transition metals were doped or co-doped with N in graphene or carbon nanotubes (CNTs) for better conductivity and stability. [6][7][8] Among these materials, the Fe and N codoped carbon (Fe-N-C) catalysts have attracted great attention due to their low costs and relatively high performance, [9][10][11][12][13][14][15] but the structure of the most active site is still controversial.

On one hand, the mononuclear Fe in coordination with four N atoms (FeN4) has been considered as the active site in most publications. [11][12][13] Several characterizations such as extended X-ray adsorption fine structure (EXAFS), 11 X-ray absorption near-edge structure (XANES), annular dark-field scanning transmission electron microscopy (ADF-STEM) 12 and aberration-corrected scanning transmission electron microscopy (ASTEM) 13 were used to shed light on the local atomic environment of these FeN4 catalysts. On the other hand, FeN2 sites 14,15 have also been shown to have high activities with experimental evidence from high-angle annular dark-field scanning transmission electron microscopy (HADDF-STEM), EXAFS and Mossbauer spectroscopy analysis. 14 These experiment observations provide important information for identifying the real active site, with a debate between FeN4 and FeN2 sites. However, none of the aforementioned characterization methods provides an atomically precise environment around the metal. Furthermore, several sites can be present on the carbon material, rendering the determination of the "real" active site responsible for the measured activity difficult. We herein show that computational chemistry provides complementary insight into the nature of active site structures.

In previous theoretical studies, several models of graphene doped with iron and nitrogen have been proposed. For instance, the ORR performance of FeN4 and FeN2-doped graphene layers was compared using Density functional theory (DFT) calculations with the Perdew-Burke-Ernzerhof (PBE) functional, indicating that the electron transfer between the 3d orbital of Fe and 2π* orbital of O is beneficial to the reaction. 16 The calculations also suggested that the FeN2+2 system, where the 2+2 index represents four N atoms that were doped in pairs on the edge of two graphene fragments, could be superior to FeN4 in extended graphene for ORR. FeN4 sites with pyridinic-type N or pyrrolic-type N were also compared and the former showed lower overpotential. 17 Furthermore, an axial ligand was shown to improve the activity, with the pyridine adsorbate strengthening the O2 adsorption 18 and the OH ligand [19][20][21] increasing the O-O bond scission or moderating the adsorption of adsorbates. In addition, existence of a metal-O bond was proposed to be essential in improving the activity. 22 In the wellaccepted four-electron mechanism, reduction of OH to H2O is considered as the potential-determining step since it is shown to be the least exothermic step. When the adsorption of OH is weakened, the formation of OOH will be the potential-determining step instead. 23 A volcano correlation has been found for the limiting potential versus the adsorption energy of OH for N-doped graphene by relating the free energy of each step and the adsorption energy of OH. 24 However, only three Fe-containing structures were computed and they were rather far from the optimum potential on the volcano.

The above mentioned computational studies rely on the computational hydrogen electrode (CHE), 25 which is the most popular approach to estimate the impact of the electrochemical potential. In that approach, the electron and proton are transferred together and only the chemical potential of exchanged electrons is assumed to be affected by the electrochemical potential, i.e., the electrode polarization is neglected and only neutral systems are considered. CHE is an easy-to-use zeroth order approximation that provides valuable information for determining the potentialdetermining step and generally compares well with experimental results. 26 Since only neutral systems are considered, most CHE studies neglect the effect of solvation as well, 16,18 even though CHE is also applicable in the presence of a solvent, e.g., Liang et al 24 included a conductor-like screening model (COSMO) model to simulate water.

As we have shown before, neglecting the explicit impact of the electrochemical potential and solvation is most dramatic when the surface dipole moment changes during a reaction, which is also reflected by a change in the workfunction of the neutral reactants, intermediates and products. 27,28 The simulation of the environment of electrocatalysts, including description of the solvent, the electrolyte and the applied potential, is also indispensable for a better understanding of electrocatalytic processes. 29 Beside the CHE method, a more detailed approach, named surface charging (SC), 30 can be applied to model electrocatalysis. In this approach, the influence of the solvent and of the applied potential are taken into consideration, which are especially important for studying adsorbates with large dipole moments. 31 A similar approach, including a continuum model of the electrolyte, 32,33 has been used on Pt(111) to understand the influence of water molecules 34,35 on ORR and also the influence of Pt particle size on the activity. 36 As for the non-noble metal catalysts, SC has been applied at a selected potential for comparing the four-electron and two-electron mechanisms on FeN4 and with a CN ligand coordinated to Fe in axial direction. 37 Understanding the influence of the potential on electrochemical reactions is important, the solvent and electrolyte can also provide refined energy profiles and more realistic insights to develop sites with high activity.

In this work, the Fe-N-C catalysts were investigated using cluster and periodic models in the framework of CHE approach, with the inclusion of Van der Waals interactions and an implicit solvent model. Several active site structures were set up by changing the number of N atoms coordinated to the Fe center and by adding an oxyl ligand (-O or -OH) to the Fe center. In addition to the investigation on ORR, the reverse reaction, i.e., the oxygen evolution reaction (OER), was studied as well since the overpotential shifts the relevant potentials in the opposite direction on the SHE scale. Although Fedoped carbon materials 38,39 have been applied in OER and experimentally found to show relatively high overpotential, 39,40 it is worth exploring the active site structure and understanding the influence of the applied potential on this reaction from a fundamental perspective. Scaling relations were obtained between the reaction energy of the elementary steps and the binding energy of the OH fragment, for various proposed sites.

Based on the obtained volcano plots for the limiting potential, FeN4, HO-FeN4 and HO-FeN2 sites were selected and recalculated by including the influence of the potential with the SC approach as a comparison with CHE. Based on these analyses, the promising active sites are identified as HO-FeN4 and HO-FeN2, with overpotentials at around 0.40 and 0.60 V for ORR and OER, respectively. The scaling relation offers key information about the catalytically active sites and generates guidelines for the design of more promising catalyst. The present comparisons between the CHE and SC approaches shed light on understanding the impact of the potential and microenvironment of Fe and N co-doped carbon in electrocatalysis.

Methods and Models

The well-accepted mechanism for ORR and OER is the four-electron mechanism, including chemisorbed OOH, O and OH as intermediates (Fig. 1). The detailed mechanism is provided in the Electronic Supplementary Information (ESI).

All calculations were performed in the framework of spin-polarized DFT by using the Vienna Ab initio Simulation Package (VASP). 41,42 The generalized gradient approximation (GGA) in the formulation of the PBE 43 functional was applied together with electron-ion interactions obtained through the Projector-Augmented Wave (PAW) 44 method. The cutoff energy was set as 400 eV. The convergence criteria of the energy and the force were set as 10 -6 eV and 0.01 eV/Å, respectively. The Van der Waals interactions were described using the dDsC method. 45,46 The influence of the solvent was modeled using an implicit dielectric model with the dielectric constant of water being 78.4. The electrolyte was represented using the Poisson-Boltzmann (PB) continuum approximation. The cavity surface tension was set as 0.00 and the electrolyte concentration within the linearized PB equation was set as 1.0 M. 47 Several structures on the graphite surface (L-FeNx, x = 2 or 4, L = O or OH) were considered, using both cluster and periodic models (Fig. S1). For the H terminated cluster models, a box of 15 × 15 × 20 Å was used and the electronic structure described at the gamma point. In the periodic models, the Fe and N were co-doped at a perfect graphite surface layer. For creating a FeNx site, six carbon atoms were replaced by one Fe, x N and 4-x C atoms. A 4 × 4 super cell with 20 Å vacuum layer in the z direction was used. The Brillouin zone was integrated by a 3 × 3 × 1 Monkhorst-Pack grid.

In the SC approach, a five-layer graphite slab was used with a 3 × 3 super cell, with the Fe and N co-doped in the top and bottom layer in a symmetric manner. All slabs are hence symmetric in order to avoid a net dipole and ambiguity in calculating the work functions. During the optimization, the middle layer was fixed and the other layers were allowed to relax. The solvent spacing between the slabs was set as 45 Å in order to avoid overlap between the electrolyte densities on each side of it. 47 The Brillouin zone was integrated by a 5 × 5 × 1 Monkhorst-Pack grid.

For computations in periodic models, we first computed the energy of each spin state for each intermediate (Fig. S2), as to select the most stable spin state for the free energy diagram. Without any adsorbate, FeN4 is most stable in the triplet state, followed by the singlet (about 0.50 eV higher) and quintet state, which is in the same sequence compared to high-level computations of iron porphyrin. 48 The most stable spin states for the intermediates were chosen to plot the free energy diagram for the ORR and OER.

(Fig. S3).

Treating Van der Waals interaction and solvation effect is essential for modeling Fe and N co-doped electrocatalysts. Van der Waals interaction and solvation effect both stabilize *OOH, *O and *OH intermediates, with a total effect of 0.30-0.50 eV, and a contribution of 40% to 60%, respectively. However, a universal scaling relations were obtained with the presence or absence of these two effects, where a similar small variation of the scaling relations was observed earlier in the case of the energy difference between OOH and OH adsorption energies on precious metal surface 49 when changing the exchange correlation functional.

Results and discussion

Active mononuclear Fe sites on periodic and cluster models with the CHE approach

The FeN4 structure (Fig. 2a) has been frequently considered as a potential active site for ORR. 16,17 In the framework of CHE, following the four-electron mechanism of ORR, the fourth step (*OH + H + + e -→ H2O + *) has usually been found as the potentialdetermining step. 50 Our results (Fig. 2a black line) show that for this FeN4 structure, the potential-determining step depends on the level of theory. If a standard PBE exchange functional is used and the influence of the solvent is not included, the third step (*O + H + + e -→ *OH) is found as the potential-determining step, with a difference of 0.03 eV in free energy with respect to the fourth step (*OH to H2O). In contrast, when Van der Waals interaction and solvation effect are included, the adsorption of OH is stabilized and the fourth step (*OH to H2O) becomes the potential-determining one. The overall stabilization energy of OH is calculated as 0.31 eV, with contributions of Van der Waals interaction and implicit solvent of 0.13 and 0.18 eV, respectively. By looking at the reaction in the reverse direction, which corresponds to OER, the most endothermic step will be potential-determining and this is the third one (*O + H2O → *OOH + H + + e -), in accordance with the earlier literature. 51 If we get back to the ORR, since the third (*O to *OH) and the fourth (*OH to H2O) steps are weakly exothermic at zero potential, and strongly endothermic at equilibrium potential for FeN4 (Fig. 2a), further adsorption and reaction of O2 on O=FeN4 or HO-FeN4 might compete with these steps. Hence, in this case, the first reaction cycle would not complete and instead O=FeN4 or HO-FeN4 would form as a new active site (structures in Fig. 3). Such catalytic cycles starting from O=FeN4 or HO-FeN4 have also been considered and depicted in red and blue respectively in Fig. 2 (see also Table S1 and Table S2). The first step (*O2 + H + + e -→ *OOH) is even more endothermic (Fig. S4a) in the absence of Van der Waals interaction or solvation effect, but if these are properly included, the formation of O=FeN4-OOH or HO-FeN4-OOH are competitive to or much more favorable than the third (*O to *OH) or the fourth to H2O) step on the bare FeN4 site (Fig. 2a). Thus, O=FeN4 or HO-FeN4 would more likely be the true stationary active sites if mononuclear FeN4 is assumed (Fig. 2a and Fig. S4a).

Apart from the FeN4 structure, we have also considered a different coordination situation where two N atoms and two graphite carbon atoms are coordinated to Fe (Fig. 3). 52 The two N atoms can be either on the same side of the Fe center (FeN2-cis) or on the diagonal (FeN2-trans). With decreasing the amount of doped N around Fe, the adsorption of OH is strengthened significantly (by 0.29 and 0.41 eV for FeN2-cis and FeN2-trans, respectively) when comparing with FeN4 (Fig. 2), and hence the fourth step (*OH to H2O) becomes clearly the potential-determining one for ORR. As before, any intermediate in the cycle can be considered as a potential starting point for another cycle of ORR. In contrast to the proceeding results on FeN4, starting from O=FeN2 as active site and forming O=FeN2-OOH is not favorable due to the weak adsorption of OOH when comparing with the continuation of the first cycle and formation of FeN2-OH (Fig. 2bc and Fig. S4bc). As for the adsorption of O on the O=FeN2 sites, since the C atom bound to Fe prefers to form a C-O bond in all possible spin states for O=FeN2-trans and in the high spin state for O=FeN2-cis (Fig. S5), O adsorption would be strong and the active site structure would be distorted. So that the formation of *OOH (*O+ H2O → H + + e -+*OOH) in OER can be largely hindered on O=FeN2 sites. Therefore, neither O=FeN2-trans nor O=FeN2-cis will be discussed further as the potential steady-state active site. Altogether, the present calculations suggest that HO-FeN2 (FeN2-cis and FeN2-trans) would most likely serve as active site for ORR (blue lines in Fig. 2bc). The increased reactivity of OH and OOH intermediates by the presence of the other OH ligand in trans of Fe could be associated to the well-established trans effect in inorganic chemistry [START_REF] Miessler | Inorganic Chemistry[END_REF] . The (HOFeN2)-OH and (HOFeN2)-OOH bonds are elongated and weakened since the trans ligands share the same d orbitals of Fe.

Our computations using the periodic and cluster models of the active sites provide similar results and conclusions that the presence of extra OH ligand and modification of N coordination number around Fe could largely influence the activity (Fig. S6). From the geometrical perspective, the FeNx (x= 2 or 4) doped carbon cluster structures with H atoms termination tend to bend when O, OH and, particularly, the OOH is adsorbed, which potentially due to relative flexibility of the structure and tension on the edge.

Nevertheless, the presence of a hydroxyl ligand HO-FeNx would help to stabilize the intermediate structures in cluster models as the same in periodic models. The adsorption of O on O=FeN2 sites in cluster models would also distort the structure and thus were not considered further as potential active sites. Quantitatively, a weaker adsorption of OH on the cluster models (by 0.2 eV) was observed, due to an approximate description of the graphene electronic structure.

Scaling relations for different proposed active sites.

The comparison of the ORR activity for the seven active sites specified in this study (FeN4, FeN2-cis, FeN2-trans, HO-FeN4, HO-FeN2-cis, HO-FeN2-trans and O=FeN4) can be made more universal by constructing scaling relations for adsorption free energies and volcano plots for the limiting potential. Since the fourth step (*OH to H2O) is potential-determining, the adsorption free energy of OH is a key parameter, and therefore used as a descriptor for the adsorption free energy of all other intermediates. 17 Fig. 3a shows the calculated potential as a function of the adsorption free energy of OH (∆G*OH) for each elementary step in ORR. The fourth step (*OH to H2O) is usually considered as the potential-determining step, however when the adsorption of OH is weakened, the formation of OOH will be the potential-determining step instead. 23 The actual output potential will correspond to the lowest envelope. ∆G*OH increases by following the sequence FeN2-trans < FeN2-cis < FeN4 < HO-FeN4 < HO-FeN2-cis < O=FeN4 < HO-FeN2-trans. In the case of highly stable OH adsorption (left part of Fig. 3a), the fourth step (*OH to H2O in magenta) appears to be potential-determining step

The highest limiting potential and hence the lowest overpotential (0.30 V) is obtained at the crossing point where ∆G*OH is 0.93 eV.

The bare FeN4 site, often proposed as the active site for ORR in the literature, 11,12 is on the left of the diagram, far from the optimum, with a large calculated overpotential of 0.83 V. The bare FeN2 sites adsorb OH even more strongly than FeN4, which is not favorable for ORR. In the preceding section, we have shown that the real steady state active site would bear an extra OH ligand. This -OH ligand significantly weakens the adsorption strength of the extra *OH intermediate and places the sites in the favorable region of the volcano. The best site is thus determined by the present calculations to be HO-FeN2, where Fe is coordinated with two N and two C atoms, and stabilized by an extra OH ligand. In the submission process on the present manuscript, other works appeared in the literature underlining the importance of existing OH ligand 19,20 , either from the dissociation of H2O or from the basic media. Our work suggests that the active sites with extra OH ligand would be formed in situ during the reaction. Also compared with the further OH adsorption on the same side of the already present OH ligand 20 , adsorption on the opposite side of the Fe center, as considered here, is weakened by the trans effect.

The same treatment can be applied to the OER, by considering the reverse steps. The third step (*O to *OOH) is generally considered as the potential-determining step in OER, even though the second step (*OH to *O) could also be limiting. 51 We selected ∆G2 = ∆G*O -∆G*OH as energy descriptor [START_REF] Man | [END_REF] to plot the scaling relation (Fig. 4). ∆G3

shows a good quality linear relation to ∆G2. Since we are now in electrolysis conditions, the diagram should be read differently where the higher potential envelope is the limiting potential. Without any additional ligand, the FeNx sites (x=2 or 4, square symbols) generate large overpotentials (Fig. 4 and Fig. S8). However, with an -OH or -O ligand, the overpotential decreases dramatically (circles and triangle stand for HO-FeNx and O=FeN4, respectively in Fig. 4). The optimal overpotential is 0.32 V for OER.

The promising sites for OER are O=FeN4, HO-FeN2-cis and HO-FeN4. HO-FeN2-cis hence appears as an efficient active both for OER and ORR, leading to small overpotential in each case.

Including the influence of the potential with the SC approach

From the proceeding scaling relations, three promising sites were chosen (FeN4, HO-FeN4 and HO-FeN2-cis) and were subjected to more detailed computation by taking explicitly the potential into account with the SC approach (Fig. S9). In order to understand the influence of the potential on each step, we first looked into the adsorption energy of intermediates. Unlike in the CHE approach, where the potential has no influence on the adsorption energies (left panel of Fig. 5, dash line), the adsorption energy varies as a function of applied potential (left panel of Fig. 5, solid line) in the SC approach. The adsorption energy of formally neutral intermediates is obtained with the same equation as in the CHE approach (equation 9-11 in the ESI), where H2 and H2O are taken as energy reference for the adsorbates at constant charge.

However, the intermediates are considered at different potentials. Generally, the change (in eV per V) is largest for the *O adsorbate on all selected sites followed by the *OOH and *OH species. This is in agreement with *O featuring the highest surface dipole moment and, thus, the most positive workfunction (about 1 V higher than the bare surface, see ESI). Adsorption of O would be destabilized by 0.30, 0.32 and 0.43 eV per V on FeN4, HO-FeN4 and HO-FeN2-cis, respectively, which are far from negligible. To be noted, the crossing point between the SC and CHE adsorption energies occurs at different potential. After determining how the potential affects adsorption energies of each intermediate, the influence of potential on the free energy of the corresponding reaction step can be obtained (Fig. 5, right panels). As formally one electron is transferred in each step, not only the adsorption energies of the intermediates, but also their charge as a function of the potential comes into play since it modulates the effective number of electrons The first step is only linked with OOH adsorption energy, and hence shows a negligible deviation for its reaction energy ∆G1 between CHE and SC approaches for HO-FeN2cis (Fig. 5c, right panel). Larger differences are obtained for FeN4 and HO-FeN4, with a decrease of ∆G1 at equilibrium potential by 0.12 eV with the SC approach. The second (*OOH to *O) and the third (*O to *OH) steps involve *O as an intermediate and hence are generally markedly affected by the potential, with, for example, a decrease of reaction energy at equilibrium potential for the second step (*OOH to *O) by 0.32 eV for HO-FeN2-cis. This is mitigated in cases where OOH and O adsorption energy evolve in the same way (FeN4, HO-FeN4). Since the fourth step (*OH to H2O) only depends on OH adsorption which is weakly potential dependent, it shows a small deviation between CHE and SC approaches. Note that for each step, there is a potential value where the CHE and the SC free energy cross, providing a zone of potentials where results from both approaches coincide.

Based on these graphs, the fourth step (*OH to H2O) is the potential-determining step on all selected sites, since it is the first one to become endothermic upon potential increase. The crossing point at ∆G4 = 0 determines the output potential, and its deviation from the ideal redox-potential gives the overpotential. Since the fourth step (*OH to H2O) is only weakly affected by the applied potential, the calculated overpotentials by the SC approach (0.79, 0.45 and 0.38 V for FeN4, HO-FeN4 and HO-FeN2-cis, respectively) are close to the CHE value (0.81, 0.53 and 0.37 V). However, SC has significant impacts on other steps, which become important if we change from ORR to OER. In this case, the third step (*O to *OOH, the reverse of the second step in ORR) becomes the potential-determining step and a large influence of the potential in the SC approach is shown, especially at high potential relevant to this reaction. The calculated overpotentials for OER with SC (0.59, 0.35 and 0.47 V for FeN4, HO-FeN4 and HO-FeN2-cis, respectively) show a strong decrease compared to CHE values (0.89, 0.69 and 0.65 V) by up to 0.34 V on the HO-FeN4 site (Table S4 and Table S5). One point to note is that the most endothermic step is changing versus potential on HO-FeN2-cis, which is reflected by the change of potential-determining step from the (*O to *OOH) step with CHE approach to the (*OH to *O) step when using the SC approach. In contrast to ORR, the potential-determining step in OER involves the *O intermediate.

The adsorption energy of *O is markedly potential dependent, which explains the important correction brought by the SC approach for OER. In a nutshell, the potential has contrasted influence on the adsorption energy of intermediates, hence on reaction steps and on overpotential.

As noted earlier, 27,55,56 the qualitative impact of SC compared to CHE can be predicted by comparing the workfunction changes for a given elementary reaction step (see Fig. S16). For quantitative predictions, the capacitance (and its change) is, however, required. 28 In particular, we find that even though the FeN4, HO-FeN4 and HO-FeN2 systems are very closely related, the capacitance of FeN4 is almost 50% smaller than for the other two systems and, thus, the impact of the change in workfunction is approximately halved for FeN4 compared to HO-FeN4 or HO-FeN2, indicating that the constant capacitance approximation would be very rough and should not be applied when detailed insights and precise energetic predictions are sought after. 57

Conclusions

Several Fe-N-C active sites for Fe and N co-doped carbon (L-FeNx, x = 2 or 4, L = nothing, O or OH) were considered for a computational study of the catalytic mechanism of ORR and OER. The adsorption free energies of OOH, O and OH are strongly affected by changes in the coordination of the Fe center, either a modification of the number of N and C atoms coordinated to Fe in the equatorial plane, or addition of O or OH ligand in the axial direction. In general, the coordination of an -OH ligand can largely decrease the adsorption of OH on the opposite side of the Fe center, and thus reduce the overpotential, which is crucial for the moderate adsorption of OH in the catalytic cycle and the activity of the FeNx site. From the scaling relations, the optimal potential values were obtained at an applied potential of 0.93 V and 1.55 V for ORR and OER, corresponding to an overpotential of 0.30 V and 0.32 V respectively. The description of Van der Waals interactions and the effect of the solvent do not change the scaling relations and limiting potential volcano plot, but considerably affect the OH

Fig. 1

 1 Fig. 1 Schematic presentation of the four-electron mechanism for the electrochemical ORR. A limited H-terminated cluster model of the site is used and grey, blue, orange, red and white spheres represent C, N, Fe, O and H atoms, respectively.

Fig. 2

 2 Fig. 2 Free energy diagrams for ORR assuming (a) FeN4, (b) FeN2-trans and (c) FeN2-cis as active sites, respectively. Computations include Van der Waals interaction and implicit solvent with the CHE approach at equilibrium potential (1.23 V) is used. Intermediate structures of O=FeNx or HO-FeNx (x=2 or 4) are also considered as active sites to start a second ORR cycle as shown with the red and blue lines and labels, respectively.

Fig. 3

 3 Fig. 3 (a) Potential scaling relations as a function of OH adsorption free energy for ORR using PBE with Van der Waals interaction and solvation effect included. The energies are based on the most stable spin state of each intermediate using periodic models. The free energy of the first, second, third and the fourth step of the four-electron mechanism, namely formation of *OOH, *OOH to *O, *O to *OH and *OH to H2O, are shown in black, red, blue and magenta lines, respectively. The square, the circle and the triangle symbols stand for FeNx, HO-FeNx (x=2 or 4) and O=FeN4, respectively. (b) The geometry of different periodic active site models, which are listed in the sequence of increasing OH adsorption free energy. Grey, blue, orange, red and white spheres represent C, N, Fe, O and H atoms, respectively. The structures are labeled I-VII in (b) and their corresponding performance are shown in (a).

Fig. 4

 4 Fig. 4 Scaling relations for the limiting potential of each step as a function of the difference between O and OH adsorption free energy for OER using PBE with Van der Waals interaction in implicit solvent. The energies are based on the most stable spin state of each intermediate using the periodic models. The square, the circle and the triangle symbols stand for FeNx, HO-FeNx (x=2 or 4) and O=FeN4, respectively (I: FeN2-trans, II: FeN2-cis, III: FeN4, IV: HO-FeN4, V: HO-FeN2-cis, VI: O=FeN4, VII: HO-FeN2-trans).

Fig. 5

 5 Fig. 5 The adsorption energy of intermediate *OOH (black), *O (red) and *OH (blue) as a function of potential (left). The solid and the dash line stand for the SC and CHE results. The free energy with CHE or SC approaches for the various steps of ORR as a function of potential (right). The potential range is from 0.00 to 2.20 V. ∆G1 to ∆G4 represent the free energy of the first, second, third and fourth step of the four-electron mechanism, namely formation of *OOH, *OOH to *O, *O to *OH and *OH to H2O, with black, red, blue and magenta lines, respectively. Considered active sites are (a) FeN4, (b) HO-FeN4 and (c) HO-FeN2-cis.

Acknowledgements

The authors would like to thank China Scholarship Council (CSC) and Tsinghua University (20131089311) for the financial support. This work used computational and storage services associated with the Hoffman2 shared cluster provided by UCLA Institute for Digital Research and Education Group.

adsorption free energy and hence the determination of the best catalyst structure.

The promising active sites with low overpotential (HO-FeN4 and HO-FeN2) were further studied by the SC approach, which includes the influence of the applied potential, solvent and electrolyte. The potential was shown to affect the adsorption energy of intermediates, especially on *O with a destabilization of up to 0.43 eV/V.

Overall, the influence of an explicit treatment of the potential appears minor for ORR, but a large influence was observed in OER with an overpotential reduced by up to 0.34 V. Thus, the explicit treatment of the potential does have an influence on the free energy steps but the impact on the limiting potential depends on the reaction.

In conclusion, this work provides key insight on the optimal structure of the catalytic active site for ORR and OER on Fe-N-C catalysts from a theoretical perspective, which is important for the design of new efficient catalyst. Our comparison between the CHE and SC approaches illustrates the difference between these two approaches and offer instructions for a method choice in future studies.

Conflicts of interest

There are no conflicts to declare.


Sarah Blanck

Sophie Loehlé

Stephan N Steinmann

Carine Michel
email: carine.michel@ens-lyon.fr

Adhesion of lubricant on aluminium through adsorption of additive head-groups on γ-alumina: a DFT study

To improve lubricant formulation, a better understanding of the wettability properties and hence of the interactions between the surface and the lubricant additives is a key parameter. Herein, we use density functional theory calculations to characterize 32 typical head groups, classified into four categories: protic, aprotic, aromatic and phosphorus-containing molecules. By comparing their adsorption energies on γ-alumina, used as a surface model of aluminium sheets, and their solvation energies in a model lubricant base oil, we found that the solvation energy was not a discriminant parameter while the adsorption energy was critical. Phosphates and carboxylic acids are the most strongly adsorbed, and thus more likely to yield to improved wettability properties of the lubricant through film formation.

Introduction

Lubricants are widely used formulated products, with a global market size valued over USD 125 billion in 2018. The major applications occur in the automotive, the marine, and the industrial (metal working and textile) sectors. They are usually composed of a base oil which can be made of either mineral, vegetal or synthetic oils depending on the usage [START_REF] Ayel | Lubrifiants -Constitution[END_REF], and a package of additives used to enhance the properties of interest of the base oil or add new properties to the lubricant, generally leading to a rather complex composition. A wide range of additives for various purposes can be found in lubricants. Empirically, they are classified into different groups depending on their function in the lubricant, e.g., surfactants, friction modifiers, anti-wear and extreme pressure additives and corrosion inhibitors [START_REF] Rudnick | Lubricant additives: chemistry and applications[END_REF][START_REF] Kajdas | Additives for metalworking lubricants-a review[END_REF][START_REF] Minami | Concept of molecular design towards additive technology for advanced lubricants[END_REF]. In this study, we focus on lubricants to be used in metalworking of aluminium sheets, aluminium being for instance more and more used as a cost-effective alternative to steel in the automotive industry where weight reduction is key to both fuel-operated and electric vehicles. In this field, friction modifiers are the most studied additives, as they directly influence the tribological properties of lubricants [START_REF] Righi | A comparative study on the functionality of S-and P-based lubricant additives by combined first principles and experimental analysis[END_REF][START_REF] Lingala | Effect of lubricant additives on the tribological behavior of aluminum alloy against steel[END_REF][START_REF] Wan | Tribological behavior of lubricating oil additives in lubricated aluminumon-steel contact[END_REF]. On the contrary, only few studies have discussed the wettability properties of lubricants and more generally the influence of additives on the wetting of a given surface [START_REF] Zhou | Dynamic wetting of rolling oil on aluminum surfaces[END_REF][START_REF] Arieta | Experimental technique to determine the wettability of surfaces to lubricating oils[END_REF], which is the main focus of the present study.

Even if lubricants are used for centuries [START_REF] Rudnick | Lubricant additives: chemistry and applications[END_REF], their development is still ongoing and the lubricants composition is constantly evolving. Indeed, changes in formulations are triggered mostly by modifications of legislations and by the constant search for better performances. The changes imposed by more stringent legislation can lead to significant effects on the behavior of the lubricant towards the surface, e.g., it can deteriorate the wettability property of the lubricant, i.e., its ability to quickly spread over the surface to achieve its optimal tribological performance [START_REF] Vakis | Modeling and simulation in tribology across scales: An overview[END_REF][START_REF] Nishi | Influence of unforced dewetting and enforced wetting on real contact formation and friction behavior between rubber hemisphere and glass plate during contacting and sliding processes[END_REF]. Therefore, it is important to be able to understand and to predict the effects of the different components of the lubricant and their link with the wettability properties. Indeed, surface engineering to achieve the desired wettability is an active field of research. These modifications can be obtained by either chemical of physical processes [START_REF] Wong | Dependence of Macroscopic Wetting on Nanoscopic Surface Textures[END_REF][START_REF] Anastasiadis | Development of Functional Polymer Surfaces with Controlled Wettability[END_REF][START_REF] Mundo | Influence of Chemistry on Wetting Dynamics of Nanotextured Hydrophobic Surfaces[END_REF] and can even lead to superhydrophobic or superhydrophilic properties [START_REF] Li | What do we need for a superhydrophobic surface? A review on the recent progress in the preparation of superhydrophobic surfaces[END_REF][START_REF] Feng | Design and Creation of Superwetting/Antiwetting Surfaces[END_REF][START_REF] Kietzig | Patterned Superhydrophobic Metallic Surfaces[END_REF][START_REF] Kumar | Development of durable self-cleaning superhydrophobic coatings for aluminium surfaces via chemical etching method[END_REF]. However, in the case of lubrication, the surface is often predetermined and cannot be cost-efficiently modified before application of the lubricant to increase its wettability. Hence, the lubricant formulations need to be adapted to get an acceptable wettability.

A large number of experimental methods have been developed to measure the contact angle of a liquid on a surface [START_REF] Yuan | Contact Angle and Wetting Properties[END_REF] but in some case, contact angles are not sufficient to describe the difference in wetting observed between two liquids [START_REF] Gim | Multiscale Simulation Method for Quantitative Prediction of Surface Wettability at the Atomistic Level[END_REF][START_REF] Kalin | The correlation between the surface energy, the contact angle and the spreading parameter, and their relevance for the wetting behaviour of DLC with lubricating oils[END_REF] and other techniques are needed to be able to predict the difference in the wettability properties. Theoretical models have been developed by Blake et al. to describe the dynamics of wetting of a liquid drop on a surface [START_REF] Blake | The influence of solid-liquid interactions on dynamic wetting[END_REF][START_REF] Bertrand | Influence of solid-liquid interactions on dynamic wetting: a molecular dynamics study[END_REF]. They showed that solid-liquid interactions affect the spreading of the liquid on the surface in two opposite directions: on the one hand, the interaction between the liquid and the solid constitutes the driving force for spreading and, on the other hand, too strong interactions result in a resistance to wetting. Therefore, very much like a volcano curve in catalysis [START_REF] Nørskov | The nature of the active site in heterogeneous metal catalysis[END_REF], a system dependent optimal value of solid-liquid interactions, corresponding to the maximal wetting of surface, should be aimed at. Furthermore, to affect the wettability, the additives should have the tendency to go to the solid-liquid interface rather than to stay in the liquid bulk, i.e., their solvation free energy should not be more important than the adsorption free energy.

Based on these considerations, we set out to theoretically compare the adsorption and the solvation energies of different types of additive head groups, aiming to classify them according to their expected influence on the spreading of the lubricant on the considered surface. DFT is the method of choice for theoretically determining the adsorption free energy of additives on surfaces. In contrast to force fields, this first principles method is applicable to any head group imaginable, while its computational cost is acceptable when using single crystal surfaces and model head-groups of up tens of atoms. Its accuracy for adsorption energies has previously been assessed by comparison with accurate experimental adsorption energies and found to be satisfactory when dispersion interactions are properly taken into account [START_REF] Gautier | Molecular adsorption at Pt(111). How accurate are DFT functionals?[END_REF]. Furthermore, previous tribological studies based on density functional theory (DFT) [START_REF] Righi | A comparative study on the functionality of S-and P-based lubricant additives by combined first principles and experimental analysis[END_REF][START_REF] Gattinoni | Adsorption of Surfactants on α-Fe 2O3 (0001): A Density Functional Theory Study[END_REF] of various lubricant additives on different surfaces obtained results in good agreement with experimental data. The major interaction of the lubricant additives with the surface comes from the adsorption of polar head groups, while the apolar tail ensures solubility and smooth interactions with the base oil to achieve low friction. For this reason, we choose to limit our study to the interactions between the additives head groups and the surface, which is equivalent to consider that the carbon chains of the additives are all the same and play the same role for all additives. The studied polar head groups cover a wide range of additives including protic and aprotic molecules, functionalized aromatics and phosphorus containing compounds. The exhaustive list is given in Table 1 to Table 4.

Regarding the surface model, since aluminium oxidizes immediately when exposed to air, we choose to model the surface by aluminium oxide. Alumina can exhibit different types of crystallographic structures depending on the temperature [START_REF] Wefers | Oxides and hydroxides of aluminum[END_REF][START_REF] Santos | Standard transition aluminas. Electron microscopy studies[END_REF]. It has been shown that the oxidation of aluminium in air leads to the formation of an amorphous alumina layer of 10 to 20 Å of thickness [START_REF] Jeurgens | Structure and morphology of aluminiumoxide films formed by thermal oxidation of aluminium[END_REF][START_REF] Olefjord | Surface analysis of oxidized aluminium. 2. Oxidation of aluminium in dry and humid atmosphere studied by ESCA, SEM, SAM and EDX[END_REF]. However, aluminium sheets are fabricated by a rolling process at temperatures between 420 and 530°C [START_REF] Deneuville | Mise en forme de l'aluminium -laminage[END_REF]. Under these conditions, the amorphous alumina layer is transformed into crystalline g-alumina [START_REF] Jeurgens | Structure and morphology of aluminiumoxide films formed by thermal oxidation of aluminium[END_REF]. Therefore, g-Al2O3 has been chosen in order to model the oxide layer present at the top of the studied aluminium sheets.

After defining the surface of g-alumina which is the most suitable to model the considered system, we will discuss the values of adsorption energies obtained for the different studied head groups of additives.

Then, we will study the solvation energies obtained for the same molecules. And finally, we will be able to compare the adsorption and solvation energies in order to assess which types of additives seem to be the most appropriate to improve the wettability of the lubricant on the considered surface.

2 Models and method

Computational details

All DFT based computations are carried out using CP2K [START_REF] Hutter | cp2k: atomistic simulations of condensed matter systems[END_REF]. The exchange correlation energies were computed with the generalized gradient approximation of Perdew, Burke and Ernzerhof (PBE) [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] with the DFT-D3 dispersion correction developed by Grimme et al. including C9 terms [START_REF] Grimme | A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu[END_REF][START_REF] Grimme | Effect of the damping function in dispersion corrected density functional theory[END_REF]. A mixed Gaussian and plane waves (GPW) [START_REF] Lippert | A hybrid Gaussian and plane wave density functional scheme[END_REF][START_REF] Vandevondele | Quickstep: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach[END_REF] approach has been adopted. A double-z basis set as provided by the Molopt basis set library is used to describe valence electrons, while core electrons were treated with the Goedecker, Teter and Hutter (GTH) pseudopotentials [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][START_REF] Vandevondele | Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases[END_REF]. A 400 Ry cutoff energy has been chosen for the plane waves basis set. During the self-consistent field process, the energies are converged to 1•10 -6 Ha.

The geometries are optimized at the G-point to forces below 4.5•10 -4 Ha/bohr. The total energies are then evaluated using a single point energy calculation with a 4x3x1 Monkhorst-Pack [START_REF] Monkhorst | Special points for Brillouin-zone integrations[END_REF] k-point grid for a p(1x1) unit cell and adapted to keep the k-point density constant for larger cells.

The g-alumina (100) surface is modeled by a periodic slab, with a finite thickness to limit the overall computational cost. This slab is non-symmetric, with one of the exposed facet being frozen in the alumina bulk position to mimic the greater thickness of the alumina found on the aluminium sheet, while the opposite facet is free to adapt upon adsorption of any molecules, to mimic the exposed surface of the aluminium sheet. In our case, we chose a p(1x1), p(2x1) or p(2x2) unit cell depending on the size of the adsorbed molecule, with a thickness of 14 atomic planes (two layers), the bottom 8 of which were kept frozen (Figure 1b). A vacuum thickness of 20 Å was set between two periodically repeated slabs to allow the adsorption of molecules on the upper relaxed side of the slab without interacting with the bottom of the periodic image of this slab. The adsorption energy is defined as:

𝐸 "#$ = 𝐸 &'()*+()@-( . / 0 -2𝐸 -( . / 0 + 𝐸 &'()*+() 4
The solvation energy evaluations are performed with Gaussian09 [START_REF] Frisch | Gaussian 09[END_REF] at the DFT level with the PBE functional and a triple-z basis set Def2-TZVP [START_REF] Weigend | Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy[END_REF][START_REF] Weigend | Accurate Coulomb-fitting basis sets for H to Rn[END_REF], using the Polarizable Continuum Model (PCM). The PCM method used here is the Integral Equation Formalism PCM (IEFPCM) [START_REF] Mennucci | Evaluation of Solvent Effects in Isotropic and Anisotropic Dielectrics and in Ionic Solutions with a Unified Integral Equation Method: Theoretical Bases, Computational Implementation, and Numerical Applications[END_REF][START_REF] Cancès | A new integral equation formalism for the polarizable continuum model: Theoretical background and applications to isotropic and anisotropic dielectrics[END_REF], which is the default one in Gaussian09. In this method, the solvent is described thanks to macroscopic properties such as its permittivity [START_REF] Tomasi | Quantum Mechanical Continuum Solvation Models[END_REF]. The solute is treated by quantum mechanics the same way as in vacuum, and the interactions between the solvent and the solute are added by an interaction potential 𝑉 6 7 . Therefore, the vacuum Hamiltonian 𝐻 9 : is replaced by an effective Hamiltonian 𝐻 9 );; given by: 𝐻 9 );; = 𝐻 9 : + 𝑉 6 7 [START_REF] Mennucci | Evaluation of Solvent Effects in Isotropic and Anisotropic Dielectrics and in Ionic Solutions with a Unified Integral Equation Method: Theoretical Bases, Computational Implementation, and Numerical Applications[END_REF].

g-alumina

g-alumina is a commonly used material, especially as a catalytic support for heterogeneous catalysts, explaining the abundance of theoretical and experimental studies involving g-alumina [START_REF] Réocreux | Structuration and Dynamics of Interfacial Liquid Water at Hydrated γ-Alumina Determined by ab Initio Molecular Simulations: Implications for Nanoparticle Stability[END_REF][START_REF] Ngouana-Wakou | An Atomistic Description of the γ-Alumina/Water Interface Revealed by Ab Initio Molecular Dynamics[END_REF].

Nevertheless, since it is not perfectly crystalline in reality, its structure is not well known and different types of models can be found in the literature [START_REF] Krokidis | Theoretical Study of the Dehydration Process of Boehmite to γ-Alumina[END_REF][START_REF] Gutiérrez | Theoretical structure determination of γ-Al2O3[END_REF][START_REF] Ferreira | Direct comparison between two structural models by DFT calculations[END_REF]. In agreement with our previous studies [START_REF] Réocreux | Structuration and Dynamics of Interfacial Liquid Water at Hydrated γ-Alumina Determined by ab Initio Molecular Simulations: Implications for Nanoparticle Stability[END_REF][START_REF] Réocreux | Reactivity of shapecontrolled crystals and metadynamics simulations locate the weak spots of alumina in water[END_REF],

we chose to work with the non-spinel model of g-alumina proposed by Krokidis et al. [START_REF] Krokidis | Theoretical Study of the Dehydration Process of Boehmite to γ-Alumina[END_REF][START_REF] Digne | Use of DFT to achieve a rational understanding of acid-basic properties of γ-alumina surfaces[END_REF] .

g-alumina mainly exposes three different facets: (100), ( 110) and (111). The (100) surface is the one having the lowest surface energy [START_REF] Digne | Use of DFT to achieve a rational understanding of acid-basic properties of γ-alumina surfaces[END_REF], and therefore the most probably exposed, and has for this reason been adopted in this study.

The p(1x1) unit cell of this surface exhibits four aluminium AlV atoms (Figure 1a), where the subscript V indicates that the atom is coordinated to five oxygen atoms. Two of the aluminium atoms are equivalent (AlVb and AlVb'). Three non-equivalent oxygen atoms can also be seen on the surface, one of them being coordinated to four aluminium atoms (O4 which is equivalent to O4'), whereas the others (O3a, O3a', O3b and O3b') are coordinated to three aluminium atoms only.

Al Vb Al Vc Al Vb' O 4' O 3b' O 3b Al Va O 4 O 3a' O 3a (a)
Flexible atomic planes

Frozen atomic planes

Set of studied molecules

Different kinds of molecules of interest have been studied. These molecules were chosen because they model the polar head group of various additives of fully formulated lubricants. These functional groups thus probe the major adsorption energy contributions of the additives. The chosen molecules have been classified into four categories depending on their characteristics.

Molecules derived from fatty acids containing a polar group, such as a carboxylic acid, alcohols or esters, are often used as friction modifiers in lubricant formulations, whereas amines exhibit good antioxidant and corrosion inhibitor properties [START_REF] Rudnick | Lubricant additives: chemistry and applications[END_REF]. We have, however, chosen to classify them according to basic physical chemical properties. Hence, the two first families gather together respectively protic (Table 1)

and aprotic (Table 2) small molecules with various functional groups. The third group (Table 3) is composed of aromatic molecules. This category serves to determine the influence of the aromatic ring on functional groups of particular interest. Indeed, the steric hindrance induced by the presence of an aromatic group can modify the adsorption of a molecule. The last category gathers phosphorous compounds (Table 4). As they can be used as anti-wear, extreme pressure additives or antioxidants, these compounds are key lubricant additives [START_REF] Rudnick | Lubricant additives: chemistry and applications[END_REF]. Both our study and literature [START_REF] Digne | Use of DFT to achieve a rational understanding of acid-basic properties of γ-alumina surfaces[END_REF][START_REF] Digne | Hydroxyl Groups on γ-Alumina Surfaces: A DFT Study[END_REF][START_REF] Zuo | Effect of surface hydroxyls on DME and methanol adsorption over γ-Al2O3 (hkl) surfaces and solvent effects: a density functional theory study[END_REF] identified the AlVa aluminium atom as the most reactive one, and therefore as the preferential adsorption site. This could have been predicted from the Lewis acidities of the different aluminium atoms of the surface [START_REF] Digne | Use of DFT to achieve a rational understanding of acid-basic properties of γ-alumina surfaces[END_REF]. Indeed, Digne et al. have shown that AlVa is the most acidic aluminium atom of the surface. Additionally, surface oxygen atoms can also be involved in the adsorption process, acting as Bronsted basic sites. The adsorption energy of each molecule is reported in Table S1.

Protic molecules

For protic molecules, the Bronsted acidity of the molecule is a key parameter to understand the strength of the adsorption.

The water molecule is adsorbed on the (100) g-alumina surface with an adsorption energy of -0.97 eV. This adsorption is dissociative: the oxygen atom of the water molecule is adsorbed on the AlVa aluminium atom and the hydrogen atom is adsorbed on the O3b oxygen atom of the surface (Figure 2a). This is consistent with previous studies by Digne et al. which obtained the same adsorption mode with a similar energy of -1.08 eV [START_REF] Digne | Use of DFT to achieve a rational understanding of acid-basic properties of γ-alumina surfaces[END_REF][START_REF] Digne | Hydroxyl Groups on γ-Alumina Surfaces: A DFT Study[END_REF]. The slight difference in energy observed between this study and the literature can be explained by variations in the computational method.

Given the adsorption mode of water, one could think that the adsorption of methanol would occur analogously. However, as previously published by Zuo et al. [START_REF] Zuo | Effect of surface hydroxyls on DME and methanol adsorption over γ-Al2O3 (hkl) surfaces and solvent effects: a density functional theory study[END_REF], we find that the adsorption of methanol is following a non-dissociative mode. This can be rationalized by the fact that methanol has a higher pKa (15.2) than water. The adsorption mode corresponds to the oxygen atom of the alcohol adsorbed on AlVa with an energy of -1.07 eV. This adsorption energy is slightly stronger than the one obtained by Zuo et al. (-0.88 eV) [START_REF] Zuo | Effect of surface hydroxyls on DME and methanol adsorption over γ-Al2O3 (hkl) surfaces and solvent effects: a density functional theory study[END_REF], because we account for dispersion interactions. With -1.00 eV, the dissociative adsorption of methanol is nevertheless not significantly less favorable than the nondissociative one. Similarly, methanethiol is adsorbed in a non-dissociative way. Despite larger dispersion interactions, methanethiol is even 0.1 eV less strongly adsorbed than water (-0.89 eV). The selected PAG monomer on the contrary, is adsorbed the same way as water in a dissociative mode on AlVa and O3b and sticks more to the surface (-1.29 eV) because of dispersion interactions.

Let us now move to the family of amines, which interact with the g-alumina surface through their nitrogen atom. NH3 is not dissociated during the adsorption, but nevertheless reaches an adsorption energy of -1.24 eV. The non-dissociative character of this adsorption was expected because of the high pKa value of ammonia, which is equal to 23. The nitrogen atom of the molecule is still adsorbed on the AlVa aluminium atom. The same occurs for primary, secondary and tertiary amines. The adsorption is stronger for the primary amine than for NH3 and even more for the secondary amine (-1.43 and -1.55 eV respectively). However, the adsorption energy of the tertiary amine (-1.29 eV) is decreasing due to the steric effects. This steric hindrance is also reflected in the geometry: the bond between the nitrogen atom of the tertiary amine and AlVa is much longer (2.21 Å) than the one observed with primary or secondary amine (2.05 Å). Now, let us consider the more complex case of ethoxylated amine. This polydentate molecule preferentially adsorbs in a bidentate mode, involving either the nitrogen and a hydroxyl group, or the two hydroxyl groups (Figure 2b and Figure 2c). In the most stable one (-1.78 eV), the N is not involved: the oxygen atom of the first hydroxyl group and the hydrogen atom of the second one interact with AlVa and O3a, respectively (Figure 2b). This is quite surprising as the adsorption of a tertiary amine is stronger as the one of an alcohol. But, the strongest adsorption obtained by bonding the nitrogen atom to the surface has an energy of -1.41 eV, demonstrating that the adsorption of complex multidentate molecules on this alumina surface cannot simply be derived from their sub-groups. The corresponding bi-dentate structure is shown in Figure 2c.

The carboxylic acid also shows a bidentate adsorption, its two oxygen atoms adsorbed on the AlVa and AlVc sites (Figure 2d). Moreover, the OH bond is dissociated and the hydrogen atom is adsorbed on the O3b' atom of the surface. This adsorption is very favorable, having an adsorption energy of -1.96 eV.

Similarly, two nitrogen atoms of triazole are adsorbed on AlVa and AlVc and the dissociated hydrogen atom on O3b, yielding an adsorption energy of -1.55 eV.

From these results, we can conclude that the adsorption of protic molecules on the (100) g-alumina surface leads in most cases to a dissociative adsorption for acidic molecules with energies up to -1.96 eV and that soft Lewis bases such as thiols only interact weakly with the surface. As explained previously, some of the studied protic molecules are exhibiting dissociative adsorptions. This also occurs for some aromatic and phosphorous compounds which are also protic but will be discussed in the following sections. For these molecules, the adsorption energy strongly correlates with the pKa value in water (Figure 3). The pKa indicates the strength of an acid: the lower the value of pKa, the strongest the acid is. Therefore, the lower the pKa value, the more the molecule will have the tendency to dissociate and the stronger the adsorption energy will be.

Aprotic molecules

The aprotic molecules we selected for this study are Lewis bases and will therefore interact with their electron pair on AlVa. The majority of the aprotic molecules under consideration are adsorbed in a monodentate way, but some of them show a bidentate adsorption on the alumina surface.

The adsorption of dimethylether occurs via its oxygen atom and features an energy of -1.15 eV. This is consistent with a previous study [START_REF] Zuo | Effect of surface hydroxyls on DME and methanol adsorption over γ-Al2O3 (hkl) surfaces and solvent effects: a density functional theory study[END_REF] that obtained an adsorption energy of -0.87 eV on the (100) surface of g-alumina, with the energy difference coming mainly from the dispersion contribution. The carbonyl oxygen of acetone is a weaker Lewis-base and the interaction energy is indeed slightly weaker (-1.00 eV) compared to the ether. For the ester, either the carbonyl or the other oxygen atom could be the anchoring site. Considering the adsorption of ketone and ether, one could first think that the adsorption of ester would occur on the ether-like oxygen. However, the carbonyl oxygen benefits from electron delocalization in minor Lewis structures and thus becomes the most basic oxygen atom. Hence, the adsorption via the carbonyl oxygen atom reaches -1.08 eV, while the other one leads to only -0.84 eV. Adsoprtion energy (eV) pKa Dimethylsulfide is adsorbed exactly the same way as the ether but slightly more strongly: the sulfur atom adsorbs on AlVa yielding an adsorption energy of -1.28 eV. The greater strength of adsorption results from the combination of a lower destabilizing deformation of the surface with greater stabilizing dispersion interactions.

As the nitrogen atom of amides is not basic, the adsorption of N,N-dimethylacetamide is analogous to the one of ester and occurs via the oxygen of the carbonyl group. It leads to a relatively strong adsorption energy of -1.40 eV presumably due to a more important delocalization of electrons from the nitrogen to the carbonyl oxygen, which is also responsible for the more rigid planarity of amide bonds compared to esters [START_REF] Macarthur | Deviations from Planarity of the Peptide Bond in Peptides and Proteins[END_REF].

1-methyltriazole adsorbs the same way than triazole, the two nitrogen atoms not linked to the methyl group are adsorbed on the AlVa and AlVc atoms of the surface, the only difference being the fact that 1methyltriazole does not dissociate as it is aprotic. For this reason, its adsorption is weaker than the one observed for triazole as it has an adsorption energy of -1.19 eV, compared to -1.55 eV for the triazole.

For the 1-methyl-1,2,4-triazole, contrary to what has been observed for the triazole and 1-methyltriazole molecules, the bidentate adsorption is not stable. Indeed, while performing a geometry optimization calculation with an initial bidentate adsorption of the molecule on the surface, the molecule moves away from the surface and adsorbs then again in a monodentate way. This can be expected as it is not favorable to adsorb the methylated nitrogen due to steric effects, and the two other nitrogen atoms are too far from each other to allow a bidentate adsorption. Thus, the most favorable geometry is obtained for a monodentate adsorption of the N4 nitrogen atom on AlVa with an energy of -1.44 eV. The 1,2dimethyl-4,5-dihydroimidazole molecule is adsorbed on the surface by the less substituted nitrogen atom. This can easily be predicted considering steric effects. The resulting adsorption energy (-1.40 eV) is quite strong.

From the results obtained it can be seen that the aprotic molecules studied here exhibit a smaller range of adsorption energies than the one observed with protic molecules and that the adsorption mode is driven by Lewis acid/base pairing.

Aromatic molecules

Aromatic groups are present in many lubricant additives. Therefore, it was essential to determine the influence of aromatic rings on functional groups of particular interest.

Benzene is the simplest aromatic compound and it is the common structure of all aromatics tested here.

But this molecule does not chemisorb on the (100) g-alumina surface and stays flat above the surface, at a distance of about 3 Å, which is characteristic of a physisorption. On the contrary, all other aromatic molecules studied here adsorb on the surface, thanks to their substituents.

The adsorption of phenol occurs in analogy to the dissociative adsorption of methanol: its OH group dissociates and the oxygen atom is adsorbed on the AlVa aluminium as the hydrogen atom is transferred to the O3b' oxygen atom. The adsorption of phenol is, however, stronger (-1.38 eV) than the one obtained for methanol. This can be explained by the fact that the dispersion contribution to adsorption energy is greater for the phenol than for the methanol (-0.54 and -0.27 eV respectively). Moreover, phenol has a higher tendency to dissociate than methanol due to its lower pKa value (10.0 for phenol compared to 15.2 for methanol), which explains the difference in dissociation behavior compared to methanol. The 18 same occurs for 4-methylphenol, with exactly the same energy than the one observed for the phenol.

Butylated hydroxytoluene, on the contrary, does not dissociate nor adsorb on the surface. Indeed, as benzene, it remains at a distance of 3.5 Å above the surface, due to the steric effects induced by the tertbutyl groups.

The adsorption mode of diphenylamine is equivalent to the one of the secondary amine: the nitrogen atom adsorbs on AlVa, with an energy of -1.28 eV. The adsorption energy is however weaker, which can be rationalized by the combined effect of electron delocalization of the nitrogen lone-pair with the aromatic cycle and the larger steric hinderance of the phenyl groups compared to the methyl groups.

Similarly, the adsorption of anisole follows the one of ether, but is slightly weaker (-1.02 eV).

In a nutshell, the adsorption of aromatic molecules on alumina is driven by the nature of the substituent.

When protic, it is driven by the pKa if the substituent can sterically access to the surface. When exhibiting Lewis base properties, it will more or less strongly bind, in a subtle compromise between the Lewis acid/base interaction, the dispersion, steric hindrance, etc.

Phosphorus-containing molecules

Phosphorous compounds exhibit a wide range of adsorption energies going from -0.70 to -2.41 eV. They globally have higher adsorption energies than the other molecules studied here, as ether, ketone or alcohol for example. This is due to the fact that phosphorus is less electronegative than carbon, so the P-O bond is more polarized than the C-O bond and, for this reason, the oxygen atoms of phosphorous compounds are stronger Lewis bases.

Methylphosphate and dimethylphosphate are adsorbed on the surface alike (Figure 4a). For each of them, one OH group dissociates, creating two equivalent oxygen atoms on the molecule. These two oxygen atoms are adsorbed on the AlVa and AlVc atoms of the surface, while the dissociated hydrogen binds to the O3b' oxygen atom. As a result, this bidentate adsorption mode is very stable, with an adsorption energy of -2.41 eV and -2.34 eV for methylphosphate and dimethylphosphate respectively, which is very high compared to all other molecules studied here. The strength of the adsorption lies not only in the formation of two Al-O bonds but also in the low pKa values (1.54 and 1.29 for methylphosphate and dimethylphosphate respectively), which yields to a dissociative adsorption.

Comparatively, as there is no hydroxyl group in the trimethylphosphate molecule, its adsorption is weaker (-1.61 eV) and occurs differently: only its non-substituted oxygen is strongly interacting with AlVa.

O,O,O-trimethylthiophosphate is very similar to trimethylphosphate, the only difference is the replacement of the non-methylated oxygen by a sulfur atom. For this reason, one could think that the two molecules adsorb in a similar way through the non-substituted atom, sulfur and oxygen, respectively. This adsorption mode is found to lie at -0.95 eV for the O,O,O-trimethylthiophosphate, much weaker than for the fully oxygenated analog (-1.61 eV). However, this is not the most favorable adsorption for the sulfurized derivative. When one of the oxygen atoms is adsorbed on the AlVa aluminium atom and a second one on AlVc (Figure 4b), an adsorption energy of -1.10 eV is reached. Then, the adsorption of O,O,O-trimethylthiophosphate resembles much more the ether adsorption than the other phosphorous compounds as the adsorption is here performed via two ether-like oxygen atoms.

Going further in the sulfurization, let us now consider O,O,S-trimethyldithiophosphate. It presents a bidentate adsorption of -1.23 eV of the methylated sulfur and of one of the oxygen atom on the AlVa and AlVc atoms respectively (Figure 4c). Indeed, despite the fact that this molecule is quite similar to O,O,O-trimethylthiophosphate as the only difference is an oxygen atom replaced by a sulfur, they exhibit different adsorptions. This can be intuited regarding the difference of adsorption strength of ether and dimethylsulfide: the latter is more strongly adsorbed than ether, which explains the fact that in the case of O,O,S-trimethyldithiophosphate, the methylated sulfur sticks to the surface instead of the methylated oxygen.

In summary, since several adsorption modes are in competition for these compounds, the best mode of adsorption cannot be easily predicted and mixtures are likely relevant under experimental conditions.

The adsorption of dimethylphosphite is similar to the one of trimethylphosphate discussed before, i.e., it is the non-substituted oxygen that sticks to the AlVa aluminium atom of the surface with an energy of -1.53 eV, slightly lower than the one observed for trimethylphosphate. The trimethylphosphite molecule is quite different from all other phosphorous compounds studied. This molecule shows a bidentate adsorption with two oxygen atoms sticking to the AlVa and AlVc aluminium atoms of the surface (Figure 4d) by total adsorption energy of -0.70 eV, which is very low, compared to the other molecules.

Solvation

Prediction of solvation effects at the solid/liquid interface is very challenging due to the quasi-absence of reliable experimental reference data. Hence, the use of implicit solvent models [START_REF] Fisicaro | Soft-Sphere Continuum Solvation in Electronic-Structure Calculations[END_REF] or including a few explicit solvent molecules [START_REF] Desai | A Periodic Density Functional Theory Analysis of the Effect of Water Molecules on Deprotonation of Acetic Acid over Pd(111)[END_REF] to capture the main influence of the solvent is a common strategy [START_REF] Schweitzer | Can microsolvation effects be estimated from vacuum computations? A case-study of alcohol decomposition at the H2O/Pt(111) interface[END_REF][START_REF] Michel | Role of water in metal catalyst performance for ketone hydrogenation: a joint experimental and theoretical study on levulinic acid conversion into gamma-valerolactone[END_REF].

Alternatively, molecular mechanics based simulations also allow to assess the solvation energy at the interface [START_REF] Steinmann | Solvation free energies for periodic surfaces: comparison of implicit and explicit solvation models[END_REF], with the issue of identifying a force field which is accurate for the solid liquid interface [START_REF] Steinmann | Force Field for Water over Pt(111): Development, Assessment, and Comparison[END_REF]. Since the base oils are rather apolar, we use here an implicit model to compute solvation energies (see 2.1 for a complete description). This approach was shown to deliver accurate values [START_REF] Cossi | Energies, structures, and electronic properties of molecules in solution with the C-PCM solvation model[END_REF]. We herein make also a drastic approximation: the bulk solution solvation energies are assumed to be a reliable proxy to study solvation effects on the adsorption energy, implicitly assuming that the entire solvation energy is lost upon adsorption.

(c) (d) (b) (a)
The solvation energies of the different molecules studied were calculated in two different solvents: nhexadecane and pentyl ethanoate. These two solvents were chosen to model two different kinds of base oil used in lubricant formulations: saturated hydrocarbons (group II; modeled by hexadecane) and esters (modeled by pentyl ethanoate). Comparing these two solvents, n-hexadecane has the least solvating power, but trends between different molecules are the same for the two chosen solvents.

Quantitatively, the solvation energies in pentyl ethanoate are 81% higher compared to n-hexadecane.

For this reason and as group II base oil are more commonly used in lubricant formulation than ester base oils such as ethylhexyl palmitate, only the solvation in hexadecane will be discussed here.

Since we are assessing the solvation energy of small polar head groups in an apolar solvent, the solvation energies are quite small. Furthermore, they are quite similar, ranging from -0.02 for the tertiary amine to -0.15 eV for methylphosphate. General tendencies can nevertheless be seen from these results. Aromatic compounds are less solvated than the other types of molecules, with solvation energies of about -0.07 eV. The phosphorous compounds, on the contrary, are on average more solvated than the other molecules. In general, the more substituted a molecule is, the less it is solvated. This phenomenon can be observed for different families of molecules studied here, among which we can find the phosphates and the amines, but also the phosphites, the triazoles… As an example, let us consider the phosphates which have quite strong solvation energies. Methylphosphate is better solvated (-0.15 eV) than dimethylphosphate (-0.13 eV) while trimethylphosphate is even less solvated (-0.10 eV). For the amines, the primary amine has a solvation energy of -0.05 eV which is weaker than the one of ammonia (-0.06 eV) but stronger than the one calculated for the secondary amine (-0.03 eV) which is also stronger than the energy obtained for the solvation of the tertiary amine (-0.02 eV).

Even more generally, the results show that more polar molecules are better solvated. This result can be quite surprising on first sight given that we are discussing a non-polar solvent. However, it can be explained by the fact that we are looking only at the solvation energy and not considering the problem of miscibility. Indeed, if we would like to assess the miscibility (or solubility), two different phenomena have to be considered. The first one is that polar molecules are always more solvated than non-polar ones, independently of the characteristics of the solvent. The second phenomenon is the competition between molecule-molecule and molecule-solvent interactions. In apolar solvents, molecule-molecule interactions are more favorable for polar molecules than molecule-solvent interactions, which is the reason for demixing in apolar solvents. Demixing and low solubility is, however, not a real issue in our case since the real additives contain a long hydrophobic tail that ensures solubility in the apolar base oil and we assume that this tail will influence all head groups in roughly the same way.

Comparison between adsorption and solvation energies

After having studied separately both the adsorption and the solvation energies of the different polar head groups of additives, we now compare the adsorption with the solvation energies. As solid-liquid interactions are favoring the wetting of the surface [START_REF] Blake | The influence of solid-liquid interactions on dynamic wetting[END_REF], we are looking for additives which have the tendency to go at the interface between the solid surface and the lubricant rather than to stay in the bulk liquid, in order to get better wettability properties, i.e. a better spreading of the lubricant over the surface. For this reason, additives with a high adsorption energy and a low solvation energy are the most promising ones.

For each family of molecules, the adsorption and the solvation energies of the different studied molecules are plotted in Figure 5. On Figure 5a, it can be observed that for the protic molecules the adsorption energies are quite similar and barely correlated with the solvation energies. For aprotic molecules (Figure 5b), the adsorption energies fall in a small range (-1.00 to -1.44 eV) and there is no correlation at all between the adsorption and solvation energy. On the contrary, for the aromatic and the phosphorous compounds the adsorption energies are strongly correlated with the solvation energies as all the points are on the same line in Figure 5c and Figure 5d respectively. In other words, for these two families, a stronger surface adsorption is only achieved for molecules that also get better solvated.

Despite the qualitative character of our computations (neglect of entropy, single molecule adsorption, solid-vacuum interface and absence of the hydrophobic tail) compared to lubrication, our analysis suggests, for instance, that secondary amines are more likely to improve the wettability of a lubricant than phenols: the adsorption energy is very similar, but the solvation energy of the phenols is slightly higher. Similarly, given their similar adsorption energy, dimethylphosphate seems preferable over methylphosphate, since the latter is more solvated than the former. Moreover, by comparing their adsorption energies, our results suggest that carboxylic acids will improve the wettability properties better than alcohols or the esters, which is in agreement with a previous experimental study by Zhou et al. [START_REF] Zhou | Dynamic wetting of rolling oil on aluminum surfaces[END_REF] concluding that acids are improving the dynamic wetting of an oil on an aluminium surface better than alcohols and esters.

Conclusion

In this study, we have performed DFT computations to compare the solvation and adsorption energies of various head groups of lubricant additives in order to investigate which ones are the most suitable to favor the wetting of an alumina surface. We show that the adsorption energies of the studied molecules cover a wide range going from -0.64 to -2.41 eV. On the contrary, their solvation energies are quite similar and are much smaller in absolute value than the adsorption energies.

Our results show that acidic protic molecules, such as carboxylic acid, tend to be dissociated when adsorbed on the (100) surface of g-alumina, whereas dissociations do not occur for non-acidic protic head groups. The behavior of aromatic molecules generally follows from their functional groups, as the aromatic part does not stick on the surface. The studied phosphorous compounds exhibit a wide range of adsorption energies. The phosphates are the most adsorbed head groups, which suggests that the lubricant additives containing phosphate head groups should be better lubricant additives to improve

Figure 1 -

 1 Figure 1 -(a) Top view of the first atomic layer of the p(1x1) unit cell of the (100) surface and the corresponding labeling of the exposed adsorption sites. The other atomic layers are shaded. (b) Side view of the p(1x1) (100) g-alumina slab.

Figure 2 -

 2 Figure 2 -Adsorption of some selected protic molecules on the (100) surface of g-alumina: (a) water, (b) ethoxylated amine, (c) ethoxylated amine with the nitrogen atom adsorbed and (d) carboxylic acid. Distances are given in Å.

Figure 3 -

 3 Figure 3 -Correlation of adsorption energy with pKa for molecules exhibiting a dissociative adsorption.

Figure 4 -

 4 Figure 4 -Adsorption of some selected phosphorous compounds on the (100) surface of g-alumina: (a) methylphosphate, (b) O,O,O-trimethylthiophosphate, (c) O,O,S-trimethyldithiophosphate and (d) trimethylphosphite. Distances are given in Å.

Figure 5 -

 5 Figure 5 -Adsorption energy as a function of solvation energy for the studied additives: (a) protic molecules, (b) aprotic molecules, (c) aromatic compounds and (d) phosphorous compounds. * indicates dissociative adsorptions.

Table 1 -

 1

	Name	Structure	Name	Structure
	Water	H	O	H	Methanol	OH
	Polyalkylene glycol
		O			Methanethiol	SH
	(PAG) monomer			OH

Chemical structures of studied protic molecules. Tertiary amine is here classified as a protic molecule due to its similarity with the protic amines.

Table 2 -

 2 Chemical structures of studied aprotic molecules

	Name	Structure	Name	Structure
	Ether		Ketone
	(dimethylether)		(propanone)
	Ester (methyl
			Dimethylsulfide
	ethanoate)
	N,N-
			1-methyltriazole
	dimethylacetamide
	1-methyl-1,2,4-		1,2-dimethyl-4,5-
	triazole		dihydroimidazole

Table 3 -

 3 Chemical structures of studied aromatic molecules

	Name	Structure	Name	Structure
	Benzene		Phenol

Table 4 -

 4 Chemical structures of studied phosphorus-containing molecules

	Name	Structure	Name	Structure
	Methylphosphate		Dimethylphosphate
			O,O,O-
	Trimethylphosphate
			trimethylthiophosphate
	O,O,S-
			Dimethylphosphite
	trimethyldithiophosphate
	Trimethylphosphite
	3 Results and discussion
	3.1 Adsorption

Adsorption of molecules on the (100) surface of g-alumina is driven by acid/base characteristics of the molecule and surface site, as the aluminium atoms of g-alumina are hard Lewis-acidic adsorption sites.

Acknowledgements

S. Blanck gratefully acknowledges Total MS for her PhD fellowship. The authors would like to express their gratitude to the PSMN data center for CPU hours and assistance.

the wetting of the aluminium surface by the lubricant than the additives containing the other studied head groups.

This study shows that the behavior of a lubricant additive at the base oil/surface interface can be unraveled by DFT. While the solvation plays only a minor role and is not a discriminating parameter, the adsorption energy appears as a key parameter. We also identified that acidic additives adsorb more strongly on alumina, with a strength that can be correlated with a typical chemical property, the pKa.

The stronger the adsorption, the more likely is the formation of a film, which will, eventually, reduce friction and wear at the macroscopic scale. This DFT study can easily be extended to other surfaces of interest and other families of additive to better rationalize the choice of lubricant additives.

Supplementary material

Supplementary material to this article can be found online at: xxx All optimized structures of adsorbed molecules are available in the xyz format.