Jack Jon Hinsch

Assil Bouzid

Jordan Clive Barker

Jessica Jein White

Fabien Mortier

Huijun Zhao

Yun Wang
email: yun.wang@griffith.edu.au

Revisiting the Electrified Pt(111)/Water Interfaces through an Affordable Double-Reference Ab Initio Approach

Keywords: Grand canonical density functional theory, electrified platinum-water interface, potential of zero charge, capacitance, hydrogen adsorption

HAL is

Introduction

Electrified solid-liquid interfaces are of great importance in many energy and electrochemistry applications, including electrocatalysis, batteries, supercapacitors, solar cells, and fuel cells 1- 3 . The interface often determines the charge transfer rate of the electrochemical redox reactions, which is overseen by the applied bias electrode potential [START_REF] Che | Elucidating the Roles of Electric Fields in Catalysis: A Perspective[END_REF] . However, understanding the atomistic properties of the electrified solid-liquid interface and the thermodynamic state of interfacial processes is still challenging [START_REF] Doblhoff-Dier | Electric double layer of Pt(111): Known unknowns and unknown knowns[END_REF][START_REF] Steinmann | Understanding electrified interfaces[END_REF] . The difficulty arises from the operando characterizations on the solid-liquid interfaces [START_REF] Zaera | Probing Liquid/Solid Interfaces at the Molecular Level[END_REF] . To address these issues, atomistic modelling is a powerful and complementary tool to access the structural and electronic properties of the electrified interface [START_REF] Wang | Numerical Simulation of Electrified Solid-Liquid Interfaces[END_REF] .

Several ab-initio simulation strategies have been developed thus far to characterize the electrified solid-liquid interface [START_REF] Wang | Numerical Simulation of Electrified Solid-Liquid Interfaces[END_REF] . The computational hydrogen electrode (CHE) model is one of the widely used approaches for studying the electrified solid-liquid interface [START_REF] Bligaard | The Brønsted-Evans-Polanyi relation and the volcano curve in heterogeneous catalysis[END_REF][START_REF] Norskov | Towards the computational design of solid catalysts[END_REF] . In this method, the standard hydrogen electrode (SHE) serves as the reference potential. A pair of protons (H + ) and electrons (e -) are introduced in each elementary step to charge transfer processes [START_REF] Zhu | Coexisting Single-Atomic Fe and Ni Sites on Hierarchically Ordered Porous Carbon as a Highly Efficient ORR Electrocatalyst[END_REF][START_REF] Xiong | Cobalt Covalent Doping in MoS2 to Induce Bifunctionality of Overall Water Splitting[END_REF][START_REF] Wang | Electrolyte Effect on Electrocatalytic Hydrogen Evolution Performance of One-Dimensional Cobalt-Dithiolene Metal-Organic Frameworks: A Theoretical Perspective[END_REF] . However, this method dismisses the influence of the applied bias potential on the properties of the solid-liquid interface [START_REF] Liu | Theory-Experiment Gap[END_REF] . To accurately reveal the dynamic nature of the electrified solid-water interface, the ab-initio molecular dynamics (AIMD) method has been employed by several groups to investigate the Pt(111)-water interface, which is the most studied electrified solid-liquid interface. Otani and colleagues introduced an effective screening medium (ESM) method to explore half-reactions at electrode-electrolyte interfaces [START_REF] Otani | First-principles calculations of charged surfaces and interfaces: A planewave nonrepeated slab approach[END_REF] . Their findings indicate that the applied bias potential can significantly affect the properties of water layers near the surface. The water hydrogen-bond network is found to be incredibly disruptive at negative potentials. These potential changes force the water molecules to reorientate to accommodate the new electrostatic forces. Cheng et al. developed a computational standard hydrogen electrode (cSHE) method to understand the mechanisms of interface water molecules used to control the potential of zero charge (PZC) [START_REF] Li | Establishment of the Potential of Zero Charge of Metals in Aqueous Solutions: Different Faces of Water Revealed by Ab Initio Molecular Dynamics Simulations[END_REF][START_REF] Le | Molecular origin of negative component of Helmholtz capacitance at electrified Pt(111)/water interface[END_REF] . The PZC is the potential at which the electrode surface is electrically neutral, and it provides valuable information about the electrochemical properties of the system. Their study revealed that the double layer capacitance shape is a result of the interaction between the water and the surface.

The orientation of the adsorbed water molecules can greatly influence the interfacial properties.

Bouzid and Pasquarelllo developed a constant Fermi-level AIMD method to study the Pt(111)water interface with the consideration of variable electrode potential referenced to the SHE.

They observed significant double-layer reorganization dependent on the applied bias potential and accurately estimated the double-layer capacitance as well as the PZC of the metal electrode 18, 19 . Furthermore, Groß et al. employed AIMD with the RPBE-D3 functional to investigate the Pt(111)/water interface 20, 21 22 . Their study demonstrated the importance of accounting for van der Waals corrections since water molecules exhibited different symmetry and density with and without van der Waals corrections. While AIMD-based simulations have been successful in studying the Pt(111)/water interface, the computational expense of AIMD limits its applicability for large-scale investigations of different solid-liquid interfaces, highlighting the need for the development of computationally efficient yet accurate methods.

To this end, Taylor and colleagues developed a static density functional theory (DFT)-based technique that enables the calculation and the manipulation of the electrostatic potential in a two-phase system, providing a cost-effective approach for studying electrified solid interfaces [START_REF] Taylor | First principles reaction modeling of the electrochemical interface: Consideration and calculation of a tunable surface potential from atomic and electronic structure[END_REF] . To calculate the applied electrode potential relative to the Fermi level (EFermi), a thin vacuum layer is incorporated within the solution region as an internal reference potential, allowing for comparison to the standard hydrogen electrode (SHE) with an absolute reduction potential of 4.44 V at 298.15 K 24 . However, this artificial inclusion of a vacuum layer does not account for the nature of the vacuum/solution interface. To address the limitations of including an artificial vacuum layer in the simulation, a model that does not require the presence of a vacuum layer is sought, allowing for a more accurate representation of the solid-liquid interface.

A recent study by Haruyama et al. proposed an implicit solvation method to determine the electrode potential in a system comprising an electrode and an electrolyte solution [START_REF] Liu | Theory-Experiment Gap[END_REF] . In their model, the solvent is treated as a continuous medium, taking into account its collective effects on the electrode-electrolyte interface [START_REF] Mathew | Implicit solvation model for density-functional study of nanocrystal surfaces and reaction pathways[END_REF][START_REF] Mathew | Implicit self-consistent electrolyte model in plane-wave density-functional theory[END_REF] . By immersing the solute in a solvent "bath" and incorporating the averaging over the solvent's degrees of freedom, the simulation approach enables the study of the electrochemical double layer and allows for the manipulation of the system's potential [START_REF] Ringe | Implicit Solvation Methods for Catalysis at Electrified Interfaces[END_REF][START_REF] Cramer | Implicit Solvation Models: Equilibria, Structure, Spectra, and Dynamics[END_REF][START_REF] Kovalenko | Self-consistent description of a metal-water interface by the Kohn-Sham density functional theory and the three-dimensional reference interaction site model[END_REF][START_REF] Tomasi | Quantum Mechanical Continuum Solvation Models[END_REF] . However, the results obtained from ab-initio molecular dynamics (AIMD) simulations highlight the significance of including explicit water layers to accurately describe the interfacial properties of the system [START_REF] Otani | First-principles calculations of charged surfaces and interfaces: A planewave nonrepeated slab approach[END_REF][START_REF] Li | Establishment of the Potential of Zero Charge of Metals in Aqueous Solutions: Different Faces of Water Revealed by Ab Initio Molecular Dynamics Simulations[END_REF][START_REF] Le | Molecular origin of negative component of Helmholtz capacitance at electrified Pt(111)/water interface[END_REF] . Therefore, the simulation using an implicit solvation model alone may not be sufficient, as evidenced by the calculated potential of zero charge (PZC) of Pt(111) being significantly higher than the experimental value [START_REF] Hammer | Improved adsorption energetics within densityfunctional theory using revised Perdew-Burke-Ernzerhof functionals[END_REF] .

Recently, Xu et al. developed an affordable grand-canonical density functional theory (GC-DFT) approach with acceptable accuracy, which can provide an approximation of reality of the widely studied Pt(111)-water interface [START_REF] Xu | Optical method for quantifying the potential of zero charge at the platinum-water electrochemical interface[END_REF][START_REF] Le | Modeling Electrified Pt(111)-Had/Water Interfaces from Ab Initio Molecular Dynamics[END_REF] . To accurately capture the solvent environment, they employed a hybrid approach that includes both explicit and implicit solvation. This hybrid approach allowed us to capture the essential features of the electrified interface while maintaining computational efficiency. Despite its simplicity, the proposed model serves as a valuable starting point for understanding solid-liquid interfaces. This method has been successfully validated by comparing the calculated PZC and capacitance with experimental data [START_REF] Xu | Optical method for quantifying the potential of zero charge at the platinum-water electrochemical interface[END_REF][START_REF] Le | Modeling Electrified Pt(111)-Had/Water Interfaces from Ab Initio Molecular Dynamics[END_REF] , demonstrating its accuracy and reliability in capturing the electrochemical properties of the electrified interface. However, only small potential window was considered in this study.

As a result, the linear relationship between the applied potential and the surface charge was proposed, which was different from the previous AIMD results with a larger potential window [START_REF] Bouzid | Atomic-Scale Simulation of Electrochemical Processes at Electrode/Water Interfaces under Referenced Bias Potential[END_REF] . Moreover, the adsorption of hydrogen atoms on the electrified Pt(111) surface is known as the Volmer step for HER, which largely determines the energy conversion efficiency of this important electrocatalytic reaction for hydrogen production [START_REF] Liu | Theoretical Understanding of Electrocatalytic Hydrogen Production Performance by Low-Dimensional Metal-Organic Frameworks on the Basis of Resonant Charge-Transfer Mechanisms[END_REF][START_REF] Dong | First-Principles Determination of Active Sites of Ni Metal-Based Electrocatalysts for Hydrogen Evolution Reaction[END_REF] . Understanding the changes in electrochemical properties and the effects of H adsorption on the Pt(111) surface is valuable for elucidating the mechanisms of electrochemical reactions and developing improved catalysts for hydrogen-related applications. However, the impact of adsorption hydrogen atom on the interfacial properties are unclear.

In this study, we used the GC-DFT method with the hybrid solvation model to reinvestigate the electrified Pt(111)/water interface. Our results reproduced the relationship between the applied bias potential and the surface charge found in the AIMD simulations. The subtle impact of the adsorbed hydrogen atom on the properties of the electrified interface was also revealed.

Computational Details

In this study, the Vienna Ab-initio Simulation Package (VASP) was utilized for all computations, employing the projected augmented wave (PAW) method 36 37 . The electron-ion interaction was described using the PAW pseudopotentials, where the valence electronic configurations of Pt, H and O are 6s 1 5d 9 , 1s 1 , and 2s 2 2p 4 , respectively. The kinetic energy cutoff value was set to 600 eV to expand the smooth part of the wave function, as suggested in previous studies [START_REF] Mathew | Implicit self-consistent electrolyte model in plane-wave density-functional theory[END_REF] . During structural optimizations, all atoms relax until the Hellmann-Feynman forces were smaller than 5 × 10 -2 eV/Å. The convergence criterion for the electronic self-consistent loop was set to 1 × 10 -5 eV. The gamma-centered k-point meshes were generated with a reciprocal space resolution of 2π × 0.04 Å -1 . The optPBE exchangecorrelation (XC) functional was chosen for this study as it showed success at predicting solids and water properties [START_REF] Yong | optPBE-vdW density functional theory study of liquid water and pressure-induced structural evolution in ice Ih[END_REF][START_REF] Avelar | van der Waals exchange-correlation functionals over bulk and surface properties of transition metals[END_REF] . The optPBE functional incorporates nonlocal correlation contributions into the XC energy, allowing for the consideration of less spherical densities and providing a better approximation of weak interactions. [START_REF] Yong | optPBE-vdW density functional theory study of liquid water and pressure-induced structural evolution in ice Ih[END_REF][START_REF] Klimeš | Chemical accuracy for the van der Waals density functional[END_REF][START_REF] Malone | A van der Waals Inclusive Density Functional Theory Study of the Nature of Bonding for Thiophene Adsorption on Ni(100) and Cu(100) Surfaces[END_REF] .

The Pt(111) surface was presented by a four-layer (2 × √3 ) Pt(111) slab. Based on previous studies, it has been shown that a single explicit water layer is sufficient to capture the characteristics of the inner Helmholtz layer [START_REF] Gauthier | Solvation Effects for Oxygen Evolution Reaction Catalysis on IrO2(110)[END_REF] . Therefore, one layer of explicit water bilayer was placed on both sides of the Pt slab here, while an implicit water layer with a thickness of 20 Å was also included. Each explicit water bilayer was composed of eight water molecules arranged in a honeycomb configuration. The implicit continuum solvation model was implemented using the VASPsol package, with a permittivity value of 78.40 and a Debye length of 3.040 Å [START_REF] Mathew | Implicit solvation model for density-functional study of nanocrystal surfaces and reaction pathways[END_REF][START_REF] Mathew | Implicit self-consistent electrolyte model in plane-wave density-functional theory[END_REF] .

To simulate an applied bias potential (U), the double layer charge of the system was varied [START_REF] Lozovoi | Ab initio simulation of charged slabs at constant chemical potential[END_REF] .

The structures were then further fully optimized. The applied bias potential was derived from the Fermi energy of the interface systems (EFermi) with respect to the SHE potential using a double-reference method, as illustrated in Fig. 1. In this approach, the Fermi level is first referenced to the electrostatic potential in the center of the electrolyte, which is systematically shifted to zero (Eshift), as required from the Poisson-Boltzmann equation. Then, the electrostatic potential of the implicit solvent is further referenced to the SHE potential using an additional empirical alignment term of 4.6 V, which was obtained as the difference between the experimental and computed PZC values using the implicit solvation method [START_REF] Mathew | Implicit self-consistent electrolyte model in plane-wave density-functional theory[END_REF] . Hence, the applied bias potential in terms of the SHE (U(V) vs. SHE) can be calculated as below: To better understand the interaction between water and the surface, we calculated the dipole orientation between the bisector of the hydrogen bonds of water and the surface normal [START_REF] Le | Modeling Electrified Pt(111)-Had/Water Interfaces from Ab Initio Molecular Dynamics[END_REF] . Therefore, the average distance of the water layer to the Pt surface as well as the bond angle between OH bond and the normal of the surface were calculated. The average distance between the Pt surface and the water layer was determined by calculating the average distance between the z positions of the oxygen atoms in the water molecules and the outmost Pt layer. The bond angle was calculated as follows,

𝑈(𝑉) 𝑣𝑠. 𝑆𝐻𝐸 = -(𝐸 𝐹𝑒𝑟𝑚𝑖 + 𝐸 𝑠ℎ𝑖𝑓𝑡 )/|𝑒| -4.6

𝜃 = 𝑎𝑟𝑐𝑐𝑜𝑠𝑖𝑛𝑒 ( 𝑍(𝐻) -𝑍(𝑂) 𝑑 𝑂-𝐻 )

where 𝑍(𝐻) -𝑍(𝑂) is the difference between z coordinates of H and O atoms with the vector from O towards H. And dO-H is the corresponding O-H bond length in water. This value can then be used to calculate the dipole orientation of the water molecule.

Three sites were considered for H adsorption (see Fig. S1). It was found that the atop site is the most energetically preferred one, which agree with previous conclusions [START_REF] Le | Modeling Electrified Pt(111)-Had/Water Interfaces from Ab Initio Molecular Dynamics[END_REF][START_REF] Le | Modeling electrified metal/water interfaces from ab initio molecular dynamics: Structure and Helmholtz capacitance[END_REF] . Consequently, only atop site will be considered here. The adsorption energy of hydrogen atom (Eads) was calculated as below:

𝐸 𝑎𝑑𝑠 (𝑉) = 𝐸 𝑇𝑜𝑡𝑎𝑙 (𝑉) -𝐸 𝐶𝑙𝑒𝑎𝑛 (𝑉) -𝐸 𝐻
Here, 𝐸 𝑇𝑜𝑡𝑎𝑙 is the total energy of the interface with adsorbed H atom at a specific electrode potential, 𝐸 𝐶𝑙𝑒𝑎𝑛 (𝑉) is the total energy of the pristine interface at a specific electrode potential, and 𝐸 𝐻 is the half of the energy of a H2 molecule in an isolated cell with the continuum method.

The hydrogen binding free energy was also calculated to better understand the binding properties. The ∆𝐺 𝐻 * was calculated as below:

∆𝐺 𝐻 * = ∆𝐸 𝑎𝑑𝑠 + ∆𝑍𝑃𝐸 -𝑇∆𝑆
Here, ZPE is the zero point energy [START_REF] Nørskov | Origin of the Overpotential for Oxygen Reduction at a Fuel-Cell Cathode[END_REF] , T is temperature according to the recent experiment and ∆𝑆 entropic contributions retrieved from the literature [START_REF] Nørskov | Origin of the Overpotential for Oxygen Reduction at a Fuel-Cell Cathode[END_REF][START_REF] Yang | The Nature of Hydrogen Adsorption on Platinum in the Aqueous Phase[END_REF] .

Results and Discussion

Atomic interface model

From Fig. 1, the double-reference model was employed to compute the properties of the Pt(111)-water interface. The interface was surrounded by a dielectric continuum representing the average electrostatic attraction of the water molecules. The adequacy of this solvation model has been a subject of prolonged debate within the scientific community. Bramley et al.

and Akinola believed that the DFT-based implicit solvent is problematic [START_REF] Akinola | Adsorption Energies of Oxygenated Aromatics and Organics on Rhodium and Platinum in Aqueous Phase[END_REF][START_REF] Bramley | Understanding Adsorption of Organics on Pt(111) in the Aqueous Phase: Insights from DFT Based Implicit Solvent and Statistical Thermodynamics Models[END_REF][START_REF] Bramley | Reconciling Work Functions and Adsorption Enthalpies for Implicit Solvent Models: A Pt (111)/Water Interface Case Study[END_REF] . However, several groups also suggested that the solvation model can effectivley predict interfacial properties 14- 17, 50 . Fig. 2a shows the calculated electrostatic potential of the Pt(111) in the absence of an explicit water layer but in the presence of the implicit solvation model. The calculated PZC value obtained is 1.09 V vs SHE, which significantly differs from the experimentally reported value of 0.4 V. [START_REF] Climent | Coulostatic Potential Transients Induced by Laser Heating of a Pt(111) Single-Crystal Electrode in Aqueous Acid Solutions. Rate of Hydrogen Adsorption and Potential of Maximum Entropy[END_REF][START_REF] Langenbach | The adsorption of water on Pt(111) studied by irreflection and UV-photoemission spectroscopy[END_REF][START_REF] Pajkossy | Double layer capacitance of Pt(111) single crystal electrodes[END_REF][START_REF] Trasatti | Structure of the metal/electrolyte solution interface: new data for theory[END_REF][START_REF] Tripkovic | Standard hydrogen electrode and potential of zero charge in density functional calculations[END_REF] . Our calculated PZC value from the implicit model is close to a similar work by Mathew et al. [START_REF] Hammer | Improved adsorption energetics within densityfunctional theory using revised Perdew-Burke-Ernzerhof functionals[END_REF] , which is 0.85 V. The small difference may be ascribed to the use of different XC functionals. Our findings support the conclusion made by Mathew et al. that the implicit solvation model alone is inadequate for accurately reproducing experimental observations [START_REF] Hammer | Improved adsorption energetics within densityfunctional theory using revised Perdew-Burke-Ernzerhof functionals[END_REF] .

Subsequently, an explicit water layer was added to the top side of the Pt(111) surface, resulting in an asymmetric system. However, the electrostatic potential of the system exhibits poor convergence in the bulk solvent region, as indicated by the inset of Fig. 2b. To address the convergence issues, one possible approach is to use a very thick implicit solvation layer, although this would significantly increase computational cost. Another method is the dipole correction, which helps counterbalance the slow convergence with respect to the size of the supercell [START_REF] Makov | Periodic boundary conditions in ab initio calculations[END_REF] . However, the dipole correction method is not suitable for charged systems, and currently, it can only be applied to cubic supercells in the VASP code implementation.

In order to address the asymmetry of the electrostatic potential, we employed a "sandwich" model by mirroring the system along a perpendicular plane. This symmetric configuration enables rapid convergence of the electrostatic potential in the implicit water region, shown in

Water orientation

Various initial configurations were assessed, each featuring distinct water orientations. These configurations were labeled based on the orientation of OH bonds pointing away or towards the surface. These models included configurations where the hydrogen in water pointed away from the surface (H-up), towards the surface (H-down), and a mixed orientation (H-mix) configuration (see Fig. 3). The H-down orientation exhibited the best performance among the studied models, as it showed the closest agreement with experimental results in terms of PZC values. The experimental potential of zero charge (PZC) is around 0.4 V vs SHE [START_REF] Climent | Coulostatic Potential Transients Induced by Laser Heating of a Pt(111) Single-Crystal Electrode in Aqueous Acid Solutions. Rate of Hydrogen Adsorption and Potential of Maximum Entropy[END_REF][START_REF] Langenbach | The adsorption of water on Pt(111) studied by irreflection and UV-photoemission spectroscopy[END_REF][START_REF] Pajkossy | Double layer capacitance of Pt(111) single crystal electrodes[END_REF][START_REF] Trasatti | Structure of the metal/electrolyte solution interface: new data for theory[END_REF][START_REF] Tripkovic | Standard hydrogen electrode and potential of zero charge in density functional calculations[END_REF] . Our PZC with the H-down configuration is 0.44 V vs SHE. Additionally, the water bilayer with the Hmixed and H-up configurations have a much lower PZC values of -0.01 and -0.50 V, respectively. Hence, even minor modifications to the explicit water layer in this interface configuration can yield substantial alterations in the overall interface properties, consistent with the findings of Xu et al. [START_REF] Xu | Optical method for quantifying the potential of zero charge at the platinum-water electrochemical interface[END_REF] . However, there are slight deviations between our calculated PZC values and the values reported by Xu et al., which were measured at 0.23 V vs. SHE [START_REF] Xu | Optical method for quantifying the potential of zero charge at the platinum-water electrochemical interface[END_REF] . This disparity can be attributed to various factors. In their research, the D3 vdW correction method and a low energy cut-off of 400 eV were employed. Additionally, their study incorporated a difference of 4.44 V between the electrostatic potential of the implicit solvent and the SHE potential.

The H-down orientation the most energetically stable one, which total energy is 0.04 and 0.10 eV lower than that of H-Mix and H-Up (see Table. S1). The H-down water orientation configuration was also demonstrated as the most energetically favorable and most likely to be observed at the interface by several groups from previous simulations [START_REF] Le | Molecular origin of negative component of Helmholtz capacitance at electrified Pt(111)/water interface[END_REF][START_REF] Xu | Optical method for quantifying the potential of zero charge at the platinum-water electrochemical interface[END_REF][START_REF] Le | Modeling Electrified Pt(111)-Had/Water Interfaces from Ab Initio Molecular Dynamics[END_REF][START_REF] Bouzid | Atomic-Scale Simulation of Electrochemical Processes at Electrode/Water Interfaces under Referenced Bias Potential[END_REF][START_REF] Le | Modeling electrified metal/water interfaces from ab initio molecular dynamics: Structure and Helmholtz capacitance[END_REF] . It further supports that the H-down configuration is the dominant one in the Pt(111)-water interface.

Building upon these findings, we proceeded to conduct calculations on the H-down configuration for the remainder of the study.

The calculations demonstrate the essential role of water molecules in the model for accurately predicting the electrochemical properties of the interface. The inclusion of water molecules allows for a more realistic representation of the system, leading to improved agreement with experimental observations. Without considering the presence of water, the predictions may deviate significantly from the actual behavior of the electrified interface. Indeed, the orientation of water molecules and hydrogen bonding play a crucial role in electron transfer and intermolecular interactions at the electrified interface [START_REF] Nishino | Direct Measurement of Electron Transfer through a Hydrogen Bond between Single Molecules[END_REF] . Implicit solvent methods, which approximate the solvent as a continuous medium, often struggle to accurately capture these short-range interactions, including hydrogen bonding [START_REF] Zhang | Evaluation of Generalized Born Models for Large Scale Affinity Prediction of Cyclodextrin Host-Guest Complexes[END_REF][START_REF] Zhang | Comparison of Implicit and Explicit Solvent Models for the Calculation of Solvation Free Energy in Organic Solvents[END_REF] . The explicit inclusion of water molecules in the model allows for a more detailed description of the intermolecular interactions, enabling a more accurate prediction of electron transfer processes and other properties of the electrified interface. Based on our results, further investigations in this work will focus on the H-down configuration as it represents the most dominant configuration for the electrified interface. The impact of the thickness of the implicit solvation layer on the PZC values was investigated, considering implicit layers of 20, 25, and 30 Å thickness (see Fig. S2). The calculated PZC values were found to be nearly identical, indicating that a thickness of 20 Å for the implicit solvation layer is sufficient to obtain converged properties of the electrified Pt(111)-water interface.

Double Layer Capacitance

The computation of the double layer capacitance based on our models revealed two linear trends in the evolution of double layer charge vs. electrode potential, as shown in Fig. 4a. The capacitance values are obtained by taking the derivative of the surface charge with respect to the potential. In Fig. 4b, a finite difference approach was used to calculate the derivative between two adjacent points, resulting in an average capacitance within that range. A small potential range from -1 to 1 V vs. SHE was considered, which is wider than that used in previous studies [START_REF] Xu | Optical method for quantifying the potential of zero charge at the platinum-water electrochemical interface[END_REF][START_REF] Le | Modeling Electrified Pt(111)-Had/Water Interfaces from Ab Initio Molecular Dynamics[END_REF][START_REF] Le | Modeling electrified metal/water interfaces from ab initio molecular dynamics: Structure and Helmholtz capacitance[END_REF] . Small ranges are often used in both computational and experimental studies to prevent unwanted electrochemical reactions such as electrolysis [START_REF] Doblhoff-Dier | Electric double layer of Pt(111): Known unknowns and unknown knowns[END_REF] . Additionally, one water layer would likely be insufficient to describe the interaction at higher potentials [START_REF] Doblhoff-Dier | Electric double layer of Pt(111): Known unknowns and unknown knowns[END_REF] .

Our results suggest a change in the capacitance at around 0.5 V. In the lower potential region (< 0.5 V) around the PZC, the double layer capacitance was calculated to be 21.64 𝜇𝐹/cm 2 , which agrees well with experimental measurements of about 20 𝜇𝐹 /cm 2 53 , while higher potential range (> 0.5 V) presented a different capacitance of 46.59 𝜇𝐹/cm 2 . Fig. 4b illustrates the relationship between the applied potential and the capacitance, and it demonstrates a close reproduction to the experimental measurements [START_REF] Bouzid | Atomic-Scale Simulation of Electrochemical Processes at Electrode/Water Interfaces under Referenced Bias Potential[END_REF][START_REF] Ojha | Double Layer at the Pt(111)-Aqueous Electrolyte Interface: Potential of Zero Charge and Anomalous Gouy-Chapman Screening[END_REF] . This indicates that our computational models effectively capture the behavior of the capacitance as a function of the applied potential.

It is worth noting that the capacitance change shown in Fig. 4b is noisy, which was expected as this approach is oversimplified and neglecting the dynamic feature of the water.

Additionally, electrolytes are not explicitly described here. And electrolyte adsorption and the formation of the outer Helmholtz plane are essential components of the electrified solid liquid interface at these potentials [START_REF] Doblhoff-Dier | Electric double layer of Pt(111): Known unknowns and unknown knowns[END_REF][START_REF] Favaro | Unravelling the electrochemical double layer by direct probing of the solid/liquid interface[END_REF] , which may lead to the theory-experiment discrepancy observed in Fig. 4b.

Structural properties

The structural properties of the interface under an applied bias potential play a crucial role in the initial steps of electrochemical processes, including oxygen evolution reaction (OER), hydrogen evolution reaction (HER), and carbon dioxide reduction reaction (CO2RR). In this study, the water position to the surface provides information about the proximity of the water molecules to the Pt(111) surface. It shows how the distance between the water layer and the surface changes with the applied potential, reflecting the interaction between the water molecules and the electrified surface. The average dipole orientation of the water molecules with respect to the surface normal indicates the orientation and arrangement of the water molecules at the interface. It provides insights into the hydrogen bonding network and the structural properties of the water layer in the presence of the electrified Pt(111) surface.

Fig. 5 shows the change of the water height and dipole orientation as a function of the applied bias potential. Our findings indicate that the double layer charge of the interface exhibits two distinct linear domains, which are influenced by the distance between the metal and the water layer as well as the orientation of the water molecules. At positive potentials, the distance between the water and metal interface decreases with increasing potential, while it increases at negative potentials. Furthermore, the analysis of the dipole orientation indicates that water molecules point away from the interface under positive potentials. This behavior can be attributed to the electrostatic interaction between the positively charged surface and the negatively charged oxygen atoms, causing the hydrogens to be pushed aside. These observations are aligned with findings from other studies [START_REF] Le | Molecular origin of negative component of Helmholtz capacitance at electrified Pt(111)/water interface[END_REF][START_REF] Le | Modeling Electrified Pt(111)-Had/Water Interfaces from Ab Initio Molecular Dynamics[END_REF] . Under positive potentials, an interesting phenomenon is observed where water molecules approaching the surface can experience collisions with each other. At a potential of 0.5 V, a majority of the water molecules are adsorbed on the surface in an H-up configuration. This leads to a disruption and decreased organization of the water layer. In contrast, at lower potentials, a more rigid ice-like layer forms as water molecules move away from the surface. The height variance of all water molecules was calculated to assess the degree of spatial organization (see Fig. S3). Variance is a statistical measure that quantifies the dispersion of a data set around the average. Using this variable, the organization of the water molecules can be quantified. At a potential of 0.2 V, the height variance of water was found to be 0.37, indicating a relatively lower level of spatial organization. However, at 0.7 V, the variance increased to 0.42, suggesting a higher degree of disorder and less spatial organization within the water layer. This change in variance further supports the observation that the water layer becomes disrupted and less organized at higher potentials. These observations highlight the dynamic behavior of the water layer and its dependence on the applied potential, resulting in changes in the overall structure and organization of the interface.

The orientation of water molecules at the electrified interface indeed plays a crucial role in determining the electrochemical properties. Specifically, at lower potentials H2O contribute to the formation of hydrogen bonding networks that enhance the stability and organization of the water layer, which leads to a denser and more organized water layer [START_REF] Li | In situ probing electrified interfacial water structures at atomically flat surfaces[END_REF] . The enhanced stability and organization of the water layer in these orientations can affect various electrochemical processes. Xu, et al. employed a comparable methodology and yielded outcomes consistent with our findings [START_REF] Xu | Optical method for quantifying the potential of zero charge at the platinum-water electrochemical interface[END_REF] . Their observations also highlighted the pronounced influence of water structures on interface behavior. This congruence underscores the reproducibility of our approach and its potential applicability in other research areas.

H adsorption

The influence of applied bias potentials on the adsorption properties of H atoms on the electrified Pt(111)-water interface were also investigated here. The atop site was the energetically preferred with the lowest energies (see Table . S2). The PZC of the system with the adsorbed H atom at the atop site is 0.38 V, which is lower than that without adsorbed H.

The lower PZC value suggests that the presence of H alters the double layer charge distribution and influences the electrochemical properties of the interface. This observation is consistent with the fact that H adsorption on Pt surfaces can affect the surface chemistry and catalytic activity, particularly in the context of the HER. The aqueous hydrogen binding free energy change (∆𝐺 𝐻 * ) was also computed at a calculated zero-point energy change (ΔZPE) of 0.24 eV, energy change (ΔE) as -0.26 eV, and a entropy change (ΔS) at 343K obtained from Yang et al. [START_REF] Yang | The Nature of Hydrogen Adsorption on Platinum in the Aqueous Phase[END_REF] . The resulting ∆𝐺 𝐻 * value was determined to be 0.18 eV, which exhibits an exceptionally close correspondence to the experimental value of 0.10 eV [START_REF] Yang | The Nature of Hydrogen Adsorption on Platinum in the Aqueous Phase[END_REF] . This convergence between the calculated and experimental values underscores the reliability of the method and its capability to predict the properties of the solid-liquid interface. Similar to the Pt(111)-water interface without adsorbed H, Fig. 7a shows that the double layer charge for the Pt(111) surface with adsorbed H exhibits two linear trends as the applied potentials increase. This suggests that there are two distinct regimes in the behavior of the double layer charge as a function of the applied potential. There are several distinctions caused by the hydrogen adsorption. The presence of adsorbed hydrogen on the Pt(111) surface leads to a reduction in capacitances at the low (< 0.5 V vs. SHE) and high applied bias potentials (> 0.5 V vs. SHE). The reduction in capacitance can be attributed to several factors. Firstly, the adsorption of hydrogen atoms on the surface may alter the electronic structure and charge distribution, leading to changes in the surface capacitance. Additionally, the presence of adsorbed hydrogen can influence the surface chemistry and the interaction of the electrode with the water, thereby affecting the charge storage capacity.

The adsorption energies of hydrogen on the Pt(111) surface are influenced by the applied potential. Similar to the double layer charge, the adsorption energy can be described by two linear relationships as the potential is varied (see Fig. 7b). At high applied potentials, the adsorption of hydrogen on the Pt(111) surface becomes significantly stronger. This can be attributed to the increase in double layer charge resulting from the higher negative potential.

When the potential is raised, the surface carries more positive charge (See Fig. 7a). The increased double layer charge enhances the electrostatic attraction between the negatively charged O atoms in water and the positively charged surface, resulting in a more stable adsorption configuration. [START_REF] Le | Modeling Electrified Pt(111)-Had/Water Interfaces from Ab Initio Molecular Dynamics[END_REF] . The height of water and its dipole orientation at the Pt(111)-water interface with the adsorbed H atoms closely resembled those observed on the interface without adsorbed H atoms (see Figs. 7c andd). The water height and dipole orientation of the water molecules varied depending on the applied potential with two distinct linear domains. At the point of zero charge, the presence of the protruding H atom at the Pt(111)-water interface caused a bending effect on the water layer. This occurred because neighboring water molecules repelled each other through hydrogen bonding interactions. As the applied potential increased, this effect became stronger, leading to a maximum height variance in the water position of 0.37 at 0.63 V.

Conversely, as the potential decreased, the water molecules exhibited a tendency to move away from the surface, resulting in the formation of a flatter water layer. The variance in the distance between the water molecules and the surface decreased with decreasing potential, indicating that the water molecules were approximately located at the same height. For example, at -1.67

V, the height variance was measured to be 0.19. The dipole orientation results are also close to previous observations without adsorbed H regarding the orientation of water molecules at the Pt(111)-water interface under different potentials (see Fig. 7d). Under positive potentials, the orientation of the water molecule's hydrogen atom is directed away from the surface, whereas under negative potentials, it points towards the surface. A distinct disparity emerged following hydrogen adsorption. The model with adsorbed hydrogen exhibited a broader dipole range compared to the clean surface, and the distance between water and the surface was greater. The presence of adsorbed hydrogen amplified the trends observed on the clean surface (see Fig. 4567). Therefore, the presence of the adsorbed hydrogen species influences the behavior of the water molecules. The orientation analysis reveals that the water molecules adjacent to the adsorbed H atoms adjust their positions to accommodate the adsorbed hydrogen into the hydrogen bonding network. The adsorbed H atom prevents the neighboring water molecules from approaching too close to the surface.

Interestingly, despite the observed changes in adsorption strength with applied potential, the bond length between the adsorbed H atom and surface Pt atom remained relatively constant at approximately 1.5 Å. Given that the trends in adsorption energy and water height show a striking similarity, it indicates that the adsorption energy is influenced not only by the Pt-H bond but also by the interactions between water and the adsorbed hydrogen. The strong adsorption energy at positive potentials can be attributed to the integration and support of the adsorbed hydrogen species within the water network. In contrast, at low potentials, the water molecules move away from the surface with H atoms in water pointing downwards, leaving no available positions for the adsorbate to interact with the water network. As a result, the adsorption energies change more slowly along the applied potential. These findings highlight the intricate interplay between the adsorbed hydrogen, water molecules, and the electrochemical environment, which affects the adsorption behavior and overall properties of the Pt(111)-water interface.

Nevertheless, it is essential for future investigations to expand the scope of this model by integrating other solvent parameters like electrolytes and pH. Furthermore, while our current study has primarily examined the interface under low potentials to avoid undesired chemical reactions, it's crucial to acknowledge that many heterogeneous electrocatalytic reactions necessitate higher voltages. In such instances, a single water layer might not suffice for accurately modeling these electrified interfaces. The consideration of water coverage and the establishment of inner and outer Helmholtz planes are likely imperative.

Conclusion

In conclusion, we used an affordable approach based on the GC-DFT method to reinvestigate the electrified Pt(111)-water model system. We used a hybrid approach with a "sandwich" design to combine explicit and implicit solvation methods, which can overcome convergence issues and provide more reliable results. The relationship between electrochemical properties including double layer charge and capacitance and the applied bias potential derived from our GC-DFT results shows correlation with two distinct linear domains, in good agreement with experimental measurements. The great match highlights the importance of considering explicit water molecules and their interactions in modeling electrified solid-liquid interfaces accurately.

Further analysis on the structural properties of the interface demonstrates that the orientation and position of water molecules play a significant role in determining the PZC and double layer capacitance of the system. At high potentials, the water layer becomes disordered and moves away from the surface, while at low potentials, the water molecules point towards the surface, stabilizing the structure. The adsorption of hydrogen on Pt also affects the system's performance, with different adsorption energies observed at different potentials.

This affordable GC-DFT approach provides a valuable and efficient means to enhance our understanding of electrified solid-liquid interfaces. By combining the advantages of implicit and explicit solvation methods, this approach offers an accurate description of the interactions and properties of these interfaces while minimizing computational costs. The agreement between the findings of this approach and those reported in AIMD simulations further validates its effectiveness. Overall, the affordable GC-DFT approach enables researchers to study electrified solid-liquid interfaces in a cost-effective manner without compromising accuracy, contributing to our understanding of these important systems.

Figure 1 .

 1 Figure 1. Illustration of the double reference method for calculating the applied bias potential

Fig. 2c .

 2c Fig. 2c. The use of dipole correction becomes not necessary in this case, allowing for the

Figure 2 .

 2 Figure 2. (a) The side view of the electrified Pt(111)-water interface using the implicit-only

Figure 3 .

 3 Figure 3. Optimized atomic structures of three Pt(111)-water interface with specific water

Figure 4 .

 4 Figure 4. (a) The double layer charge as a function of the calculated applied bias potential. The

Figure 5 .

 5 Figure 5. (a) Water-surface distance and (b) Average dipole orientation of the water molecules

Figure 6 .

 6 Figure 6. H adsorption on the electrified Pt solid-liquid interface. The number underneath the

Figure 7 .

 7 Figure 7. The dependence of (a) double layer charge, (b) H adsorption energy, (c) height of

Acknowledgments: This research was conducted on the supercomputers in National Computational Infrastructure (NCI) in Canberra, Australia, which is supported by the Australian Commonwealth Government, and Pawsey Supercomputing Centre in Perth with the funding from the Australian government and the Government of Western Australia. Calculations were performed by using resources from Grand Equipement National de Calcul Intensif (GENCI, grant no. A0120913426). Computational resources provided by the computing facilities Mésocentre de Calcul Intensif Aquitain (MCIA) of the Université de Bordeaux and of the Université de Pau et des Pays de l'Adour. We acknowledge funding from CNRS Institute of Chemistry through the "International Emerging Actions 2022" mobility grant (2DH2 project) and from Agence national de recherche under the RECIFE ANR-DFG project (Grant Number ANR-21-CE08-0036-01).

Supplementary Information

Additional supporting information can be found attached to accommodate the manuscript.


Evgenii M Roginskii
email: e.roginskii@mail.ioffe.ru

Raghvender Raghvender

Olivier Noguera

Philippe Thomas

Olivier Masson

Assil Bouzid
email: assil.bouzid@unilim.fr

Ab Initio Study of the Raman Spectra of Amorphous Oxides: Insights into the Boson Peak Nature in Glassy TeO 2

Keywords: Raman spectra, DFT, Boson peak, TeO2, glass

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

Despite this progress in achieving a good understanding of the structure of glassy TeO 2 , research on the vibrational and non-linear optical properties of this glass remains scarce. At a more general level, the compu- tational study of the spectroscopic properties of amorphous systems is limited to the study on atomic clusters that are extracted from the amorphous or crystalline models. While instructive, this procedure is not straightforward and provides an inherently limited understanding of the structural correlations behind the Raman response of the glassy material. From a technical standpoint, the main stumbling block in computing spectroscopic and optical properties is the huge computational effort that is required when dealing with amorphous systems within a periodic DFT calculation. On the one hand, a proper description of these glasses requires an atomistic model that contains several hundreds of atoms. On the other hand, the present methodologies are either suitable for small size systems or only work well in the case of crystalline materials [18,19,20,21,22,23]. Consequently, there is a need for a method that is able to compute accurate Raman properties on large periodic amorphous models at a reasonable computational cost. This method would enable us to access the important quantities that are required to properly correlate the macroscopic properties of the material to its microscopic structure, thereby setting the scene to better design the roots of functional oxide materials.

In this letter, we present a scheme that can compute the Raman spectra of amorphous systems that contain several hundreds of atoms in a periodic simulation cell, which is based on the linear response perturbation theory and the finite differences method. This scheme yields a Raman spectra that is in excellent agreement with the experimental measurements. Furthermore, our scheme allows di- rect access to the atomistic nature of the normal modes, which facilitates the interpretation of the observed Raman bands, especially in the case of the Boson peak. In this letter, we focus on the Raman properties of amorphous TeO 2 by taking advantage of the achieved atomic scale models in Ref. [17], and we then compare the results to those on amorphous SiO 2 . Computational details related to the TeO 2 amorphous model generation protocol can be found in Ref. [17]. The details of the SiO 2 model are provided in the supplementary materials (SI).

At the microscopic level, Raman scattering phenomenon is related to photon-phonon interaction. The result of this inelastic interaction is the change of the energy according to the energy conservation law. In particular, when the photon energy decays after the scattering hω s = hω i -hΩ (where ω s and ω i are the frequencies of the scattered and incipient light, respectively, and Ω is the frequency of the phonon), then this process is called Stokes. When the photon energy increases after scattering, the process is then called Antistokes. The intensity of the Raman scattering depends on the change of the polarizability with respect to atomic displacements which correspond to a selected vibrational mode.

In practice, we start by focusing on the dynamical properties that are computed within the Born-Oppenheimer approximation, which itself is based on the double-harmonic approximation. Given an amorphous configuration, we first optimize the geometry (i.e., variation of the lattice parameters and atomic positions at fixed angles) to achieve a ground state at the minimum of the total energy. Phonon calculations are then performed to extract the vibration frequencies and the normal vectors using the PHONOPY software package [24] to compute the dynamic matrix using the finite difference method.

In a Raman spectrum, the positions of the bands are determined by the vibrational mode frequencies near the Γpoint and in the nonresonant Raman scattering processes. The intensities of the bands can be calculated within a Placzek approximation [START_REF] Born | Dynamical Theory of Crystal lattices[END_REF]. For the Stokes process, the intensity of the m th vibrational mode is given by [START_REF]Light scattering in solids II[END_REF]:

I m = (ω 0 -ω m ) 4 ω m |e i • A m • e s | 2 (n m + 1), (1)
where e i and e s are the unit polarization vectors of the incident and scattered radiation, respectively; ω 0 is the frequency of incipient radiation; ω m is the frequency of the phonon mode; n m is the Bose-Einstein distribution given by n(ω m ) = 1/(e hωm/k B T -1); and A m is the Raman tensor, which is defined by:

A m ij = √ Ω � α,β ∂χ ij ∂r α,β U m α,β , (2)
where Ω is the unitcell volume, r α,β is α coordinate of β-th atom in the unitcell, and U m α,β corresponds to normal vector of the m-th vibrational mode. The linear dielectric susceptibility tensor χ ij is obtained within density functional perturbation theory [START_REF] Souza | [END_REF] by obtaining the derivative of the macroscopic polarization P with respect to an electric field E:

χ ij = ∂P i ∂E j = ∂ 2 E tot ∂E i ∂E j (3)
Therefore, the Raman tensor A m ij is a third order derivative of the total energy with respect to the electric field twice and atomic displacements. While one can use the (2n+1) theorem in quantum mechanics to calculate the third order derivative of the total energy, this technique requires considerable computational effort and implementation of theoretical formalism because one has to access the perturbed wave functions up to the first order. Alternatively, the finite differences method can be applied to calculate the Raman tensor as defined in Eq. ( 2). This method is based on the generation of structures with atomic displacements along each normal mode in positive and negative Copyright line will be provided by the publisher directions. This is followed by calculations of susceptibility tensor χ ij for each distorted structure using the linear response approach as implemented in CP2K code [START_REF] Souza | [END_REF]. Finally, the values of the Raman tensor components for each vibrational mode are obtained as derivatives of susceptibility with respect to normal vectors obtained by the finite difference methods. Technical details of the implementa- tion of the finite differences method are provided in the supplementary materials.

We note that the following rotational invariants [START_REF] Long | The Raman effect: A unified treatment of the theory of Raman scattering by molecules[END_REF] can be used in the case of an amorphous sample:

ζ i = 1 3 � α i xx + α i yy + α i zz � , (4)
γ 2 i = 1 2 � (α i xx -α i yy ) 2 + (α i yy -α i zz ) 2 + (α i xx -α i zz ) 2 � + + 3 4 � (α i xy + α i yx ) 2 + (α i xz + α i zx ) 2 + (α i yz + α i zy ) 2 � , (5)
Consequently, the intensity of the Raman bands for parallel and crossed polarization can be written as:

I ⊥ i = (ω L -ω i ) 4 ω i [n(ω i ) + 1] γ 2 i 15 , (6)
I � i = (ω L -ω i ) 4 ω i [n(ω i ) + 1] 45ζ 2 i + 4γ 2 i 45 . (7)
The total Raman activity of powders (i.e., the case of unpolarized scattering) is a sum of both invariants

I tot i = I � i + I ⊥ i .
Finally, the Raman intensity is calculated according to Eq. ( 6) and (7), where each peak of the simulated spectra is obtained as a sum of Lorentzian shape lines constructed using the Raman activity of each vibrational mode with a fixed halfwidth. The codes and workflows that are used to achieve these calculations are available at [START_REF]Raman codes and workflows are[END_REF].

Following this procedure, we calculated the Raman spectra on glassy TeO 2 and glassy SiO 2 as an average of 10 spectra computed on 10 amorphous configurations that were sampled from the T=300K trajectory of the glass for each system. The results are presented in Fig. 1 and are compared to experimental results.

By inspecting the experimental Raman spectra of TeO 2 glass, one can identify several ranges with particular signatures that have been correlated in the literature to particular vibrational modes. In the low wavenumbers region, one can associate the main peak centered around 30 cm -1 with the Boson peak [33], while the shoulder around 140 cm -1 can be attributed to the intra-chain vibrations of Te-Te bonds [34]. In the central wavenumbers range, a broad peak coincides with the typical signatures of symmetric stretching vibrations in nearly symmetric (425 cm -1 ) and asymmetric (500 cm -1 ) Te-O-Te bridges [35,36]. The asymmetric stretching vibrations in nearly symmetric and asymmetric Te-O-Te bridges lead to typical signatures around 605 cm -1 and 660 cm -1 , respectively [35,36]. Refs. [START_REF] Zaki | [END_REF], [31], and [32], respectively. A 15% blue shift was applied to the Raman shift in the case of TeO 2 . The Raman intensity is normalized to the maximum band intensity.

Finally, the highest frequency band in the high wavenumber region can be attributed to the asymmetric stretching vibrations in asymmetric Te-O-Te bridges (710 cm -1 ) and asymmetric stretching of essentially covalent Te eq O bonds (770 cm -1 ) [35,36].

By following the calculated Raman spectrum (see Fig. 1), one can distinguish all of the observed features in the experimental spectra. Nevertheless, for wavenumbers higher than 300 cm -1 , all of the positions of the bands are underestimated by about 15%. This discrepancy is due to the small overestimation of the Te-O bond length in our models by about 0.04 Å [17]. This is a well-known issue and is related to the generalized gradient approximation family of exchange and correlation functionals, which slightly underestimate the lattice parameters and interatomic bond lengths. In addition, the large asymmetry around O atoms might influence the values of the cal- culated phonon frequencies [17]. Therefore, applying a 15% blueshift leads to modeled spectra that have an unprecedented good agreement with the measured spectra over the whole range of frequencies (see Fig. 1).

Similarly, the calculated Raman spectrum of glassy SiO 2 shows a good reproduction of the experimental spectrum. In this case, we note that there is a very limited blueshift thanks to the very well-defined structure of SiO 2 . Compared to TeO 2 , the Raman spectrum of glassy SiO 2 features a small peak at around 50 cm -1 that corresponds to the Boson peak. Given that the origins of the Boson peak in amorphous TeO 2 remain largely unknown, a comparative study of the Raman spectra between the TeO 2 and SiO 2 can be instructive. It is noteworthy that the simulated Raman spectra at zero temperature reveals no Boson peak feature for both oxides. At zero temperature, the Bose-Einstein term is negated and the pre-factor G simply acts as a damping function (1/ω) for the squared Raman tensor components. Interestingly, we observe that the "natural" Raman activity (i.e., squared Raman tensor components) is very small in the low frequency range (the highest intensity band in range [0-50] cm -1 is related to the highest intensity band in the whole spectra by 1:100). As temperature increases, the impact of the Bose-Einstein distribution function becomes stronger and one can find a significant response in the Raman spectra in the low frequency range. At T=100 K, one can clearly see that the Bose-Einstein distribution overpowers the damping function in the low frequency range, which leads to an "amplification" of the Raman spectra in the low frequency range (or, in other words, damping intensity and higher frequencies). Therefore, the pre-factor G takes a key role in the low frequency range Raman spectra. In general, at finite temperature, it acts as an amplifier for low frequency range. The strength of this amplification depends on a competition be- tween the two terms in the pre-factor G according to the temperature and the frequency of the considered Raman band.

When the nature of the low frequency Raman spectra amplification is established, it is worth comparing the spectra of TeO 2 and SiO 2 . At zero temperature, the Raman spectra of TeO 2 shows a steady growth of the intensity, starting rapidly from zero frequency. Actually, the lowest Raman-active vibration is found at 17 cm -1 . In the case of SiO 2 , the spectra shows a zero intensity at very low frequencies and only starts to grow weakly at approximately 50 cm -1 (the lowest frequency Raman-active vibration is found at 32 cm -1 ). At finite temperature, the band at low frequencies is clearly observed in the case of TeO 2 . Nevertheless, in the case of SiO 2 glass, even at T=300 K this band is not as strong as in the Raman spectra of TeO 2 . This can be explained by a frequency dependence of the prefactor G. This dependence follows a simple rule: the amplification of the Raman spectra is higher when the frequency of the band is lower. As mentioned earlier, the smallest phonon frequency of TeO 2 is found at 17 cm -1 , while the lowest frequency phonon in the case of silica glass is found at 32 cm -1 . This difference in ca. 15 cm -1 leads, in part, to a significant difference in the intensity of the Boson peak.

The advantage of the method that is presented in this letter is that it allows us to make a detailed analysis of the Raman bands because it gives access to the normal modes (atomic displacements of vibrational modes) of the Raman active vibrations. By analysing the structural fragments that are involved in vibrations related to the Boson peak, one infers that these modes are the result of collective vibrations of nanoclusters composed by structural units (TeO 4 , TeO 3 , TeO 5 ) as a whole. This result is consolidated by investigating the overall weight (impact) of the chemical species (Te and O) in a given Raman active mode, which is calculated as a ratio of the summary displacements made by atoms of a given chemical species with respect to the sum of displacements of all of the atoms in a given vibrational mode. In Figure S3, the weight for each chemical species is plotted as a function of the Raman frequency. It is noteworthy that the weight of O atoms will always be at least twice as much than the weight of the Te atoms because there are twice as many O atoms as Te atoms in TeO 2 glass. This weight proportionality holds true at the low frequency range, where the impact of Te and O atoms is nearly the same (taking a factor of 2 into account). In contrast, the weight of the oxygen atoms becomes dominant in the high frequency range, which indicates the key role of O vibrations in this range of frequencies. The low frequency vibration can be visualized thanks to the obtained eigenvectors of the dynamic matrix. An example of structural unit vibrations is plotted in Figure S4 in the supplementary materials. Understanding the nature of the frequency vibrations provides an explanation of the origin of the weak "natural" Raman activity that is observed at low frequency vibrational modes. The vibra-Copyright line will be provided by the publisher tions of structural units as a whole lead to a weak change of the polarizability of the whole system, and therefore to a small Raman response, which in turn is proportional to the derivative of the polarizability with respect to atomic displacements. However, the dominant Te-O vibrations in the middle and high frequency range, including bending, libration and stretching, leads to a significant impact on the polarizibility changes, and consequently to a high Raman activity, which is damped by prefactor G, as discussed earlier.

Finally, the origins of the strong Boson peak in TeO 2 , which were explained earlier by the low frequency vibrational states, compared to SiO 2 can be also attributed to the physical properties of both systems. In particular, assuming that the low frequency states are related to displacements of the whole structural units, one can put forward a correlation between their low frequency and the mass of the structural units involved in these vibrational states. In a simple harmonic approximation, the frequency is proportional to the inverse square root of the mass as ω = � k/m. Considering as a first approximation that ω is mostly affected by the mass of the structural units, and bearing in mind that the mass of the structural unit TeO 4 is almost twice that of SiO 4 , one can easily establish that the ratio of Boson peak frequency for TeO 4 and SiO 4 systems is equal to √ 2. The experimental frequency of Boson peak of TeO 2 is 35 cm -1 . By taking the ratio √ 2 into account, one obtains a frequency shift up to 49.5 cm -1 for SiO 2 amorphous silica. This is very close to the experimental data plotted in Figure 1b (see Table S1 in the SI for more details). This correlation is of a particular interest because it is only related to the mass of the structural units involved in the vibrational modes. Therefore, it is possible to predict the approximated position and the strength of the Boson peak in Raman spectra in other amorphous materials.

The Boson peak that is observed in the Raman spectra of amorphous TeO 2 (see Fig. 1) is significantly stronger than that in SiO 2 . This can be explained, in part, by the high concentration of phonon states in the low-frequency range. The density of phonon states (VDOS) is plotted in Fig. 3. The strong low frequency band in the VDOS of TeO 2 , when compared to the much weaker band in the VDOS of SiO 2 , rises rapidly and reaches a maximum at frequency ∼50 cm -1 . In contrast, the low frequency band in VDOS of SiO 2 is very broad and has a maximum at ∼130 -1 . It is noteworthy that the prefactor G acts as a damping factor at higher frequencies and, in addition to the lower density of states in the VDOS of SiO 2 , this is the reason why the Boson peak in the Raman spectra of SiO 2 is weaker than that of TeO 2 . These findings show that the intensity and the frequency of the Boson peak can be qualitatively estimated by VDOS analysis.

The Boson peak is widely interpreted as a consequence of the absence of long-range order and the frequency of the peak maximum is related to a medium-range order- size [31], which remains not very well defined. Indeed, an amorphous system can be considered as a crystal system with breaking of the periodic boundary conditions, which leads to many times Brillouin zone (BZ) folding. A good example to explain the appearance of low frequency peaks in Raman spectra due to BZ folding can be found in supercells [37]. Many times BZ folding leads to a complete disorder, which is accompanied by a degeneracy of whole BZ at the Γ point. In the case of ordered struc- tures, the frequency of acoustic vibrations near the Γpoint of the BZ is zero and starts to grow when k-vector shifts along the BZ boundary. The calculated phonon dispersion branches of the γ-TeO 2 crystal accompanied with VDOS compared to those of crystobalite α-SiO 2 are plotted in Figure S5 (see SI). One can find a sharp and strong peak in the VDOS of γ-TeO 2 at a frequency of ∼55 cm -1 , which is close to the Boson peak frequency in the experimental Raman spectrum. This band mostly consists of acoustic branches, which have a small dispersion due to low sound velocity in TeO 2 . In contrast, in the case of α-SiO 2 , the low frequency peak in VDOS is weak and broad in the low-frequency range, and has a maximum at ∼75 cm -1 . The acoustical phonon branches dispersion in α-SiO 2 is higher, and therefore the low frequency band in VDOS consists of a mixture of acoustical and optical phonon branches. This mixture triggers the formation of broad band from 0 up to 500 cm -1 with local maxima. The results of crystalline systems VDOS analysis are in line with those made on amorphous systems. Therefore, they can be used for qualitative analysis of the Boson peak nature. In addition, the dispersion of acoustic branches, which is related to sound velocity, is a key feature to predict the strength and position of the Boson peak. Our microscopic approach has shown that the Boson peak corresponds to collective displacements of small size nanoclusters composed from the glassy structural building blocs (e.g., TeO 4 ). Therefore, one can consider the mean Te-Te (Si-Si) distance as a Boson peak characteristic. In the case of TeO 2 , this distance corresponds to the mean shortest Te-Te (Si-Si) distance obtained in Te-

Copyright line will be provided by the publisher O-Te bridges, which is about 3.75 Å [15]. In the case of SiO 2 glass, the corresponding distance is equal to ∼3 Å (see Table S2 in the SI for more details). The minimum distance is responsible for the force constant k value in the simple harmonic approximation, namely the force constant is harder when the distance is shorter. This regularity explains why the Boson peak in case of silica has a higher frequency.

In summary, we present an efficient scheme to calculate the Raman spectra on amorphous periodic glassy oxide models containing several hundreds of atoms by combining the density functional theory and the finite differences method. Our scheme has the advantage of a quite reasonable and time efficient computational cost thanks to the efficient scalability of the CP2K code. The obtained Raman spectra in the case of TeO 2 and SiO 2 demonstrate a high level of agreement with the experiments. In addition, we study the origins of the Boson peak observed in glassy TeO 2 by focusing on two aspects. First, we investigate the reasons behind the strong Boson peak observed in glassy TeO 2 in close comparison to the less pronounced Boson peak reported in glassy SiO 2 . We show that the strength and position of the Boson peak is a function of the phonon density of the state's maximum (frequency and intensity). Furthermore, we find that as the frequency increases, the Boson peak is damped by inverse frequency and the Bose-Einstein distribution terms. Consequently, because the TeO 2 shows a lower frequency of the Boson band, it undergoes a higher amplification of its Raman spectra in the low frequency range. Second, we inspect the atomistic nature of the Boson peak by analysing the normal modes vectors and attribute the Boson peak to vibrations of Te-centered structural units as a whole. This leads to characterization of Boson peak by defining a medium-range ordersize that is equal to the first maximum in the partial cation-cation pair distribution function in oxide materials, which is very probably the case for the big family of binary (ternary) systems. Overall, we demonstrate that our computational scheme is able to compute accurate Raman spectra of large glassy systems, which paves the way towards a systematic computational study of the vibrational properties of oxide glasses. and crystobalite α-SiO 2 (bottom) calculated using supercell approach. The dynamical matrix was calculated for 4×4×3 supercell using finite difference method. The density of vibrational states (VDOS) is calculated as formally defined by equation: g(ω) = 1 N ν,q δ(ω -ω ν (q)), where N is the number of unit cells, ν is the band index and q is the q-point. The Dirac delta function δ(x) is substituted by Gaussian function with constant standard deviation. Appendix 1: Finite differences method According to Ref. 13 the Raman tensor of vibrational mode is defined as derivative of linear dielectric susceptibility χ with respect to normal mode ξ:

VDOS of crystalline oxides

Boson peak frequency and cation-cation distances

A ij = dχ ij dξ = α,β ∂χ ij ∂τ αβ u αβ , (1)
where u αβ = ∂τ αβ ∂ξ is a vector of atomic displacement of α'th atom along β Cartesian direction.

The calculation of derivative in equation ( 1) by finite difference method can be implemented in two general approaches. The first one corresponds to the distortion of relaxed structure according to normal modes, i.e. for each Raman active vibrational mode atomic positions are shifted by vector of atomic displacements ξ m with respect to equilibrium atomic positions r 0 in positive and negative direction, then for each disturbed structure the susceptibility tensor is calculated. Then the derivative for each vibrational mode is calculated as S6

Figure 1

 1 Figure1Calculated Raman spectra on glassy TeO 2 (a) and glassy SiO 2 (b) compared to experimental results from Refs.[START_REF] Zaki | [END_REF],[31], and[32], respectively. A 15% blue shift was applied to the Raman shift in the case of TeO 2 . The Raman intensity is normalized to the maximum band intensity.

Figure 2

 2 Figure 2 Simulated Raman spectra of glassy TeO 2 (a) and glassy SiO 2 (b) at different temperatures. The spectra are normalized to the highest frequency band, where the effect of pre-factor is minimal.

Figure 3

 3 Figure 3 Calculated phonon density of states of amorphous TeO 2 and SiO 2 .

Figure S2 :

 S2 Figure S2: Comparison of total X-ray and partial O-O, Si-O, Si-Si g(r) PDFs for various simulation sizes.

Figure S3 : 4 Figure S4 :

 S34S4 Figure S3: Frequency dependence of displacement weights of Te and O atoms in Raman active vibrational modes.

Figure S5 :

 S5 Figure S5: Phonon dispersion branches and vibratrional density of states of γ-TeO 2 (top)and crystobalite α-SiO 2 (bottom) calculated using supercell approach. The dynamical matrix was calculated for 4×4×3 supercell using finite difference method. The density of vibrational states (VDOS) is calculated as formally defined by equation: g(ω) = 1

Table S1 :

 S1 Boson peak frequency of common oxides. This work b Ref. 6 (Lowest frequency peak in vibrational density of states).

	Compound	Frequency, cm -1
	TeO 2	34 a	75 b	35 c
	SiO 2	65 a	50 d
	B 2 O 3	30 e
	GeO 2	40.5 f
	a

1 S5 c Ref. 7 d Ref. 8 e Ref. 9 f Ref. 10

Table S2 :

 S2 The most probable nearest neighbor distances in glassy TeO 2 and glassy SiO 2 .

	Property distance value,	Å
		TeO 2
	d(Te-Te)	3.75 a	3.7 b
	d(Te-O)	1.96 a	1.9 b
	d(O-O)	2.79 a	2.9 b
		SiO 2
	d(Si-Si)	3.1 c	3.2 d
	d(Si-O)	1.6 c	1.6 d
	d(O-O)	2.6 c	2.5 d
	a Ref. 11

b Ref. 6 c This work d Ref. 12

Copyright line will be provided by the publisher

Acknowledgements This work was supported by the ANR

via the projects TRAFIC (ANR-18-CE08-0016-01) and AMSES (ANR-20-CE08-0021) and by région Nouvelle Aquitaine via the F2MH project AAP NA 2019-1R1M01. This work was also supported by institutional grants from the French National Research Agency under the Investments for the future program with the reference ANR-10-LABX-0074-01 Sigma-LIM. Calculations were performed by using resources from Grand Equipement National de Calcul Intensif (GENCI, grants No. A0080910832, A0080910832, A0100910832 and A0120910832). We are grateful to TGCC for the generous cpu allocation within the special session project spe00019. We used computational resources provided by the computing facilities Mésocentre de Calcul Intensif Aquitain (MCIA) of the Université de Bordeaux and of the Université de Pau et des Pays de l'Adour.

Supporting Information

Ab initio study of the Raman spectra of amorphous oxides:

Insights into the Boson peak nature in glassy TeO 2 Evgenii M. Roginskii, * , † Raghvender Raghvender, ‡ Olivier Noguera, ‡ Philippe Thomas, ‡ Olivier Masson, ‡ and Assil Bouzid * , ‡ †Ioffe Institute, Polytekhnicheskaya 26, 194021 St. Petersburg, Russia ‡Institut de Recherche sur les Céramiques (IRCER) -UMR CNRS 7315, Université de Limoges, Centre Européen de la Céramique, 12 rue Atlantis, 87068 Limoges Cedex, France E-mail: e.roginskii@mail.ioffe.ru; assil.bouzid@unilim.fr

DFT setup

The DFT calculations presented in this work, were carried out within generalized gradient approximation to density functional theory (DFT) as implemented in CP2K software package. 1 The Perdew-Burke-Ernzerhof (PBE) 2 functional was used in all the calculations.

The Gaussian Plane Waves method 3 was used to form wavefunctions basis set with cutoff energy equal to E g = 1000 Ry for planewave part and using triple-zeta valence basis for local combination of molecular orbitals part. Goedecker-Teter-Hutter pseudopotentials are used to describe core-valence interactions. 4 The self consistent calculations were performed with electronic convergence criteria of 10 -8 Ha. Geometry optimization is performed by relaxing S1 Draft version, January 19, 2023 both cell parameters and atomic positions at fixed angles for a given atomistic configuration.

During the relaxation, atomic displacements and lattice parameters were varied until forces on atoms become below 10 -5 Ha/Bohr and pressure below 10 -2 GPa. Dynamical properties of relaxed structures were carried out using finite difference method.

SiO 2 Models

We have performed ab-initio MD to obtained a structural model of SiO 2 glass with 270 atoms at a constant density of 2.2 gm/cm 3 . The starting ab-initio model was produced after a classical MD run at T = 2000 K for 1000 ps using BKS type potential with potential parameters taken from Ref. 5

The second approach is devoted to calculations of derivatives with respect to atom α displacements along Cartesian directions τ β by finite difference method followed by multiplication on vector u αβ (right part of equation ( 1)).

The choice of the method depends on symmetry of the structure under study.