Javier Navarro-Ruiz

Camila Rivera-Cárcamo

Bruno Machado

Philippe Serp

Iker Del Rosal

Iann C Gerber

Computational Design of Pd Nanocluster and Pd Single-Atom Catalysts Supported on O-Functionalized Graphene

ABSTRACT: A crucial step in the preparation of supported metal catalysts is related to the choice of support, since it regulates the anchoring of the metal species in an environment of strong interaction. In this sense, the binding on the support, which may present defects, as well as induced modifications of the electronic structure, help to determine efficient metal/support combinations for catalysis applications. In this work, first-principles studies have been carried out to model and describe the geometric and electronic properties of an O-functionalized graphene as a model carbon support, in addition to single Pd atom and Pd 13 nanocluster (experimentally observed on various carbon supports) as supported metal catalysts. The carbon-based nanomaterial includes experimentally probed abundant oxygen functional groups and point defects, demonstrating high thermal stability at room temperature. The characterization confirms the structural motifs presented in the model, distinguishing the organic functionalities and various defects by spectroscopy techniques. Interestingly, the oxygenated carbon-based support presents a metallic character when the formal charges of the embedded Pd atom and nanocluster have been determined to be [Pd] + and [Pd 13 ] 3+ . The hydrogenation process of the adsorbed cluster is also presented, with the determination of a ratio between the adsorbed hydrides and the surface metal atoms close to unity at room temperature, thanks to ab initio molecular dynamics. KEYWORDS: Palladium, nanocluster, single atom, supported catalyst, oxygen functional group, graphene, hydrogenation 1. INTRODUCTION Supported-metal catalysts are among the most important materials in heterogeneous catalysis. [START_REF] Ma | Encyclopedia of Inorganic Chemistry[END_REF] Consisting of single atoms (SA), clusters, and/or nanoparticles (NP) anchored onto solid supports, these materials can maintain the advantages of molecular catalysts (often highly active/selective) and minimize their disadvantages (expensive catalyst separation from products and corrosionrelated damage), thus allowing their use in industrial processes. [START_REF]Applied Homogeneous Catalysis with Organometallic Compounds[END_REF] To do this, the system under consideration (usually metal NP) must be in a sufficiently high dispersed form. [START_REF] Rothenberg | Heterogeneous Catalysis[END_REF] The active metal component is usually deposited on the surface of a porous solid support, a thermally/mechanically stable material, which is able to maintain the catalyst life, but also to disperse the metal and prevent its agglomeration. [START_REF] Rothenberg | Heterogeneous Catalysis[END_REF] The catalyst supports can be inert or active in the catalytic reactions and can cooperate with the metal, for example via spillover. [START_REF] Xiong | Spillover in Heterogeneous Catalysis: New Insights and Opportunities[END_REF] Importantly, electronic and geometric effects induced by the support contribute to control catalyst performances (activity, selectivity, and stability), especially in highly dispersed catalysts, for which the metalsupport interaction is decisive. [START_REF] Van Deelen | Control of Metal-Support Interactions in Heterogeneous Catalysts to Enhance Activity and Selectivity[END_REF][START_REF] Gerber | A Theory/Experience Description of Support Effects in Carbon-Supported Catalysts[END_REF] Typical support materials include metal oxides, zeolites, and various carbon nanomaterials. [START_REF] Serp | Comprehensive Inorganic Chemistry II[END_REF] Carbon-supported metal catalysts are employed in a large number of industrial reactions including oil refining, hydrogenation, and oxidation, as well as for electrocatalysis. [START_REF] Xi | [END_REF] Several carbon-based (nano)materials can be used as catalyst support: (i) 0D (NP); (ii) 1D (carbon nanotubes, CNTs); (iii) 2D (graphene sheets); and (iv) 3D (mesoporous carbons, MC). 9

Due to the unique lattice structure and physicochemical properties of (few-layer) graphene, this material has been particularly studied as catalyst support, 10 both for metal NP and metal SA. 11 Nevertheless, most of the produced graphene sheets are not perfect, presenting defects (intrinsic or due to the presence of a foreign atom) that influence their physical, mechanical, and even electronic properties. Besides the intrinsic defects [Stone-Wales (SW) and vacancies], 12 chemical doping can modify graphene properties and induce specific sites, which contribute to different catalytic performance. Among the various potential dopants, oxygen is considered as an excellent element for chemical doping of carbon-based (nano)materials. These O-doped materials are usually produced with strong oxidizing agents in acidic media to introduce oxygenated functional groups (OFG) comprising carboxylic acids, esters, ketones, phenols, lactones, lactols, or quinones. 13 O-functionalized graphene can be represented with a simplified model consisting of individual graphene sheet (sp 2 carbon atoms) 14 decorated with OFG on both the basal planes and edges, 15 resulting in a carbon network having both sp 2 -and sp [START_REF] Rothenberg | Heterogeneous Catalysis[END_REF] -hybridized carbons bearing OFG. 16 However, currently the precise atomic and electronic structure of O-functionalized graphene sheets remains largely unknown. This is mainly due to the uncertainty related to the amount and types of OFG, [START_REF] Compton | Graphene Oxide, Highly Reduced Graphene Oxide, and Graphene: Versatile Building Blocks for Carbon-Based Materials[END_REF] and their nonstoichiometric atomic composition. Owing to these OFG, the materials exhibit an amphoteric (acid-base) character that improves metal adsorption and dispersion, resulting in mixtures of metal NP and SA, which have been observed on commercial catalysts. [START_REF] Petek | Nanoparticles and Single Atoms in Commercial Carbon-Supported Platinum-Group Metal Catalysts[END_REF] When such catalysts are used for hydrogenation reactions, hydrogen spillover can also operate. It is also promoted by the presence of OFG on the carbon surface, and examples of cooperative catalysis between metal NP and metal SA thanks to H-spillover on the carbon support have been reported. [START_REF] Serp | Cooperativity in Supported Metal Single Atom Catalysis[END_REF] The understanding of the chemical reactivity of such complex systems definitively requires robust models.

Theoretical simulations can provide remarkable insight into the possible thermodynamic and kinetic mechanisms required for a greater understanding on the structure of OFG found on carbon nanomaterials and on their influence on the reactivity of metal species. [START_REF] Xi | Confined-Interface-Directed Synthesis of Palladium Single-Atom Catalysts on Graphene/Amorphous Carbon[END_REF] Nowadays, endeavors are directed towards simple schemes that allow for basic structural motifs using common oxygen groups, specifically hydroxyl and 1,2-ethers (epoxide). [START_REF] Zhou | Density Functional Theory Modeling of Multilayer "Epitaxial" Graphene Oxide[END_REF][START_REF] Guilhon | Chemically Tunable Properties of Graphene Covered Simultaneously with Hydroxyl and Epoxy Groups[END_REF] Thus, the aim of the present work is: i) to design, based on ab initio electronic structure methods, a more realistic Ofunctionalized support mainly based on experimental data; and ii) to provide a detailed description of different physicochemical features of the modified support and suggested supported palladium catalysts (Pd/C is a very common hydrogenation catalyst). In particular, analysis of the structural (average bond length, mean coordination number, and mean interatomic distances), electronic (density of states, crystal orbital Hamilton population, and Mulliken population analysis), and vibrational (infrared and Raman spectra) properties are presented. Special attention to the differences between single-atom and cluster-containing catalysts is also addressed.

COMPUTATIONAL DETAILS

Periodic density functional theory (DFT) calculations were performed using the ab initio planewave pseudopotential approach as implemented in the Vienna Ab initio Simulation Package (VASP; version 5.4). [START_REF] Kresse | Efficiency of Ab-Initio Total Energy Calculations for Metals and Semiconductors Using a Plane-Wave Basis Set[END_REF][START_REF] Kresse | Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set[END_REF] The Perdew-Burke-Ernzerhof 25 exchange-correlation functional within the generalized gradient approximation was chosen, and van der Waals interactions were taken into account through the D3 correction method of Grimme et al. [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] The innermost electrons were replaced by a projector-augmented wave (PAW) approach, [START_REF] Blöchl | Projector Augmented-Wave Method[END_REF][START_REF] Kresse | From Ultrasoft Pseudopotentials to the Projector Augmented-Wave Method[END_REF] while the valence monoelectronic states were expanded in a plane-wave basis set with a cut-off energy of 450 eV. Partial occupancies were estimated with a Gaussian smearing (σ) of 0.05 eV during all relaxations and extrapolating the energies to σ = 0.00 eV, except for the determination of the electronic density of states (DOS) for which the tetrahedron method with Blöchl corrections was used. The carbon support was modelled as a graphene slab, showing that a (7 × 7) supercell was found to be large enough to include a broad representation of experimentally observed OFG. The (7 × 7) unit cell contains 98 C atoms and a Γ-centered (3 × 3 × 1) k-point grid generated using the Monkhorst-Pack method [START_REF] Monkhorst | Special Points for Brillouin-Zone Integrations[END_REF] was employed as a good compromise between accuracy and computational cost. A higher density of k-points and energy cut-off were used in obtaining DOS calculations for all the structures considered for the present study. All of the systems considered are neutral. Spurious interactions between the modelled slab and its periodic images perpendicular to the slab were eliminated by introducing a vacuum region by at least 10 Å, and by applying a dipole correction to the total energy along the z-direction. Iterative relaxation of atomic positions proceeded until the change in total energy between successive steps was less than 10 -6 eV per cell and the residual forces on any direction acting on the atoms (evaluated using the conjugate-gradient algorithm) were less than 0.015 eV Å -1 .

To probe the thermal stability of some of the systems, first-principle Born-Oppenheimer ab initio molecular dynamics simulations (AIMD) were performed within the canonical NVT ensemble (that is, constant number of particles, the system's volume, and the temperature). Starting from the optimized structure, each slab was propagated at 298 K for 20 picoseconds at time steps of 0.5 femtoseconds for integrating the atomic equations of motion. The level of theory was the same as described above and the k-mesh was reduced to a gamma point. In all cases, the electronic structures were well converged to 10 -6 eV. The temperature was regulated by the Nosé thermostat. [START_REF] Nosé | A Molecular Dynamics Method for Simulations in the Canonical Ensemble[END_REF] To provide additional insight on the modelled support, as well as the interactions between the metal catalysts and the support itself, the Crystal Orbital Hamilton Population (COHP) analysis [START_REF] Dronskowski | Crystal Orbital Hamilton Populations (COHP): Energy-Resolved Visualization of Chemical Bonding in Solids Based on Density-Functional Calculations[END_REF] is performed with the LOBSTER package. [START_REF] Deringer | Crystal Orbital Hamilton Population (COHP) Analysis as Projected from Plane-Wave Basis Sets[END_REF][START_REF] Maintz | Analytic Projection from Plane-Wave and PAW Wavefunctions and Application to Chemical-Bonding Analysis in Solids[END_REF] Based on contributions from orbital pairs, COHP is an energy division of the band structure, so it acts on a local basis. In this regard, from a plane-wave DFT output, the LOBSTER package allows calculating the projected COHP in a local atomic base set (pCOHP) as well as a consistent density of states projected by atoms (pDOS). The COHP analysis is performed for pairs of atoms at distances according to the nearest neighbor, and subsequently the COHP is plotted by accumulating all the pairs with the same elements. To demonstrate the contributions of different elements, the pDOS and pCOHP plots are shown together. The charge spilling (a measure that assesses the quality of the projection) for all the compounds of this work reaches a maximum value of 2%. Furthermore, for particular cases the integrated pCOHP (IpCOHP) was calculated, providing qualitative data on the average bond strength in covalent entities.

Vibrational frequencies of the ground-state geometries were determined from phonon calculations by diagonalizing the numerical Hessian matrix (VASP calculates phonons at the zonecenter -Γ-point-within the supercell approach). The elements of the Hessian matrix (harmonic approximation) were calculated by finite differences of the analytical first derivatives of the energy with respect to atomic displacements using the central-difference formula, in which each atom is displaced using small steps of 0.015 Å in both positive and negative directions of each coordinate to estimate the harmonic potentials. The value of the infrared (IR) intensities for each normal mode was also computed by performing linear response calculations within density-functional perturbation theory (DFPT) formalism. From the Born effective charge (BEC) tensor, which refers to the change in polarizability of atoms with respect to an external electric field, the vibrational intensities of the modes were calculated using the formula by Giannozzi and Baroni. [START_REF] Giannozzi | Vibrational and Dielectric Properties of C60 from Density-Functional Perturbation Theory[END_REF] Additionally, for both the pristine and O-modified supports, the off-resonance Raman activity is calculated for each of the normal modes obtained, taking advantage of the derivative of the polarizability (or macroscopic dielectric tensor) with respect to the coordinates of the normal modes. [START_REF] Porezag | Infrared Intensities and Raman-Scattering Activities within Density-Functional Theory[END_REF] Therefore, from the phonon calculation at Γ-point and the macroscopic dielectric tensor of the previous DFPT, the non-resonant Raman scattering activities are computed.

RESULTS AND DISCUSSION

This section is organized as follows. First, results based on experimental data, related to the computational modelling of the oxygen-functionalized support from pristine graphene (G) are shown. Then, analyses on the structure-property, electronic structure, and vibrational of the support are provided from the density of states, partial atomic charges, and the IR and Raman spectra. Afterward, the anchoring of both a palladium single atom (Pd SA ) and a Pd cluster (Pd NP ) is discussed, paying special attention to the structural and electronic changes that occur both in the carbon support and in the metallic systems before and after metal deposition. Finally, the hydrogenation of the ultra-small Pd NP is considered, discussing how this process affects the nature of the catalyst (the overall structure and chemical composition is provided in Scheme 1).

Design of the support. The first step in the development of supported metal catalysts is the design of the support. To this end, graphene was used as model carbon support (Figure 1a), given its excellent properties shown in supported catalysis. [START_REF] Gerber | A Theory/Experience Description of Support Effects in Carbon-Supported Catalysts[END_REF] Formation and reconstruction of vacancy defects. Lattice imperfections are unavoidably introduced into graphene during its growth or processing. SW defect, [START_REF] Stone | Theoretical Studies of Icosahedral C60 and Some Related Species[END_REF] single and multiple vacancies (MV), [START_REF] Krasheninnikov | Bending the Rules: Contrasting Vacancy Energetics and Migration in Graphite and Carbon Nanotubes[END_REF] carbon, foreign adatoms, [START_REF] Lehtinen | Magnetic Properties and Diffusion of Adatoms on a Graphene Sheet[END_REF] and substitutional impurities 39 are typical stable graphene point defects that have been predicted. As a result, carbon atoms with a sp 2 hybridization could undergo patterns other than hexagonal, leading to different systems. From a thermodynamic point of view, double vacancies (DV) appear as the most favored vacancy defects in graphene over single vacancies (SV). [START_REF] Krasheninnikov | Bending the Rules: Contrasting Vacancy Energetics and Migration in Graphite and Carbon Nanotubes[END_REF] The reconstruction in the former's atomic network allow a coherent defective lattice without under-coordinated atoms, and therefore with no dangling bonds. Furthermore, the removal of more than two atoms may lead to larger and more complex defect configurations.

As a consequence of the combination of two SV or elimination of two contiguous atoms, a DV can be formed, in which completely reconstructed DV with a pair of pentagonal and a single octagonal cells of carbon atoms [V 2 (5-8-5) defect] appears. However, it is further reconstructed and transformed into a new defect structure by a SW type transformation [V 2 (555-777) defect], consisting of three pentagonal and three heptagonal cells of carbon atoms in the graphene lattice. [START_REF] Lee | Coalescence, and Reconstruction of Vacancy Defects in Graphene Layers[END_REF] From an energy point of view, the V 2 (555-777) defect is much more stable than two isolated SV, its formation energy being lower than in the case of V 2 (5-8-5), of the order of 1 eV. This type of defect is also commonly noted in electron microscopy analyses. [START_REF] Kotakoski | Stone-Wales-Type Transformations in Carbon Nanostructures Driven by Electron Irradiation[END_REF] For all these reasons, the introduction of two V 2 (555-777) defects were firstly considered on pristine G, as revealed in Figure 1b. These defects did not cause any deformation in the surface plane.

Addition of oxygen functional groups (OFG).

The reactivity of the carbon support is increased locally due to the presence of these point defects and allow surface functionalization.

Therefore, the next step was to include OFG that can be experimentally detected and quantified. [START_REF] Figueiredo | Functionalization of Porous Carbons for Catalytic Applications[END_REF] For this purpose, Temperature Programmed Desorption (TPD) experiments in conjunction with Mass Spectroscopy (MS) measurements can be applied. [START_REF] Romanos | Controlling and Quantifying Oxygen Functionalities on Hydrothermally and Thermally Treated Single-Wall Carbon Nanotubes[END_REF] TPD-MS study is certainly a frequent technique for the description of the OFG on carbon materials, decomposing into CO/CO 2 and, on certain occasions, H 2 O at different temperatures. In our case, multi-walled carbon nanotubes (MWCNT) or few-layer graphene (FLG), which have been oxidized with nitric acid, were used as support. This procedure is commonly used on carbon supports prior to metal deposition to provide anchoring sites (OFG) for the metallic phase. [START_REF] Gerber | Theoretical and Experimental Studies on the Carbon-Nanotube Surface Oxidation by Nitric Acid: Interplay between Functionalization and Vacancy Enlargement[END_REF] Full experimental details are described in the ESI and the progress of CO 2 and CO throughout TPD-MS studies under helium for the oxidized MWCNT (CNT ox ) are shown on Figure S1. Taking into account the spectra of CO 2 and CO, it can be seen that the CNT ox support comprises carboxylic (CO 2 at 217/291 ºC, peaks #1 and #2), carboxylic anhydride (CO 2 at 437 ºC, peak #3 and CO at 491 ºC, peak #1'), lactone (CO 2 at 576 ºC, peak #4), phenol (CO at 717 ºC, peak #2'), and carbonyl/quinone (CO at 882 ºC, peak #3') surface groups (the deconvolution of CO 2 and CO profiles of CNT ox using a multiple Gaussian function are shown in Figure S1 and the quantification in Table S1). Considering the specific surface area of CNT ox (S BET = 220 m 2 g -1 ) and the fact that the hexagonal unit cell of graphene contains two carbon atoms having an area of 0.052 nm 2 , so 38.5 C nm -2 , it was possible to calculate the density of OFG, as illustrated in Figure 2. For comparative purpose, the same procedure was followed for an oxidized FLG support. The comparison of the density of OFG obtained on CNT ox and FLG ox supports is shown on Figure S2, observing the same trends for both supports. It is logic to find few differences since the two supports do not present the same structure and should not have exactly the same reactivity towards the oxidant (HNO 3 ).

In order to take into account each of the OFG, the modelling of the different groups was carried out by means of a sequence of serial optimizations on the defective graphene. The addition is made by functional groups, so that the support gradually adapts against oxidation. Firstly, the carboxylic acid groups are introduced. Since a carbon edge model was never considered, the acid function is irremediably located on top of the surface. Thus, different positions were tested (that is, on pentagonal, hexagonal and heptagonal cell carbons), the most stable being the one in which the -COOH motif is adsorbed on a carbon belonging to a five-membered ring. Due to the large number of motifs to take into account according to the high density, they were distributed in several successive sets (relaxing all from one to another). Secondly, the carbonyl groups were added in an analogous way (distribution in successive relaxations). Contrary to the acidic functional group, the carbonyl groups must be located at the edges. As mentioned, the support model does not explicitly consider carbon atoms at the edge of graphene, such as zigzag or armchair edges, which necessarily causes an increase in the size of the V 2 (555-777) defects, going from two DV to two larger vacancies commonly known as MV.

Third and last, the rest of the functional groups of lower density were included at the same time, that is, phenol/ether, and lactone. For reasons of steric hindrance, the anhydride functional group was not taken into account during the design of the functionalized support. These last three functional groups must also be located on the edges of the graphene support. The final result of which can be seen in the left inset of Figure 3 (top and side views, a). Unsurprisingly, the largest MV (center of the unit cell) predominates in the support, although completely saturated by carbonyl and quinone groups. This high concentration of oxygenated groups in the MV induces the loss of the graphene planarity due to the electronic repulsion between the lone pairs of oxygen atoms, leading to the formation of surface corrugations. Finally, many experimental carbon-based catalysts use FLG and/or CNT as support. Thus, in order to validate the simplified support modelled by a single layer, a rhombohedral (ABC stacking order) trilayer graphene was carried out as a more realistic system (the atomistic structure is depicted in Figure S3). No significant structural or electronic changes were appreciated, apart from the new existing π-π stacking between layers and the lifting of some OFG that pointed towards the nearest graphene layer.

Therefore, the single layer model appears to be a realistic model of carbon material as a support.

Thermal stability. The thermal stability of the optimized O-functionalized graphene (G ox ) was investigated through an AIMD simulation at a temperature of 298 K, which can guide the design of carbon support with higher stability. With the initial geometry of the last relaxed structure, the simulation analysis revealed a very fast energy convergence over time, as shown in Figure 3. The final geometry of the simulation after relaxation (top and side views, right inset, Figure 3b) is very similar to the geometry of the starting point and the observed structural changes are small, essentially rotations of carboxylic groups and out-of-plane deformations in the graphitic region, suggesting that at room temperature the designed G ox support is stable.

Characterization of the support. Support characterization plays an important role in the process of catalyst design, and many techniques are now available for a complete characterization of the structure. To evaluate the validity of the computational model, analysis of the electronic and vibrational structure has been performed for both G and G ox . The electronic behavior of graphene can be clearly seen in Figure 4a. Around the Fermi energy, pDOS exhibited a bandgap energy close to zero, characteristic of gapless semimetals. In contrast, G ox differed from this behavior and presented a metallic character, as shown in Figure 4b, since it has a conduction band that is partially filled. The electronic band structure of these type of materials is not yet clearly understood due to nanoscale inhomogeneities in the structure, but optical studies indicate that it can be tuned by controlling the degree of functionalization. [START_REF] Mkhoyan | Atomic and Electronic Structure of Graphene-Oxide[END_REF][START_REF] Luo | Photoluminescence and Band Gap Modulation in Graphene Oxide[END_REF] On the other hand, more details can be obtained by analyzing the COHP function, [START_REF] Dronskowski | Crystal Orbital Hamilton Populations (COHP): Energy-Resolved Visualization of Chemical Bonding in Solids Based on Density-Functional Calculations[END_REF] which provides an energy-based conception of chemical bonding.

The COHP(ε) profiles for the nearest-neighbor contacts in graphene and G ox are also plotted in Figure 4a and Figure 4b, respectively. The predominant interactions of both graphene correspond to C-C interactions, which showed a bonding character of the states lying up to the Fermi level (π molecular orbitals, MO), followed by an antibonding character of the conduction states (π* MO).

However, in the case of G ox , antibonding states appeared even below the Fermi level, around -6 eV, mainly due to C-O interactions. The IpCOHP, a bond strength index, is 10 eV per C-C interaction in the case of pristine graphene, and 7-12, 11-17, and 8 eV per C-C, C-O, and O-H interaction, respectively, in the case of G ox .

The calculation of effective atomic charges plays an important role in understanding the electrostatic interactions that are established in a system. By Mulliken population analysis obtained by integrating the pDOS, the corresponding pMPA charges can be plotted. A representation of the charge distribution in both materials is found in Figure 4. While pristine G (Figure 4c) clearly exhibits no electronic variation between all carbon atoms, G ox (Figure 4d) expressed a wide variety of charge ranging from positively charged in the case of carbon and hydrogen atoms attached to oxygen atoms, to negatively charged in the case of the oxygen atoms. This range of electrostatic interactions concentrated under the same material certainly opens up the possibility of using G ox in electron exchange applications. 15 Furthermore, the identification of the OFG and the chemical composition in the G ox material is corroborated by means of IR (Figure 4e) and Raman (Figure 4f) spectroscopy. The IR spectrum of G mainly presented a single intense band at 1526 cm -1 (stretching of C=C aromatic bonds), while G ox presented a more complex IR spectrum. A detailed analysis on the vibrational modes of the most significant IR bands is collected in Table 1. The typical broad band between 2400-3400 cm -1 associated to the stretching vibration of O-H groups is not present in the oxidized material. Two explanations appear as plausible for this mode to appear inactive: the absence of solvent and, therefore, of possible H bonds with the acid group that highlight the IR signal, or the harmonic approximation used to determine the vibrational frequencies, since the presence of light atoms such as hydrogen causes the incorporation of anharmonic corrections to be almost essential for a proper description. The bands located at 1750 and 1540 cm -1 have been assigned to stretching vibration of C=O carboxyl/lactone and carbonyl groups, respectively. Compared with G spectrum, the aromatic C=C band showed a red-shift centered at around 1480 cm -1 .

Finally, the presence of C-O bonds is confirmed, and the corresponding IR bands can be observed at around 1260, 1120 and 1010 cm -1 , which have been assigned to vibrations from ether, carboxyl and lactone functional groups, respectively. In overall, simulated IR spectra reproduced rather well the characteristic IR band positions previously assigned in the literature [START_REF] Si | Synthesis of Water Soluble Graphene[END_REF][START_REF] Bagri | Structural Evolution During the Reduction of Chemically Derived Graphene Oxide[END_REF] and the ones obtained for CNT ox (see Figure S4). It should be noted that, while the O-H hydroxyl band appeared at relatively lower values, both the C=O carbonyl and C=C aromatic bands appeared at slightly higher values, in this case the C=O from carboxyl and globally all the C-O bands of the OFG agreeing relatively better.

The Raman spectrum of G (Figure 4f) displayed no disorder-induced D-band, as no defects are present, and a narrow G-band originated from the bond stretching of sp 2 atoms at 1577 cm -1 , in line with the value of 1585 cm -1 from experiments in single graphene sheets, [START_REF] Graf | Spatially Resolved Raman Spectroscopy of Single-and Few-Layer Graphene[END_REF] and 1580 cm -1 in FLG (Figure S5). Importantly, the inclusion of OFG (G ox in Figure 4f) caused a broadening of the G-band, leading to a series of new high-intensity bands red-shifted with respect to the G-band of pristine graphene, located in a range similar to the broad D-band in the experimental data (around Lactone a IR bands refer to those bands shown in Figure 4c.

b See Figure S4 for the IR spectrum.

1300 cm -1 ). Nonetheless, due to the closeness in the frequencies of the sp 3 carbon Raman-active band (-COOH moieties attached to graphitic areas) and the D-band of graphitic carbons, both types of bands could merge under the same broad spectral feature. Lastly, in structures with welldefined double bonds, active blue-shifted Raman bands can be detected, although their intensity is usually low. [START_REF] Kudin | Raman Spectra of Graphite Oxide and Functionalized Graphene Sheets[END_REF] In the case of G ox , there are two Raman peaks above the G-band of graphene, one of higher intensity at 1623 cm -1 , and the other at 1752 cm -1 . After inspecting the atoms involved in these normal modes, these vibrations come from the displacements of the carbon atoms in fiveand seven-membered rings simultaneously, as a consequence of the reconstructed DV. Therefore, the defects can also induce blue shifts. Such spectrum modifications are not visible in the case of FLG ox or CNT ox (see Figure S5), probably because of the large number of graphene layers (around 10) in these materials.

Anchorage of Pd SA and of Pd 13 cluster. Once the pristine graphene has been functionalized, it can be used as support for metal deposition. In the case of Pd/C catalysts, it has been shown that both Pd NP and Pd SA can coexist on different types of carbon materials. Thus, mixtures of Pd NP and Pd SA were observed in commercial Pd/C catalysts prepared on activated carbon, [START_REF] Petek | Nanoparticles and Single Atoms in Commercial Carbon-Supported Platinum-Group Metal Catalysts[END_REF] but also for Pd/CNT and Pd/FLG catalysts, [START_REF] Rivera-Cárcamo | Control of the Single Atom/Nanoparticle Ratio in Pd/C Catalysts to Optimize the Cooperative Hydrogenation of Alkenes[END_REF] as revealed in Erreur ! Source du renvoi introuvable.. The surface oxygen contents determined by XPS for the CNT ox and FLG ox supports used for Pd deposition were 5.7 and 3.8 at. %, respectively. [START_REF] Contreras | Effect of Mesoporous Carbon Support Nature and Pretreatments on Palladium Loading, Dispersion and Apparent Catalytic Activity in Hydrogenation of Myrcene[END_REF] The Pd particle size distributions (particle size distribution based on total particles and on total atoms) of these three catalysts are given in Figure S6. The mean Pd NP size in these samples varies between 0.9 (Pd/CNT) and 2.1 nm (Pd/FLG), and the percentage of Pd SA between 48% (Pd/CNT) and 1.8% (Pd/FLG).

For the model catalyst, the inclusion of first a Pd SA and then a small Pd NP (i.e., less than 1 nm, similar to the ones experimentally observed particularly in the Pd/CNT sample) was carried out, to finally encompass both species on the same support. As the atomic radii of transition metals are larger than the one of carbon, metal adatoms are generally located above the graphene plane, a phenomenon that is more pronounced on SV than on DV. [START_REF] Krasheninnikov | Embedding Transition-Metal Atoms in Graphene: Structure, Bonding, and Magnetism[END_REF] In this regard, DV are suitable for metal adatom adsorption, which leads to thermodynamically more favorable situations. [START_REF] Rivera-Cárcamo | Single Atom Catalysts on Carbon-Based Materials[END_REF] This greater structural stability makes them less reactive (since they are less exposed on the surface) than the SA stabilized by SV. Therefore, the Pd atom adsorbed on an SV of the G ox as single-atom catalyst (SAC) was chosen in this study, as can be seen in Figure 6a. The SAC is thus modelled by removing a carbon atom in the graphene lattice, and subsequently adsorbing a Pd atom. It is found that the vacancy formation energy on G ox is of about 5.6 eV, whereas according to the literature around 7.5 eV is needed to form a SV in graphene. [START_REF] Skowron | Energetics of Atomic Scale Structure Changes in Graphene[END_REF] Thus, the removal of a carbon atom from the graphene lattice appears to be slightly more favorable from a thermodynamic view-point in the case of G ox . Different configurations were tested, the most stable being the one in which the carbon belonged to a six-membered ring of a graphitic region and far from any oxygen functionality. The metal anchored on support forms three covalent bonds with the neighboring carbons in a Pd-C 3 coordination structure (average bond length of 1.95 Å) and is displaced outward by 1.01 Å from the graphitic surface, probably due to the intrinsic corrugation of the support itself. [START_REF] Krasheninnikov | Embedding Transition-Metal Atoms in Graphene: Structure, Bonding, and Magnetism[END_REF] The binding energy, calculated as the difference between the SAC and the energy of the SV oxygenated support plus the energy of the isolated Pd, is 4.1 eV, pointing towards a favorable energy process. 4d 10 5s 0 . However, when the metal is anchored to the SV, the computed average occupation of the bands became [Kr] 4d 8.8 5s 0.2 , undergoing an oxidation process (giving an electron to the carbon support) that makes it formally have an atomic charge of +1, corroborated by the intense blue color in the corresponding pMPA charges of SAC depicted in Figure 6c. This presents the possibility of addressing a Pd + -based chemical reactivity. In the particular case of alkene hydrogenation reactions, having a Pd + can strengthen the Pd-alkene interaction resulting in less active catalysts. [START_REF] Huang | Toward Understanding of the Support Effect on Pd1 Single-Atom-Catalyzed Hydrogenation Reactions[END_REF] Experimentally, it was shown that the X-ray Adsorption Near Edge Structure (XANES) spectrum of Pd SA /graphene is close to those of PdO, which implies that the Pd SA on graphene surface were mostly present in a +2-valence state. [START_REF] Huang | Toward Understanding of the Support Effect on Pd1 Single-Atom-Catalyzed Hydrogenation Reactions[END_REF] In that case, the first shell peak for Pd SA /graphene was contributed to a mixture of Pd-O and Pd-C coordination.

As expected, the SAC IR profile is very similar to that of the support (represented in Figure 6d, for a detailed analysis on the assignment of the IR spectrum of the SAC see Table S2). The changes appear in the regions associated with the metal atom. Because of the presence of the newly formed Pd-C bonds, a new band assigned to the ν(Pd-C) appeared; in particular, the band around 920 cm -1 . Pd SA anchoring caused the IR band corresponding to ν(C=C) to be displaced towards smaller values, going from 1477 to 1395 cm -1 . In view of these results, IR spectroscopy can be considered as a useful technique for the study and identification of metal-carbon bonds during the preparation of SAC.

Regarding the Pd NP , a cluster model of 13 atoms is considered, which has been proven to be robust and enough efficient for describing theoretically the hydrogenation of alkenes. [START_REF] Rivera-Cárcamo | Control of the Single Atom/Nanoparticle Ratio in Pd/C Catalysts to Optimize the Cooperative Hydrogenation of Alkenes[END_REF] Such clusters (according to their size) are representative of highly dispersed palladium catalysts (see Erreur ! Source du renvoi introuvable.b and Figure S6 for Pd/CNT catalyst. For this, an isolated Pd 13 cluster according to geometric shell model was first optimized, obtaining a classical symmetric icosahedral (I h ) structure. It was then deposited on the surface containing the Pd SA , optimizing the entire supported catalyst, the final structure of which can be found in Figure 6e. As it can be seen, the anchoring of the cluster took place in the carbonyl-decorated multiple vacancy.

In addition, twelve of the thirteen Pd atoms are accessible, so the theoretical metal dispersion is close to 100%. One of the main reasons is that the cluster is located at a reasonable distance from the Pd SA , minimizing the interaction, and thus preventing the agglomeration of the latter (the shortest distance between a Pd atom on the surface of the cluster and the Pd SA is greater than 5 Å).

The Pd SA maintained the same structural properties (average Pd-C bond length of 1.96 Å and elevation height of 1.02 Å) as before anchoring the cluster, which endorses the weak interaction between both species. Moreover, a very strong interaction (adsorption energy of 8.8 eV) is denoted between the Pd 13 cluster and the support, in agreement with the values reported in the literature for the same cluster size. [START_REF] Yang | A Fundamental Analysis of Enhanced Cross-Coupling Catalytic Activity for Palladium Clusters on Graphene Supports[END_REF] The size of the vacancy is of paramount importance to accommodate and tightly binding the metal clusters, as large vacancy defect sites allow the formation of multiple Pd-C bonds that increase the binding. This inevitably leads to a final cluster slightly distorted compared to the initial icosahedral structure, where the Pd-C bonds are in the range between 2.1 and 2.7 Å, and those of Pd-O between 2.0 and 2.3 Å. For all these reasons, the DV motif turns out to be another way to prevent the progressive leaching of the cluster during catalysis.

Through the structure-property analysis of the supported cluster, few meaningful structural parameters are investigated. For example, making use of the nearest neighbor search, the average first nearest-neighbor (NN) Pd-Pd bond length (R) is of 2.62 Å, significantly shorter than the distance obtained by EXAFS for Pd/C catalysts with larger NP size (around 2.74 Å). [START_REF] Zlotea | Size-Dependent Hydrogen Sorption in Ultrasmall Pd Clusters Embedded in a Mesoporous Carbon Template[END_REF][START_REF] Bugaev | Hydride Phase Formation in Carbon Supported Palladium Hydride Nanoparticles by in Situ EXAFS and XRD[END_REF] The average NN Pd-Pd coordination number (CN) is of 4.92, although with clearly different values in the case of the Pd atom located in the core of the cluster (CN = 10) and those located on the surface (CN = 4.50). Additionally, the mean interatomic distance (MIAD) was calculated, which reflects the effective radius of a cluster and allows, for example, an analysis of morphological change. [START_REF] Wang | Shear Instabilities in Metallic Nanoparticles: Hydrogen-Stabilized Structure of Pt37 on Carbon[END_REF] In this case, MIAD had a value of 3.77 Å. These numbers are in excellent agreement with those of isolated Pd 13 reported in the literature, [START_REF] Wang | Density Functional Study of Structural Trends for Late-Transition-Metal 13-Atom Clusters[END_REF] with R = 2.64 Å, CN = 6.46, and MIAD = 3.62 Å, highlighting that the fact of being supported caused no excessive structural change, only a relative decrease in the CN mainly due to a partial cluster opening on the face closest to the support. To provide a more precise information on the structure of the cluster during the adsorption process, the pair distribution functions of both the isolated and anchored species were analyzed (see Pd 13 black and Pd 13 @G ox blue curves of Figure S7, respectively). While the peaks in both situations appeared at approximately the same Pd-Pd bond lengths, there is an increase in metallic interactions at higher distances in the case of the anchored Pd 13 , which would explain the differences observed for the CN.

In Figure 6f has scarcely influenced the local electronic structure of the Pd SA . Regarding the pCOHP plot, the presence of bonding states between -8 and -4 eV and antibonding states between -4 and -1 eV with respect to Pd SA is remarkable, related to the deposition of the Pd 13 cluster. The highest occupied states are essentially established on metal atoms (Pd-Pd interactions), exhibiting an antibonding character (negative values of -pCOHP). The latter specifically bears similarity to the electronic structure of small organometallic clusters. [START_REF] Del Rosal | Ligand Effect on the NMR, Vibrational and Structural Properties of Tetra-and Hexanuclear Ruthenium Hydrido Clusters: A Theoretical Investigation[END_REF] The bond strength index IpCOHP is found to be 0.9 eV per Pd-Pd interaction, somewhat higher than the value of 0.7 eV obtained for the isolated nonanchored cluster, which suggests a metal-metal strengthening by the Pd 13 deposition on G ox .

An important change can also be deduced in the electronic distribution of the cluster, in line with what is observed when the Pd SA is anchored. The mean valence electron occupation of each of the palladium atoms is 4d 9.2 5s 0.6 , so taking into account the 13 atoms there is a total difference of three electrons less with respect to the expected value. This means that apparently the cluster is formally [Pd 13 ] 3+ . However, among all the Pd atoms of the cluster, the one located at the base had a lower electronic occupation than the rest, reminiscent of a Pd 2+ in a square-planar environment according to the optimized geometry, so formally the cluster is more close to

[(Pd 1 ) 2+ (Pd 12 ) + ].
G ox is then acting as an electron acceptor when it is in direct contact with a metal such as palladium, regardless of whether it is presented as SA or ultra-small NP, which is in agreement, according to the work function of both the metal and carbon or oxygen from the support. This outcome can also be confirmed together with the charge density map in Figure 6g.

Clearly, Pd SA continued to have a very positive charge density, while Pd 13 offered a wide range of densities.

Finally, the IR spectrum with both forms of metallic species is shown in Figure 6h. The clearest evidence that would reveal the presence of the Pd 13 cluster corresponds to ν(Pd-Pd), but unfortunately, they appear at frequencies below 300 cm -1 , characteristic of the fingerprint region, often the most complex and confusing region to interpret. However, other evidence in more favorable regions is observed, as is the case of the ν(C=O) of carbonyl groups sharing electronic density with Pd 13 atoms. For example, C=O stretching vibration of carboxylic acid groups in G ox usually appear in a range of about 1710 and 1740 cm -1 , but when they are coordinating to an atom of the metal cluster a significant red-shift is denoted, in this case between 1560 and 1630 cm -1 .

Similar behavior with the stretching vibration of ketone (and quinone) groups is appreciated, since in their free coordination state they appear between 1620 and 1680 cm -1 (1640 and 1700 cm -1 for quinone), and when coordinated to a Pd atom the peak is red-shifted to values of about 1230 and 1410 cm -1 (approx. 1320 cm -1 for quinone). Therefore, while direct measurements of metal-metal bonds are difficult to characterize by IR spectroscopy, indirect measurements such as surface adsorbates coordinated to metal appear as a plausible solution.

Cluster hydrogenation. Supported Pd catalysts are attractive for hydrogenation reactions, since palladium is considered the most selective among Group 10 metals. High dispersion is a prime asset of these catalysts, where the small metal NP have the ability to be hydrogenated. This also results in an important class of materials with significant potential for energy storage and energy-related applications. [START_REF] Zhou | Hydrogen Storage Performance in Pd/Graphene Nanocomposites[END_REF] By means of a systematic analysis, the hydrogenation of the supported Pd 13 cluster is evaluated on the G ox . Since H 2 dissociates easily at the cluster surface, it is basically covered with surface hydrides. According to what has been previously reported for hydride Ru NP , 63,64 a hydrogen to metal ratio (H/M) of approximately 1.5/1 has been used to ensure complete saturation, so the final stoichiometry turned out to be Pd 13 H 20 , as illustrated in Figure 7a.

It is important to highlight the outcropping towards the surface that metal cluster experienced with respect to the situation prior to hydrogenation, where the binding to the support was very noticeable. Even so, this prominent elevation has already been observed for γ-alumina-supported Pt nanoparticles. [START_REF] Mager-Maury | H2-Induced Reconstruction of Supported Pt Clusters: Metal-Support Interaction versus Surface Hydride[END_REF] The metallic cluster went from being anchored by means of a square planar molecular geometry formed by four carbonyl groups, in addition to numerous interactions between the rest of Pd atoms and graphene carbons, to being anchored by only two CO groups and three H atoms. This could cause potential stability problems under H 2 environments as is often observed for Pd NP , 66 resulting in sintering of the metal, although in general it occurs at relatively elevated temperatures. The Pd SA is not affected by the hydrogenation process, with an elevation height of 1.00 Å in line with the previous values. Regarding the nature of hydrogen atoms, they sit in high coordination sites, showing preference for hollow adsorption sites (coordinated to three Pd atoms), and to a lesser extent for the bridge ones (coordinated to two Pd atoms), revealing no hydrogen adsorption at top sites (coordinated to a single Pd atom), as occurs on extended Pd surfaces. [START_REF] Dong | Hydrogen Adsorption on Palladium: A Comparative Theoretical Study of Different Surfaces[END_REF] Since the twenty H atoms can be interpreted as ten fully dissociated H 2 molecules, the hydrogenation energy of Pd 13 can be calculated as the difference between Pd 13 H 20 and the energy of Pd 13 plus the energy of ten H 2 , obtaining a favorable value of -4.6 eV. On the other hand, the adsorption energy directly from the hydrogenated cluster can also be calculated, and it is found to be -4.2 eV, so from an anchoring point of view, the process appears to be less thermodynamically favorable compared to the bare cluster. A structural analysis of the cluster using the NN algorithm revealed an average Pd-Pd bond length and an effective radius of 2.74 and 3.95 Å; around 0.1 and 0.2 Å higher, respectively compared to the non-hydrogenated situation, showing a generalized swelling. The average Pd-H bond length is found to be 1.75 Å, and the mean Pd-Pd and Pd-H coordination number is 4.92 and 2.75 respectively, a value in the case of the former CN identical to the bare cluster. Therefore, there has not been a loss of local coordination of Pd atoms by the dissociative adsorption of the H 2 molecules. Considering the pair distribution function of the hydrogenated cluster (see Pd 13 H 20 @G ox red curve in Figure S7), the appearance of a small intensity peak around 3.65 Å (almost 1 Å larger than the main peak) would also explain the slight increase in observed size. Furthermore, the hydrogenation of the anchored Pd SA is also evaluated in the G ox . According to the calculations, H 2 can be easily adsorbed on the metal center with an adsorption energy around 0.4 eV. However, both homolytic and heterolytic H 2 activation is clearly endothermic (reaction energies greater than 0.4 eV relative to bare G ox ), so it will not take place. It appears, therefore, that the only plausible way to have an H atom on Pd SA forming the Pd SA -H species is through an H-spillover mechanism.

The corresponding pDOS and pCOHP plots are represented in Figure 7b. The presence of hydrogen atoms in the Pd cluster caused the total pDOS to be dominated with less intensity by the metal in the states around the Fermi level compared to the situation prior to hydrogenation. The dband center for both the Pd NP and the Pd SA is calculated to be 2.6 and 4.5 eV, slightly higher for the case of the former, and exactly the same for the case of the latter compared to the prehydrogenated situation. This stabilization in the metallic cluster was also computationally observed for ruthenium NP. [START_REF] Cusinato | Theoretical Characterization of the Surface Composition of Ruthenium Nanoparticles in Equilibrium with Syngas[END_REF] It is worth highlighting the appearance of states around -9 eV corresponding to the new H atoms, clearly differentiated from the ones belonging to the carboxylic acid groups, located around -15 eV. Below the Fermi level, a subtly more accentuated antibonding character than in the case of the bare Pd 13 cluster is revealed from the pCOHP plot for the highest Pd states, which further weakens the metal-metal bond, going from an IpCOHP index of 0.9 to only 0.4 eV per intermetallic interaction. However, a value of 2 eV is calculated for the newly exposed metal-H interaction. These results provided evidence on the existing polarization in the new hydrogenated clusters, as pointed out by the pMPA maps in Figure 7c. Comparing both bare and hydrogenated metallic NP, a clear tendency to oxidation of the palladium atoms is observed, turning an increasingly intense blue color in the case of the surface atoms and a whitish color in the case of the core atom (it should be noted that the core atom had a slightly negative atomic charge). Instead, hydrogen atoms could be considered as hydrides, as indicated by their reddish colorations. A more in-depth analysis on the electronic structure of the cluster was carried out to elucidate the charge balance between metals and adsorbates. Because the metal is in a low oxidation state (it is originally a d 10 , so rich in electrons) and the hydrogen adsorbate possesses low-lying empty orbitals, then a metal-to-adsorbate charge transfer transition may occur.

Specifically, up to four and one electron were transferred to the hydrogen adsorbates and also to G ox , respectively, so finally the charge balance elucidated a [Pd 13 H 20 ] + cluster.

The simulated IR spectrum is depicted in Figure 7d. The most outstanding difference compared to the spectrum of the bare cluster appeared around 1171 and 1340 cm -1 , corresponding to the ν(Pd-H) stretching vibrations. The normal modes of the metal-hydrogen interactions were observed in a wide range of frequencies due to the heterogeneity of palladium and hydrogen atoms forming the metallic nanostructure, ranging from 1100 to 1450 cm -1 . Therefore, the IR bands at 1243 and 1420 cm -1 , related to the ν(C-O) of lactone and ether functional groups, and ν(C=C) of aromatic carbons, respectively, could have experienced a relative increase in intensity. A decay of the infrared intensity appearing at 995 cm -1 , assigned to the ν(Pd-C) stretching mode, can also be detected mainly as a consequence of the elevation experienced by the Pd 13 H 20 cluster, and therefore the gradual disappearance of metal-carbon interactions.

Thermal stability. The thermal stability of the supported metal catalyst is also evaluated by means of a MD simulation at room temperature, in a procedure analogous to that carried out for the G ox support (the evolution of the energy of the system over the simulation time can be found in Figure S8). As expected, the data analysis of the simulation trajectory mostly revealed changes in the cluster region, which underwent a gradual deformation, inevitably leading to a nonreversible process. Although it maintained the distorted icosahedron shape as shown in Figure 7e, it is interesting to note the creation of local octahedra that were previously unappreciated. The most striking feature of this Pd atom rearrangement is the hydrogen recombination on the surface of the cluster, probably due to the high coverage of H atoms originally introduced. The appearance of dihydrogen took place from the first picoseconds of simulation, after which up to three new molecules appeared. Once formed, the H 2 molecules remained coordinated at the Pd 13 surface, suggesting a higher H 2 adsorption energy than the thermal energy to desorb. This resulted in a stoichiometry of Pd 13 H 14 (H 2 ) 3 , decreasing the effective H/M ratio to practically unity. Dihydrogen complexes M(η 2 -H-H) are notable in palladium chemistry, definitely acting as σ complexes in reaction intermediates. [START_REF] Andrews | Matrix Infrared Spectra and Density Functional Calculations of Transition Metal Hydrides and Dihydrogen Complexes[END_REF] Taking a more detailed analysis of the structure-property of the renewed catalyst, the energy difference between the supported structures Pd 13 H 20 and Pd 13 H 14 +3H 2 is calculated to be - S7) revealed a significant reduction in the number of closest metallic interactions (peak around 2.65 Å), and the appearance of a small peak at Pd-Pd bond lengths of approximately 6.5 Å, which would explain the decrease in intermetallic coordination observed. For all this, the stability achieved by the combined H/H 2 species is largely driven by the rearrangement of the hydrogen atoms, contributing to a lesser tendency towards sintering processes, operating at elevated temperatures in hydrogen-rich atmospheres.

Since the composition of these species is essentially the same, subtle changes in the electronic structure are appreciated. As the pDOS and pCOHP plots of Figure 7f Finally, in Figure 7h the simulated IR spectrum of Pd 13 H 14 +3H 2 is represented. Due to the recombination of some H atoms into molecular hydrogen, a small peak at 3080 cm -1 corresponding to ν(H-H) is originated. It is experimentally observed that the H-H stretching mode for the sidebonded Pd(H 2 ) complex occurs at 2971 cm -1 with a very low relative integrated intensity, [START_REF] Andrews | Observed and Calculated Infrared Spectra of Pd(H2)1,2,3 Complexes and Palladium Hydrides in Solid Argon and Neon[END_REF] consistent with what has been computationally observed so far, the small differences probably coming from the size of the metallic entity. As a consequence of the new metal-dihydrogen interaction, the corresponding ν(Pd-H 2 ) vibrational mode is perceived at frequencies significantly higher than those of Pd-H δ-band. Specifically, it can be found from 1450 to 1650 cm -1 , being extended the entire band in values between 1100 and 1650 cm -1 . Therefore, the appearance of possible changes in the nature of the supported catalyst could be diagnosed via IR spectroscopy by comparing with the reported computational data. with the metal, in addition to a significant reduction in the electron density of the supported metal species. Finally, the nanocluster has been fully hydrogenated, as a relevant state in hydrogenation processes, being less exposed to the surface due to the absorbed hydrogen atoms. However, a second thermal equilibration of the system has shown the recombination of up to three H 2 molecules, specifying that the H/M ratio for the studied system is approximately unity. § ASSOCIATED CONTENT

Supporting Information

Highly Active, Selective, and

Scheme 1 .

 1 Scheme 1. Structure and chemical composition of the nanomaterials system.

Figure 1 .

 1 Figure 1. Pristine graphene (a) and doubly V2(555-777) defect (b) created in graphene.

Figure 2 .

 2 Figure 2. Schematic representation of the different oxygen functional groups (OFG) experimentally observed and their corresponding surface density.

Figure 3 .

 3 Figure 3. AIMD of the Gox support, with the corresponding top and side views of initial (a) and final (b) structures shown as insets.

Figure 4 .

 4 Figure 4. Characterization of graphene and Gox. (a and b) Projected DOS and COHP profiles for the nearest-neighbor interactions. The position of the Fermi level (brown horizontal line) is also reported. pCOHP plots are calculated for nearest neighbors only. (c and d) Top and side views of projected MPA atomic charges. (e and f) Simulated IR and Raman spectra.

Figure 5 .

 5 Figure 5. STEM-HAADF micrographs showing mixtures of PdNP and PdSA in: (a) 5% Pd/C; (b) 1.2% Pd/CNT; and (c) 1.5% Pd/FLG catalysts. (d) Mean PdNP size and % of PdSA in the Pd/C, Pd/CNT, and Pd/FLG catalysts.

Figure 6 .

 6 Figure 6. Characterization of anchored PdSA and PdSA+PdNP catalysts. (a and e) Top and side views of atomistic structures. (b and f) Projected DOS and COHP profiles. (c and g) Top and side views of projected MPA atomic charges. (e and h) Simulated IR spectra. Conventions are the same as in Figure 4. The d-band center for SA (dashed red line) and NP (dotted black line) are also reported.

  pDOS(ε) and pCOHP(ε) of Pd SA are plotted in Figure 6b. Around the Fermi level, the total pDOS (in black) is still influenced by the G ox (in blue), since the proportion of dispersed metal with respect to the support is very low. The Pd d-band center, ε d , calculated as the normalized, energy-weighted integral of the DOS, and projected onto all d AOs, is calculated to be 4.6 eV (in red on the pDOS plot). Whereas below the Fermi energy pCOHP (C-C and O-H, ε) are always positive, that of C-O showed a substantial back-donation in line with G ox , so its origin must come from the π* MO. In addition, the highest metal states have an antibonding Pd-C character. It was found that IpCOHP has exactly the same values per C-C, C-O, and O-H interactions as in the case of G ox , while a value of 3 eV is taken for the newly formed Pd-C interactions. Furthermore, a change in the partial atomic charge of the Pd atom is observed by analyzing the electronic distribution. The electronic configuration of the ground state of gaseous neutral palladium is [Kr]

  , the pDOS(ε) and pCOHP(ε) of the Pd SA and Pd NP species are shown. Contrary to what happens in Pd SA , a greater contribution by metal atoms in the electronic states around the Fermi energy is observed. In this case, the metallic character of the cluster can be clearly appreciated, mainly by the Pd d component (in yellow) towards values below the Fermi level and by the Pd s component (in green) at values above the Fermi energy. The Pd d-band center for both the cluster (dashed red line) and the single atom (dotted black line) anchored on the support is calculated to be 2.0 and 4.6 eV, respectively. In the case of the former, the d-band center stabilizes 0.3 eV with respect to the bare cluster by means of bonding interactions with surface C atoms, widening the d band. On the other hand, with an identical value of ε d for Pd SA , the presence of Pd 13

Figure 7 .

 7 Figure 7. Characterization of hydrogenated Pd13H20 and Pd13H14+3H2 catalysts. (a and e) Top and side views of atomistic structures. (b and f) Projected DOS and COHP profiles. (c and g) Top and side views of projected MPA atomic charges. (e and h) Simulated IR spectra. Conventions are the same as in Figure 4 and Figure 6.

  pointed out, the emergence of states between -11 and -10 eV (always positive below the Fermi energy) is the most relevant difference compared to Pd 13 H 20 , belonging to the H 2 coordination (light blue curve in the pCOHP plot). These results give an interesting indication that co-adsorbed species capable of mobilizing d AOs from the Pd cluster surface could favor sideway H₂ coordination. Regarding the pMPA charges represented in Figure7g, a general withdrawal of the positive electron density of the Pd cluster atoms is appreciated, promoted by the migration of hydrogen atoms (i.e., by the negative charge density migration due to partial hydride character) during recombination into molecular hydrogen. In turn, the latter all elucidated a symmetry in the electric charge of the H-H bond, expressing a dominant homolytic component.

  functionalized graphene (carbon support) as well as supported Pd SA and Pd 13 nanocluster (supported metal catalysts) have been computationally modelled by means of periodic DFT simulations. Structure-property characteristics, as well as electronic properties and vibrational features have been analyzed in detail. For the functionalized carbon support, the inclusion of point defects and experimentally probed abundant OFG have been taken into account, demonstrating high thermal stability at room temperature. The characterization has confirmed the structural motifs introduced in the model, distinguishing the organic functionalities and the vacancy defect by IR and Raman spectroscopy, respectively; and showing the nature of the new interactions according to their electronic structure, which revealed a material of metallic character. On the other hand, the simultaneous inclusion of both classes of metallic species has been carried out in clearly differentiated regions of the support surface and directly on the defects to prevent processes such as agglomeration. It is important to highlight the red shift of the IR bands governed by interactions

Table 1 .

 1 Obtained IR bands (in cm -1 ) from the simulated spectrum of G ox , a assignment of the vibration mode of the different functional groups, and experimental IR frequencies for the CNT ox support. b

	IR band	Vibrational mode	Functional group	CNT ox
	3648	ν(O-H)	Carboxyl	3260
		ν(C=O)	Lactone	-
	1748
		ν(C=O)	Carboxyl	1724
	1540	ν(C=O)	Carbonyl	1680
	1477	ν(C=C)	Aromatic	1568
		ν(C-O) [acyl]	Lactone
	1260
		ν(C-O)	Ether
		ν(C-O)	Carboxyl	1180
	1123
		ν(C-O)	Ether
	1009	ν(C-O) [alkoxy]

  1.3 eV, so clearly the system obtained after thermal equilibration was more stable. To ensure if the energy variation came mainly from the structural reorganization of the cluster framework or from the new arrangement of the H atoms, again parameters such as average bond length, mean CN, and MIAD were calculated. For the Pd 13 H 14 +3H 2 species, the intermetallic and metal-H distances, as well as the MIAD, were 2.66, 1.73, and 3.92 Å, respectively; all values very close to those calculated for the fully hydrogenated species. In fact, the Pd-Pd distance returned to the same value as the isolated gas-phase cluster (around 2.64 Å). Coordination numbers were the ones that caused the most difference: around 4.31 (CN Pd-Pd ) and 2.05 (CN Pd-H ) compared to 4.92 and 2.75 in the case of the Pd 13 H 20 stoichiometry. This decrease was probably caused by the local opening of the bottom of the metallic structure, causing a Pd-C interaction to compensate for the loss of coordination. The pair distribution function of the thermalized cluster (see Pd 13 H 14 +3H 2 @G ox green curve in Figure

  Stable Pd Single-Atom Catalyst Anchored on N-Doped Hollow Carbon Sphere for Electrochemical H2O2 Synthesis Under Acidic Conditions. J. Machado, B. F. Nanostructured Carbon Materials for Catalysis; Catalysis Series; The Royal Society of Chemistry, 2015. (10) Machado, B. F.; Serp, P. Graphene-Based Materials for Catalysis. Catal. Sci. Technol. 2012, 2 (1), 54-75. (11) Xi, J.; Jung, H. S.; Xu, Y.; Xiao, F.; Bae, J. W.; Wang, S. Synthesis Strategies, Catalytic Applications, and Performance Regulation of Single-Atom Catalysts. Adv. Funct. Mater. Angew. Chemie Int. Ed. 2008, 47 (20), 3696-3717. (14) Lui, C. H.; Liu, L.; Mak, K. F.; Flynn, G. W.; Heinz, T. F. Ultraflat Graphene. Nature 2009,

		Catal. 2021, 393, 313-323.
	(9)	Serp, P.; 2021, 31 (12), 2008318.
		462 (7271), 339-341.
	(15) Chen, D.; Feng, H.; Li, J. Graphene Oxide: Preparation, Functionalization, and
		Electrochemical Applications. Chem. Rev. 2012, 112 (11), 6027-6053.

(12) Boukhvalov, D. W.; Katsnelson, M. I. Chemical Functionalization of Graphene with Defects. Nano Lett. 2008, 8 (12), 4373-4379. (13) Liang, C.; Li, Z.; Dai, S. Mesoporous Carbon Materials: Synthesis and Modification. (16) Szabó, T.; Berkesi, O.; Forgó, P.; Josepovits, K.; Sanakis, Y.; Petridis, D.; Dékány, I. Evolution of Surface Functional Groups in a Series of Progressively Oxidized Graphite

Notes

The


Mirella Aoun
email: mirella.aoun@univ-rouen.fr

Clément Berger
email: clement.berger@ens-lyon.fr

Jean Cazalis
email: cazalis@ceremade.dauphine.fr

Thierry Gonon
email: thierry.gonon@ec-lyon.fr

David Lassounon
email: enagnon-david.lassounon@insa-rennes.fr

Characterization of the fluctuations of an ultrasonic wave passing through a complex environment in order to simplify the modeling

Keywords: ultrasound, Gaussian processes, kriging, speckles, noise, wave fluctuations

à 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.

Introduction

This document summarizes a work that has been realized in four days for the CEA (Commissariat à l'énergie atomique et aux énergies alternatives) during the SEME (Semaine d'Etudes Mathématiques-Entreprises) organized by IP Paris (Institut Polytechnique de Paris) and AMIES (Agence pour les Mathématiques en Interaction avec l'Entreprise et la Société) in February 2022. We would like to thank our CEA supervisors, Vincent Dorval and Edouard Demaldent, for their support during this week. We also thank Samuel Amstutz as well as the whole staff of the SEME for providing us with this opportunity.

Ultrasonic testing

Many efforts are being made by the electric power industry in order to ensure security and to extend the life of electric power stations such as nuclear plants or hydroelectric dams. For this purpose, it is crucial to improve the detection of defects (inclusions, blisters, defects in sticking, etc.) or damages in these structures. However tests must not affect the integrity of the materials and thus require the use of so-called nondestructive testing methods.

Among these methods, ultrasonic testing is frequently used because it has many advantages: ease of implementation, possibility of working on a single side of the part to be inspected (no need for access to the second side), and the ability to pass through large thicknesses of material depending on the working frequency. Moreover, the existence of relationships between the propagation of ultrasound and the characteristics of the material allows its characterization. Therefore, the use of ultrasound offers the possibility, without any deterioration, on the one hand of characterizing materials in order to know their elastic properties, and on the other hand of inspecting parts to check their integrity. A challenge of this method comes from the heterogeneity of the materials which are being put to the test. Typically, concrete walls contain multiple aggregates at random places (see figure 1) which disturb the propagation of the ultrasonic waves. Numerical simulations of wave propagation through complex materials by finite-element methods are introduced by the ATHENA code developed by EDF and INRIA. ATHENA simulates the propagation of ultrasonic waves in a plane of symmetry of heterogeneous complex anisotropic materials with applications on stainless welds [START_REF] Chassignole | Influence de la structure métallurgique des soudures en acier inoxydable austénitique sur le contrôle non destructif par ultrasons[END_REF][START_REF] Chassignole | Ultrasonic and structural characterization of anisotropic austenitic stainless steel welds: Towards a higher reliability in ultrasonic non-destructive testing[END_REF][START_REF] Ploix | Étude de l'atténuation des ondes ultrasonores. Application au contrôle non destructif des soudures en acier inoxydable austénique[END_REF]. The CEA uses a more recent code developed in 2019 by Imperiale and Demaldent [START_REF] Imperiale | A macro-element strategy based upon spectral finite elements and mortar elements for transient wave propagation modeling. application to ultrasonic testing of laminate composite materials[END_REF]. However these simulations require the exact positions of the heterogeneities, which are in general not available. In addition, these methods are computationally expensive which prevents their use in a large scale. A way out consists in focusing on the averaged ultrasonic field, which makes it possible to reduce the problem to the easier case of a homogeneous material. Nevertheless, because this approach does not take the noise induced by the heterogeneities into account, it is not precise enough for industrial purposes. One could improve the procedure by modeling the fluctuations around the average wave.

Our framework

In this work, we propose to randomly simulate wave profiles which are similar to simulations given by a finite-element based solver.

The simulations are performed on a 2D rectangular grid. We will denote space coordinates by x and z and the time by t. We follow a planar wave u(x, z, t) ∈ R 2 evolving along the z axis, parallel to the x axis. This work is dedicated to the simulation of the deformation along the z axis, denoted by u z . There are 100 points in z and 200 points in x. Because the finite-element solver has an adaptative time discretization, the number of time steps depends on the simulation. There are roughly between 100 and 200 time steps in the data.

As mentioned earlier, the wave averaged along the x axis can be easily simulated. Thus we will focus only on the remaining fluctuations caused by the aggregates. Finally, at each time t, the space can be roughly divided into three different areas. Starting from the top, we have a zone with u z = 0 since the wave has not reached this zone yet. Getting down we meet the front wave, corresponding to the wave we would have if the material was homogeneous. Close to this area, the fluctuations are composed of many small disorganized patches. Finally, under this area lies the so-called coda, which corresponds to the residues of the wave. This setup is described in Figure 2. In this work, it is assumed that the coda is the main contributor to the noise measured during the testing procedures. As a consequence, Figure 2: Scheme of the setup. The two upper images show the decomposition of the wave between its average and the noise resulting from the aggregates. The decomposition of the space between the three areas is also presented.

we mainly focus on its modeling.

This paper is presented as follows: in Section 2, we present a first method based on Gaussian processes. Another method based on a dynamic update is presented in Section 3. Finally we propose different indicators in Section 4 in order to get quantitative evaluations of the performances of the methods proposed in this report.

Gaussian processes based approach 2.1 Initial model

The first method is inspired from the Gaussian process regression formalism [START_REF] Santner | The design and analysis of computer experiments[END_REF][START_REF] Williams | Gaussian processes for machine learning[END_REF]. The z component of the residual fluctuations u z is supposed to be the realization of a Gaussian process U z (x, z, t) depending on the space and time variables:

U z (x, z, t) = Y (x, z, t) • 1 ]0,+∞[ (z -vt) .
The coda is modeled by the process Y . We have denoted by v the group velocity in the direction z of the ultrasonic wave and we compute it using the simulations provided by the CEA. The function 1 ]0,+∞[ (z -vt) is equal to 1 if z -vt ∈]0, +∞[ and 0 otherwise. In other words, we set the noise to zero in the area {z > vt} that the wave front has not reached yet.

The process U z is simulated by generating a path of Y and multiplying it by the Heavyside function of the input space (x, z, t). The law of Y is chosen equal to:

1 ]0,+∞[ (z -vt). The path of Y is generated on a grid G = {(x i , z i , t i ), 1 ≤ i ≤ N }
Y ∼ GP 0, k((x, z, t), (x ′ , z ′ , t ′ )) ,
where we denote by k the covariance kernel. We choose k as a Matérn kernel with parameter 5 2 (a usual choice in the Gaussian regression community):

k((x, z, t), (x ′ , z ′ , t ′ )) = σ 2 1 + √ 5|x-x ′ | θx + 5(x-x ′ ) 2 3θx exp - √ 5|x-x ′ | θx • 1 + √ 5|z-z ′ | θz + 5(z-z ′ ) 2 3θz exp - √ 5|z-z ′ | θz • 1 + √ 5|t-t ′ | θt + 5(t-t ′ ) 2 3θt exp - √ 5|t-t ′ | θt .
The generation of the Gaussian vector Y (G) of law:

Y (G) ∼ N (0, k(G, G)) ,
is done by simulating a standard Gaussian vector ϵ ∼ N (0, I) and multiplying it by the square-root of the covariance matrix k(G, G) (see Algorithm 1):

Y (G) = k(G, G) • ϵ.
Notice that the computation of k(G, G) can be done once and for all before performing multiple simulations. An illustration of a simulation is given in Figure 3, for a 25 × 25 × 25 grid in [0, 1] 3 where the parameters are arbitrarily set to

σ 2 = 1, θ x = θ z = θ t = 0.1, v = 1.

Refinements

• To smooth the wave front, the Heavyside function used for the cutoff can be replaced by a function of the form: 1 2 -1 π arctan(α(z -vt)) with α a parameter to adjust.

• The model parameters may be adjusted numerically, instead of manually, by maximiz-Algorithm 1 Simulation based on Gaussian processes.

1: Set G of size n, σ 2 , θ x , θ z , θ t . 2: Generate ϵ ∼ N (0, I n ). 3: Compute k(G, G) = (k((x i , z i , t i ), (x j , z j , t j ))) i,j . 4: Decompose: k(G, G) = P        λ 1 0 . . . 0 0 . . . . . . . . . . . . . . . . . . 0 0 . . . 0 λ n        P T . (2.1)
▷ Eigenvalue decomposition 5: Take the square root:

k(G, G) = P        √ λ 1 0 . . . 0 0 . . . . . . . . . . . . . . . . . . 0 0 . . . 0 √ λ n        P T . (2.2) 6: Compute: Y (G) = k(G, G)ϵ. ( 2
H(i) = 1 13:
end if 14: end for 15: Final simulation:

U z (G) = Y (G) T H. (2.4)
ing the loglikelihood for the observations of one simulation:

l(σ 2 , θ x , θ z , θ t ) = -log |k(G, G)| - u z T k(G, G) -1 u z 2 .
• To take into account the fact that the fluctuations may vanish with time and far from the wave front, we can multiply the simulation by exp(β(z -vt)) with β a parameter to adjust.

Drawbacks

While the results are promising, the size of G makes the computations of k(G, G) (for the simulations) and k(G, G) -1 (for fitting the parameters) computationally expensive. Truncations of the covariance matrix have been attempted in order to obtain sparser covariance matrices. However, even with these adjustments, the calculation time of k(G, G) for a 25 × 25 × 25 grid takes hours. In conclusion, it is not realistic to use this method for bigger grids and multiple different parameters.

Dynamic simulation with smoothing

In this section we build upon the work of [START_REF] Federico | Simulation of dynamic speckle sequences and its application to the analysis of transient processes[END_REF] to simulate the noise. In [START_REF] Federico | Simulation of dynamic speckle sequences and its application to the analysis of transient processes[END_REF], the authors present a numerical simulation to generate a temporal sequence of dynamic speckle patterns which is based on a model on a superposition of waves from discrete scattering centres.

A first implementation

Algorithm 2 Dynamic simulation of the noise with space smoothing. 1: Set ϕ(., ., t = 0), c > 0 decreasing, T max > 0. Compute:

ϕ(x, z, t + ∆t) = ϕ(x, z, t) + N (0, 1) log c(t + ∆t) -log c(t)). (3.1)
▷ Computation of the new phase 5:

Compute:

ũz (x, z, t) = R F -1 (H ω 0 (F(e iϕ(.,.,t) )))(x, z, t) . (3.2)
▷ Application of low pass filter

6:
end for 7: end for 8: Compute:

u z (x, z, t) = π 2 -arctan(α(z -vt)) π ũz (x, z, t). (3.3)
▷ Application of arctan to recreate the no-wave area Let's briefly describe the algorithm proposed in [START_REF] Federico | Simulation of dynamic speckle sequences and its application to the analysis of transient processes[END_REF]. The idea is to work on the phase ϕ of the wave, with a time iterative update. We start with ϕ(t = 0) given by sampling uniform (over [-π, π]) random variables on the space grid. Then at each time, we update the phase according to the formula (3.1) in which the phase at each time consists in the previous phase shifted by a centered random variable. The quantity log c(t + ∆t) -log c(t) allows us to control the variance of the update, where the time function c represents the correlation between the image at time t and the image at time 0. It is according to the CEA's simulations using the formula:

c(t) = E[I(0)I(t)] -E[I(t)]E[I(0)] V[I(0)]V[I(t)] , ( 3.4)
where E (resp. V) denotes the mean operator (resp. the variance operator).

The resulting image is a white noise evolving with time, with a degree of time correlation. To add spatial correlation, a low pass filter is applied on each image, see (3.2). More specifically, we perform a discrete Fourier transform on the image at each time step then we apply a cutoff function H ω 0 on both space variables: given a threshold ω 0 in the frequency domain, for each frequency ω, the corresponding component is set to 0 if ω > ω 0 and is let unchanged otherwise. We then perform an inverse Fourier transform to get our result.

To finish we have to apply, as in Section 2, a cutoff on each image in order to set u z to 0 in the area that the wave has not reached yet. As suggested earlier, we use the formula (3.3) based on the arctan function rather than a Heavyside in order to obtain a more regular front. The complete algorithm is summed up in Algorithm 2.

An example of simulation using this algorithm is shown in Figure 4. The resulting image is close to the one shown in Figure 2. Also note that this time it is possible to work on the full space-time grid, within reasonable time (approximately 0.3s).

However when displaying the evolution of the noise, we observe that in our simulations the patches are being irregularly deformed and have no visible motion whereas in the real simulations, the noise behaves more like smudges coming from the front wave and going down.

Modifications to enhance temporal behavior

To mimic this behavior we made some adjustments to the original algorithm.

We start by applying a low pass not only in space, but also in time. This means that we perform a three-dimensional discrete Fourier transform on the whole sequence of images before applying a filter B (for consistency we changed the notation as it is now a threedimensional function), as shown in (3.6). This smooths out the time evolution of the noise, making its appearance closer to the ones in the CEA simulations. Since this filter already manages the time correlation, we also remove the function c and simply use a fixed variance σ 2 for the Gaussian random variables used in the phase update.

Another modification has been implemented in order to recreate the smudges going down, as described earlier. In the phase update, we take a combination between the previous update (3.1) and the value of the phase two space steps above. This allows the information coming from above to partially propagate downwards. This is done by a convex combination of both Algorithm 3 Dynamic simulation of the noise with space-time smoothing and added space correlations.

1: Set ϕ(., ., t = 0), corr ∈ [0, 1], σ corr ≥ 0, σ ≥ 0, T max > 0.

2: for 0 ≤ t ≤ T max -∆t do 3:

for all x, z do 4:

Sample: corr ef f = N (corr, σ 2 corr )

5:

Compute:

ϕ(x, z, t + ∆t) = (1 -corr ef f ) ϕ(x, z, t) + N (0, σ 2 ) + corr ef f × ϕ(x, z + 2∆z, t). (3.5)
▷ Computation of the new phase 6:

Compute: ũz (x, z, t) = R F -1 (B(F(e iϕ )))(x, z, t) . (3.6)
▷ Application of low pass filter

7:
end for 8: end for 9: Compute:

u z (x, z, t) = π 2 -arctan(α(z -vt)) π ũz (x, z, t). (3.7)
▷ Application of arctan to recreate the no-wave area quantities, using a coefficient corr ef f , see (3.5). This coefficient is chosen to be random, which in practice provides better results. A heuristic explanation is that the smudges do not always propagate as they can be cut out by other smudges and aggregates, so it seems natural for the combination to be stochastic. This algorithm is detailed in Algorithm 3.

The Figure 4 displays a comparison between two images simulated by Algorithm 2 and Algorithm 3. At fixed time, both give similar results. However, Algorithm 3 seems to produce better time evolving signals, with more time regularity and the exhibition of smudges going downwards. A simulation takes roughly 0.5 seconds to be done with this algorithm. While it is still long, it is a significant improvement compared to the finite-element solver. Note that there are many different free parameters that we did not explicit. As a first approach, these are chosen by trials and fails, based on the visual appearance of the results.

Quantitative criteria

The method proposed in the previous section seems promising. This model have many parameters and until now, the fitting has only been made by visualization. Of course, this is not satisfying and we would like to dispose of quantitative criteria to evaluate the accordance of the signals generated by our model with respect to the simulations provided by the CEA.

In this section, we propose many dissimilarity measures and we compare them through their capacity to discriminate wave profiles provided by the finite-element solver according to the input parameters (frequency f of the incoming wave, proportion p and size r of the aggregates). Afterward, we use the best dissimilarity measure in the implementation of an optimization algorithm, namely the Nelder-Mead method, in order to fit the parameters of Figure 4: Two noise simulations using the dynamic method. The left image has been generated using Algorithm 2 and the right image using Algorithm 3. Both codes have been written in python, using scipy.fftpack for the discrete Fourier transforms. the generative model.

Dissimilarity measures

For a simulation U i (x, z, t) and a direction κ ∈ {x, z, t}, we denote by V κ i the (biased) variance along the axis κ. We also denote V x,z,t i the variance along the temporal and both spatial directions. Finally, we denote by V x i (z 0 ) the variance along the axis x at fixed z 0 . In this report, we propose and study the following dissimilarity measures:

d 1 (i, j) = |V x,z,t i -V x,z,t j | , d 2 (i, j) = ∥V t i -V t j ∥ 2 , d 3 (i, j) = ∥(V t i ) T V t i -(V t j ) T V t j ∥ 2 , d 4 (i, j) = ∥V t i (V t i ) T -V t j (V t j ) T ∥ 2 , d 5 (i, j) = ∥V x i (z 0 ) -V x j (z 0 )∥ 2 ,
where ∥•∥ 2 stands for the Frobenius norm in the appropriate vector space and (V x i ) T for the transpose of V x i . Because they only take the variances and covariances into account, these functions measure the difference of dispersion between two simulations. We could have considered other combinations of variances and covariances distances but it was not possible because of insufficient time.

We want to use one of these dissimilarity measures to fit our generative model. As a consequence, one should choose the measure with respect to its capacity to correctly discriminate the simulations provided by the finite-element solver according to the input parameters f , p and r. In this direction, we construct for each dissimilarity measure a dissimilarity matrix where the entry (i, j) is given by d(i, j), see an example in Figure 5. For each input parameter, there are ten simulations which are listed as follows:

• i ∈ {1, . . . , 10} correspond to parameters f = 5 kHz, p = 20 and r = 2; • i ∈ {11, . . . , 20} correspond to parameters f = 5 kHz, p = 5 and r = 4;

d 1 (i, j) d 2 (i, j) d 3 (i, j) d 4 (i, j) d 5 (i,
• i ∈ {21, . . . , 30} correspond to parameters f = 5 kHz, p = 10 and r = 4;

• i ∈ {31, . . . , 40} correspond to parameters f = 5 kHz, p = 20 and r = 4;

• i ∈ {41, . . . , 50} correspond to parameters f = 5 kHz, p = 40 and r = 4;

• i ∈ {51, . . . , 60} correspond to parameters f = 10 kHz, p = 20 and r = 4;

• i ∈ {61, . . . , 70} correspond to parameters f = 15 kHz, p = 20 and r = 4.

In addition, we compare the simulations to white noises filtered by a low pass filter, which are labeled by i ∈ {71, . . . , 80}. We convert the dissimilarity matrix into a similarity matrix by applying the operator x i,j → -x i,j + max k,ℓ x k,ℓ to each entry. To evaluate the capacity of a measure to discriminate the input parameters, we apply spectral clustering to these similarity matrices and we compute the corresponding adjusted rand index (ARI) and adjusted mutual information (AMI), see Table 1

. Except for d 1 (i, j) = |V x,z,t i -V x,z,t j
|, the scores obtained for the proposed measures are alike (around 0.55 for the ARI and around 0.75 for the AMI). This suggested that the d 2 , d 3 , d 4 and d 5 are more or less equivalent. In the sequel, we choose d 5 for fitting our generative model. This is motivated by the fact that it will be easier to extract the coda with this dissimilarity measure.

In Figure 6 is drawn the graph of the function t → V x (z 0 , t) averaged on the ten simulations samples for each input parameters. An interpolation procedure has been applied in order to normalize the time scale. The graph exhibits a long-time behavior corresponding to the coda whereas small times correspond to the wave front propagation which generates a totally different noise regime. The long-time behavior seems sufficient to distinguish different input parameters.

Parameter optimization

In this subsection, we explain how to fit the parameters of the generative model according to the simulations provided by the CEA with the input parameters f = 5 kHz, p = 20 and r = 4.

Because we are only interested on modeling the coda, we fit only with respect to the long-time behavior that is, for times between 50 and 90 in Figure 6. Then we use the L 2 -norm on this zone to perform an optimization on the parameters of Algorithm 3 using the Nelder-Mead algorithm, a classical method for non smooth problems. Because of the duration of the simulations, we had to limit ourselves to 600 steps. While this does not 7: Graphs of the function t → V x (z 0 , t) for real data and two simulations provided by the algorithm 3, one with non-optimal parameters and the other one after a numerical optimization.

allow for the convergence of the algorithm, Figure 7 shows that it provides a significant improvement in terms of curve-fitting. As expected, we observe an important discrepancy between our simulations and the data in small times. Also, in the last times, we observe a slope of the curves corresponding to our generative model. We believe that it is a boundary effect related to numerical issues, probably linked with the low pass filter.

This optimization procedure could be improved with another choice of optimization algorithm (Nelder-Mead algorithm was chosen for implementation convenience) or with more data corresponding to the coda. However, this work indicates that it should be possible to optimize the parameters of our algorithm to adjust any of the parameters provided by the CEA.

Conclusion and further topics of research

In this work, different approaches to sample random fluctuations of ultrasonic wave within heterogeneous materials have been proposed. A first approach based on Gaussian processes gave promising results but has been discarded for computational cost issues. A second approach, inspired by speckle simulations, gave better and faster results. The addition of a low pass filter in time and an enforced space correlation along the z axis have been a key to obtain realistic time-evolving simulations. Finally different metrics have been proposed in order be able to quantify the quality of the simulations. They have been shown not only to be able to discriminate between different wave and material parameters, but also to identify noise that is not the result of an ultrasound propagation. Thus this allowed us to perform an optimization on the parameters of our model in order to get a better fit on the CEA simulations.

Multiple directions can be taken for future works. We hereby mention some of them.

A first work would have to be done on the code in order to get faster simulations. The current code has been written in python, but other languages could be explored in order to fasten the Gaussian sampling as well as the computations of the discrete Fourier transforms, which are the two components which take most of the simulation time. Note that by nature the algorithm has a high potential for parallel computations. This would allow for a significantly faster parameter optimization and faster statistical studies afterwards.

For the moment we only considered one neighbour in the update (3.5), but studies have to be made in order to find which exact neighbours should be taken into account in order to properly retrieve the behavior of the downward smudges.

A Fourier transform could be performed on the data in order to identify some characteristic modes or other properties that could nourish the pool of quantitative criteria.

Other wavelet transforms could be considered, not only to establish metrics, but also to modify the low pass filter. A suitable wavelet basis could act as a selector on the simulation in order to recreate with more fidelity the smudges observed in the data.

Finally, while the front wave area is not the most useful for detection, this is where the coda comes from. Thus, working on a proper modelization of this area could lead to significant improvements in the quality of the rest of the simulation. In this area the noise is composed of small patches created by the shocks between the front wave and the aggregates. These are very short-time patches which means that a dynamic procedure as described here might not be the best lead for this part. In particular, the use of a Generative Adversarial Network (GAN) should be considered. To briefly describe GANs, they are Neural Networks (NN) composed of one NN which has the role of creating a random vector and another one which is supposed to establish if it comes from a real set of data or if this is a fake. The alternative optimizations of both NN allows for high quality creations by the first one. This process is commonly used in image or text generation. In our case, as we have multiple time steps to simulate, it would not be realistic to train a GAN to generate the coda for the full simulation (since we have very limited data). However, successive generations of a few images of front wave noise would be possible since one CEA simulation provides multiple data of this nature. Some work would then have to be done in order to smoothly go from one small sequence to another.

Figure 1 :

 1 Figure 1: Surface of concrete showing multiple aggregates (image from wikipedia.org).

Figure 3 :

 3 Figure 3: Example of simulation using Gaussian processes at one time step. Coding has been done in R.

2 : 3 :

 23 for 0 ≤ t ≤ T max -∆t do for all x, z do 4:

Figure 5 :

 5 Figure 5: Dissimilarity matrice for d(i, j) = ∥V x i (z 0 ) -V x j (z 0 )∥ 2 .There are ten samples for each input parameter. For i ∈ {1, . . . , 50} (resp. i ∈ {51, . . . , 60}, resp. i ∈ {61, . . . , 70}) the frequency is set to f = 5 kHz (resp. f = 10 kHz, resp. f = 15 kHz). The samples i ∈ {71, 80} correspond to white noises filtered by a low pass filter.

Figure 6 :

 6 Figure 6: Graphs of the function t → V x (z 0 , t) averaged on samples for each input parameter.

Figure

  Figure7: Graphs of the function t → V x (z 0 , t) for real data and two simulations provided by the algorithm 3, one with non-optimal parameters and the other one after a numerical optimization.

Table 1 :

 1 Comparison of dissimilarity measures. The ARI and AMI scores are averaged over 50 iterations of the spectral clustering algorithm.

	j)