Yann Bulteau

Nathalie Tarrat

Nadine Pébère

Corinne Lacaze-Dufaure

8-Hydroxyquinoline complexes (Alq3

published or not. The documents may come

J o u r n a l Na me 8-hydroxyquinoline complexes (Alq 3 ) on Al(111): atomic scale structure, energetics and charge distribution Yann Bulteau, a Nathalie Tarrat, b ‡ , Nadine Pébère, a and Corinne Lacaze-Dufaure a ‡ * 8-hydroxyquinoline (8Hq) is known to efficiently inhibit the corrosion of aluminium by forming metal-organic layers (8Hq forms complexes with aluminium atoms). In the present work, the atomic scale structure and the energetics of 8-hydroxyquinoline complexes (Alq 3 ) adsorbed on an aluminium surface are investigated by dispersion-corrected DFT computations. Two scenarios are considered: (i) an Alq 3 complex, previously formed in vacuum, is deposited on a flat Al(111) surface or (ii) three deprotonated 8Hq molecules (q) directly adsorb on a defective Al(111) surface presenting Al adatoms (Al-Al(111)). For the Alq 3 formation in vacuum, each addition of a q molecule on the Al atom stabilises the system, the oxidation state of the Al atom evolving from Al I in Alq to Al III in Alq 2 and Alq 3 . The subsequent deposition of Alq 3 on Al(111) leads to a strong bonding between the q molecules of the complex and the Al(111) surface, with a significant electron transfer occurring from the surface to the complexes (0.73 to 1.57 e). The formation on the metal surface of Alq 3 complexes via the adsorption of q molecules on an Al adatom leads to more stable structures than the ones obtained from direct adsorption of Alq 3 on Al(111). For the most stable adsorption conformation, the three q molecules are bonded to the Al adatom but only two are bonded to the aluminium surface. In that case, the total electron transfer from the Al-Al(111) surface to the q molecules is of 4.40 e and the electron transfer from the Al(111) surface to the Alq 3 -like species is 2.04 e. The structure, energetics and charge distributions data demonstrate an iono-covalent bonding between the q molecules and the Al atoms, in the complex as well as on the aluminium surface.

The 8-hydroxyquinoline (8Hq, presented in Figure 1) is a chelating compound commonly used in analytical applications for the detection and removal of metal ions 1 and in other technological fields such as optoelectronics where the deprotonated 8hydroxyquinoline species (q) is combined with aluminium atoms to form highly fluorescent complexes [2][3][4][5] (Alq 3 ) employed for manufacturing organic light emitting diodes (OLEDs) [6][7][8] . The 8Hq molecule has also been the purpose of numerous studies focused on the corrosion inhibition of aluminium and aluminium alloys and it has been shown to be an highly effective inhibitor. [9][10][11][12][13] In these latter studies, the authors demonstrated that the molecule reacts with aluminium to form metal-organic layers. Indeed, Lamaka et al. 10 identified by X-ray photoelectron spectroscopy (XPS) nitrogen-containing layers adsorbed on the Al surface while Soliman et al. 12 evidenced, on the basis of Fourier transform infrared spectroscopy (FTIR) and energydispersive X-ray spectroscopy (EDX) analyses, the formation on the metallic surface of species combining 8Hq and Al cations. However, the structure of these metal-organic species deposited on an aluminium surface has not been yet totally characterised at the atomic scale. This is the main aim of the present work.

Experimentally, the deposition under ultrahigh vacuum conditions of Alq 3 complexes on transition metal surfaces (Ag, Cu, Au and Co) has been the subject of several studies. [14][15][16][17][18] Hill et al. 14 showed by ultraviolet photoelectron spectroscopy (UPS) and ellipsometry that the deposition of Alq 3 complexes on Ag surfaces leads to the formation of polarised molecular layers with an interface dipole. Ino et al. 15 investigated the electronic structure of Alq 3 /Cu and Alq 3 /Au interfaces by UPS and two-photons photoelectron spectroscopy (2PPE) and showed that Alq 3 in its ground state interacts weakly with these metal surfaces. Wang et al. 16 used scanning tunneling microscopy (STM) to resolve the structure of the first monolayer of Alq 3 on Cu. At very low coverage, individual Alq 3 complexes were identified with significant surface mobility. When increasing the surface coverage, pairs and chain-like structures could be resolved, which stability was at-tributed to direct dipolar and van der Waals interactions. Still on Cu, a combined density functional theory (DFT) and non-contact atomic force microscopy (nc-AFM) study identified intermolecular interactions between 8Hq in molecular assemblies on the copper surface. 17 Droghetti et al. 18 combined desorption, 2PPE studies and electronic structure calculations to conclude that on Co surfaces, one deposited monolayer of Alq 3 complexes is composed both of chemisorbed (first layer) and physisorbed species, the latter filling the empty spaces.

The adsorption of 8Hq molecules or Alq 3 complexes on aluminium surfaces has given rise to only a few experimental and theoretical studies, to our knowledge. Still in the context of ultrahigh vacuum deposition, UPS investigations 19 and nondispersion-corrected DFT calculations 20 claimed that the electronic properties of free Alq 3 complexes are preserved upon adsorption on aluminium, suggesting a physisorption of the complexes. Dealing with the interactions between the Alq 3 complexes and aluminium, Mason et al. 21 showed that Al deposition on an Alq 3 thin film modifies the UPS spectra of the Alq 3 thin film. The authors concluded that there was a destructive reaction between the Alq 3 layers and the Al atoms, showing strong chemical interactions. Subsequent DFT work, this time including a van der Waals correction, has highlighted the chemisorption of the Alq 3 complexes on an Al surface 22 , later confirmed by Yanagisawa et al.. 23,24 A chemical bonding was also reported from DFT investigations on Alq 3 complexes adsorbed on cobalt 18,25 and magnesium 26 . In these studies, the presence of a dipole at the organic/metal interface has been evidenced whether the Alq 3 species is deposited on Al, Co or Mg.

In the framework of corrosion inhibition of aluminium in aqueous solutions, metal-organic species can be formed from 8Hq crystal dissolved in water. This formation can occur either in solution, the complex already formed subsequently adsorbing on the surface, or directly on the aluminium surface. The direct formation of metal-organic species on metallic surfaces has been previously studied for other corrosion inhibitors, e.g. by Kokalj et al. for imidazole species 27 or benzotriazole 28 deposited on copper and by Wang et al. for alkanethiols adsorbed on gold surfaces. 29 Concerning the 8Hq species, we demonstrated in previous computational studies that all the species (native 8Hq molecule, 8Hq tautomer, deprotonated and protonated 8Hq species) can chemisorb on Al(111). 30,31 The deprotonated 8Hq form is the most reactive, 30 leading to the formation of stable and compact organic layers that can inhibit the O 2 reduction reaction on the Al(111) surface. 31 In the present study, we investigate the structure of the metal-organic complexes on an aluminium surface. To this aim, we consider two scenarios: an Alq 3 complex previously formed in vacuum adsorbs on Al(111) or three q molecules directly adsorb on an Al adatom present on the aluminium surface.

In the first part of the paper, computational details are given with a description of the models and the formulas of the different energies discussed in this study. In a second part, the structure, energetics and charge distribution are investigated for Alq n (n = 1, 2 and 3) complexes in vacuum and for Alq 3 deposited or formed on Al(111). The adsorbed complexes issued from the two scenarios are compared. Finally, the conclusion is presented.

Computational details

To investigate the adsorption of 8-hydroxyquinoline species or complexes on aluminium, we have considered the adsorption of dehydrogenated 8-hydroxyquinoline molecules on an Al(111) surface. The term 'dehydrogenated' means hereafter that the O-H bond is dissociated and that the H atom of the hydroxyl group is removed. The choice of the Al(111) surface comes from the fact that, although aluminium develops a protective aluminium oxide thin film on its surface in neutral aqueous solution, this protective film is often defective or partially dissolved due to pH variations, leading to a direct interaction between the species in solution and the bare metallic Al(111) surface. Concerning the molecule, the dehydrogenated 8-hydroxyquinoline radical (Figure 2) was used to model the deprotonated molecule, these latter corresponding to the 8Hq form leading to the formation of the complexes. This method is commonly used to investigate the adsorption of anionic molecules on metallic surfaces 27,[32][33][34][35][36] in the framework of periodic calculations because it leads to similar results independently of the initial charge on the molecule. In the following, the dehydrogenated 8-hydroxyquinoline molecule will also be noted "q". In this work, we investigate the interaction between q molecules and Al atoms within different systems: (i) Alq n complexes: n q molecules (n = 1, 2, 3) in interaction with one Al atom, (ii) Alq 3 /Al(111): one Alq 3 complex adsorbed on a flat Al(111) surface, (iii) n(q)/Al-Al(111): n q molecules (n = 1 and 3) in interaction with one Al adatom on an Al(111) surface.

In the following, the Al-Al(111) notation refers to an Al(111) surface with an Al adatom on one hcp site. Total energies were calculated in the framework of the dispersion-corrected density functional theory (DFT-D). In order to be able to compare our results in vacuum and on surfaces, all the calculations presented in this study were performed using periodic boundary conditions, the latter being necessary to model surfaces. When studying isolated systems, we used simulation boxes large enough to simulate a molecule in vacuum, i.e. not affected by its periodic im-ages. Complexation, adsorption, interaction and deformation energies were computed. Negative values for these energies reveal an exothermic process whereas positive values mean an endothermic process. The z direction is the direction normal to the Al(111) surface. ∆z value corresponds to the distance between one atom and the average surface plane.

Isolated Alq n complexes (n = 1 , 2, 3).

The Alq n (n = 1 , 2, 3) complexes are composed of one Al atom and n q molecules, noted from A to C (see Figure 3). To study the isolated q molecule and Alq n (n = 1 , 2, 3) complexes, the species were placed in orthorhombic simulation boxes whose sizes were 22x21x20 Å 3 for q and Alq, 32x31x30 Å 3 for Alq 2 and Alq 3 .

In the following equations, E X is the total energy of the optimised system X in vacuum, where X is one Al atom (X = Al), the q molecule (X = mol) or an Alq n complex (X = cx). The energy E nmol SP is the total energy of the isolated system containing n molecules (n = 1, 2, 3) at the geometry in the complex (without the Al atom). Moreover, the energy E moli SP is the total energy of the isolated q molecule i (i = A, B, C) at its geometry in the complex.

The complex formation energy is determined from the total energy of the complex E cx from which are subtracted n times the energy of an isolated q molecule, n being the number of q molecules in the complex, and the energy of an Al atom:

E cx f orm = E cx -nE mol -E Al (1)
To quantify the interactions in the Alq n (n = 1 , 2, 3) complex, two interaction energies are calculated:

i) The interaction energy between the n q molecule and the Al atom:

E nmol/Al int = E cx -E nmol SP -E Al (2)
ii) The interaction energy between the n q molecules.

E inter-mol int = E nmol SP -∑ i=A,... E moli SP (3)
For instance, in the case of the Alq 3 complex, the interaction between the three q molecules in the complex is equal to:

E inter-mol int = E 3mol SP -E molA SP -E molB SP -E molC SP ( 4
)
The deformation of the n q molecule during the formation of the Alq n (n = 1 , 2, 3) complex is quantified by:

E nmol de f = ∑ i=A, ... E moli SP -nE mol (5)
For instance, in the case of the Alq 3 complex, the deformation energy of the three q molecules is equal to:

E 3mol de f = E molA SP + E molB SP + E molC SP -3E mol (6)
1.2 Adsorption of the Alq 3 complex on Al(111).

Once the isolated complexes studied, we investigated the adsorption of Alq 3 on an Al surface. To that aim, surfaces were modelled by an asymmetric slab, composed of three Al(111) layers in a 8 × 8 supercell. Periodic images were separated by about 21 Å of vacuum height in the direction normal to the surface to avoid interaction between periodic images in that direction. The Alq 3 complexes were deposited on the top layer of the slab with an empty space of at least 12 Å between the complexes on the surface. The atomic positions of the complex and of the two top layers of the Al slab were unconstrained. The bottom layer of the Al(111) slab was fixed at the bulk positions with a calculated bulk lattice parameter of 4.037Å (experimental value 37 : 4.05 Å).

In the following equations, E X is the total energy of the optimised system X in vacuum, where X is the Alq 3 complex (X = cx), the Al(111) slab (X = Al(111)), or the Alq 3 /Al(111) system (X = cx/Al(111)). In the latter, the Alq 3 complex is adsorbed on the Al(111) surface. E X SP (X = cx, Al(111)) is the total energy of the isolated system X constrained at its geometry in the Alq 3 /Al(111) system.

To characterise the adsorption on the Al(111) slab of Alq 3 complexes previously optimised in vacuum, the following energy are calculated :

i) The adsorption energy of the complex on the Al(111) slab:

E cx/Al(111) ads = E cx/Al(111) -E cx -E Al(111) (7)
ii) The interaction energy between the complex and the Al(111) slab:

E cx/Al(111) int = E cx/Al(111) -E cx SP -E Al(111) SP (8)
iii) The deformation energy of the adsorbed Alq 3 complex and of the Al(111) slab upon adsorption:

E cx de f = E cx SP -E cx (9)
E Al(111) de f = E Al(111) SP -E Al(111) (10)
1.3 Adsorption of n q molecules (n = 1 and 3) on Al-Al(111): n q/Al-Al(111) species. Two different systems were investigated, differing by the number n of q molecules (1 or 3) which were deposited on the Al adatom of the Al(111) slab. The distances between the periodic images of the adsorbed molecules were large enough to avoid any interaction between periodic images.

In the case of the adsorption of three q molecules on the adatom, in order to explore the conformations space around an energy minimum, ab initio molecular dynamics were performed, followed by DFT optimisations of selected geometries extracted from the MD trajectories. All the details of the MD simulations are given below.

In the following equations, E X is the total energy of the optimised system X, where X is the q molecule (X = mol), the Al-Al(111) slab (X = Al-Al(111)), or the system with n molecules adsorbed on the Al-Al(111) slab (X = nmol/Al-Al(111)). E X SP (X = nmol, Al-Al(111)) is the total energy of the isolated system X at the geometry in the nmol/Al-Al(111) system.

To characterise the adsorption on the Al-Al(111) slab of n q molecules, the following energy are calculated:

i) The adsorption energy of n q molecules on the Al-Al(111) slab:

E nmol/Al-Al(111) ads = E nmol/Al-Al(111) -nE mol -E Al-Al(111) (11)
ii) The interaction energy between the n q molecules and the Al-Al(111) slab:

E nmol/Al-Al(111) int = E nmol/Al-Al(111) -E nmol SP -E Al-Al(111) SP (12)
iii) The interaction energy between two subparts of the system to locally quantify the interactions:

E part1/part2 int = E nmol/Al-Al(111) -E part1 SP -E part2 SP (13)

DFT calculations

DFT calculations were performed using the VASP (Vienna Abinitio Simulation Package) program. [START_REF] Kresse | [END_REF][39][40] The PBE exchange correlation functional 41,42 was used. The interactions between electrons and ion cores were described by potentials generated by the Projector augmented-wave method (PAW potentials). 43,44 Dispersion corrections were included, with the addition of a dispersion term to the DFT energy in the Grimme's D2 approach. 45,46 . To perform the Brillouin zone integrations, a Monkhorst-Pack sampling 47 of the k-points in the reciprocal space was used, i.e. a gamma k-point sampling for the isolated species and a 2x2x1 kpoints sampling for the systems containing an aluminium surface and for the smearing of the Fermi surface, the Methfessel-Paxton method 48 with a width = 0.1 eV was employed. The electronic convergence criterion was 10 -6 eV with a cutoff energy of 450 eV. For geometry optimisations, the convergence was considered to be reached when forces on atoms were smaller than 5.10 -3 eV Å -1 .

Geometry exploration by ab initio molecular dynamics

In order to explore the conformations space of three q molecules adsorbed on Al-Al(111), ab initio molecular dynamics were performed (T = 500 K, Langevin thermostat [START_REF] Allen | Computer Simulation of Liquids[END_REF][START_REF] Hoover | [END_REF][51] , timestep = 1.5 fs, runs of 10 ps, E cut = 300 eV, electronic convergence criterion = 10 -4 eV, gamma k-point sampling). Along the trajectories, only the atoms of the q molecules were free to move and the atomic positions of the Al atoms of the Al-Al(111) slab were fixed.

Two initial geometries were chosen. The first one included direct interactions between the molecules and the surface while keeping the O atoms bonded to the adatom (MD1). In the second one, the molecules mainly interacted with the adatom (MD2).

The MD simulations were followed by DFT optimisation of snapshots taken every 320 steps (E cut = 300 eV, electronic convergence criterion = 10 -4 eV, convergence criterion on forces for optimisation = 5.10 -3 eV Å -1 , gamma k-point sampling). The most stable minima (geometry 10 in MD1 and geometry 17 in MD2) were reoptimised with more accurate computational parameters (E cut = 450 eV, electronic convergence criterion = 10 -6 eV, convergence criterion on forces for optimisation = 5.10 -3 eV Å -1 , 2x2x1 k-points sampling).

Charge distribution analysis

Atomic net charges were calculated using the Bader population analysis 52 . For a given atom (O, N, C or Al atoms), the atomic net charge is calculated as:

q atom = q X atom -q isolated atom ( 14
)
where q X atom is the number of electrons on the atom in the X species (X = mol, cx, cx/Al(111), mol/Al-Al(111), 3mol/Al-Al(111)) and q isolated atom is the number of electrons on the isolated atom.

The electron transfer during the complexation or adsorption is also calculated on the atoms:

∆q atom = q X atom -q mol atom ( 15
)
where q mol atom is the number of electrons on the atom in the isolated molecule.

On a X species (with X = cx, mol), the total charge transfer upon adsorption ∆Q X is calculated as:

∆Q X = Q X -Q isolated X (16)
where Q X is the charge on X species in the optimised X/slab system (X = mol or cx and slab = Al(111) or Al-Al(111)) and Q isolated X is the charge on the optimised isolated X species. For one selected case (most stable complex formed on the Al-Al(111) surface), to get more insights into the bonding of the q molecules on the Al-Al(111) surface, we plot the electron density variation ∆ρ:

∆ρ = ρ nmol/Al-Al(111) -ρ nmol SP -ρ Al-Al(111) SP (17)
where ρ nmol/Al-Al(111) is the electron density distribution of the optimised nmol/Al-Al(111) system, and ρ nmol SP and ρ

Al-Al(111) SP the electron density distributions respectively of the n molecules and the Al-Al(111) slab in vacuum at their geometry after adsorption.

Results

In the following, the structure and energetics of the Alq n complexes in vacuum are first presented. Then, the structure of Alq 3 complexes on an aluminium surface is investigated within two scenarios: an Alq 3 complex previously formed in vacuum is deposited on a flat Al(111) surface or three q molecules directly adsorb on an Al adatom present on the aluminium surface.

Isolated Alq n (n = 1 , 2 ,3) complexes

Alq n complexes can be formed by interaction between n q molecules and one Al atom. For n = 3, the Alq 3 complex exists in two isomers: meridional, noted Alq mer 3 , and facial, noted Alq fac 3 . 2,53,54 In the meridional isomer, the O atoms form a plane that is almost perpendicular to the plane formed by the N atoms.

In the facial isomer, O atoms and N atoms belong to two planes almost parallel and the O atoms are directed to one side of the complex while the N atoms are directed to the opposite side.

It has been shown that the meridional isomer is more stable in aqueous solution that the facial one. 2,53,54 In order to investigate the interactions in these complexes, we performed calculations on systems composed of one Al atom and one to three q molecules, i.e. Alq, Alq 2 and Alq 3 complexes. The optimised geometries are presented in Figure 3. It was not possible to place four or more molecules around the aluminium atom due to steric hindrance. Energies calculated for the optimised complexes are presented in Table 1, O-Al and N-Al bond lengths in Table 2, and net charges on the atoms and charges variation upon complexation in Table 3.

The formation energy E cx f orm takes into account the interaction energy E nmol/Al int between the n q molecules and the Al atom, the variation of the internal energy of the n molecules due to the deformation occurring during the complexation E nmol de f and the interaction energy between the n q molecules E inter-mol int in the resulting complex, i.e.

E cx f orm = E nmol/Al int +E nmol de f +E inter-mol int
. The Alq n (n = 1 , 2, 3) complexes formation is highly stabilising, with complex formation energies per molecule ranging from -4.54 eV for Alq fac 3 to -4.99 eV for Alq. The q molecules are strongly bonded to the Al atom, proof of this being the interaction energy between the q molecules and the Al atom which ranges from -5.02 eV per molecule for Alq fac 3 to -5.39 eV for Alq 2 . The van der Waals contribution to these energies is low (1.4-5.2% for the formation energies and below 1.5% for the interaction energies), this contri-bution increasing with the number of q molecules. The calculated interaction energies between the q molecules in the complexes E inter-mol int are positive (0.24 to 0.61 eV) and show globally repulsive interactions between them, despite the attractive van der Waals contribution. Concerning the structure of the complexes, the O-Al bonds are shorter than the N-Al bonds (1.844 ± 0.074 Å and 2.051 ± 0.154 Å respectively) in all the complexes, and these values are in agreement with previous theoretical 24 The oxidation state of the aluminium atom varies when evolving from one to two complexed q molecules but remains stable when a third molecule is added. Indeed, the charge on the Al atom in the complexes is +0.86 e in Alq and +2.43 ± 0.02 e in Alq 2 and Alq 3 (Table 3). It shows that the Al atom (3s 2 3p 1 ) tends to be in Al I or Al III oxidation states in the Alq and Alq n (n = 2, 3) complexes, respectively. This significant change in the aluminium electronic structure occurring when n goes from 1 to 2 is reflected by a higher variation of the interaction energy upon the addition of an other q molecule when going from Alq to Alq 2 (-5.44 eV) than when going from Alq 2 to Alq 3 (-4.27 and -4.41 eV for the f ac and mer isomers respectively). The electron transferred by the aluminium atom are located more on the oxygen atoms (0.29 to 0.39 e) than on the nitrogen ones (0.18 to 0.27 e).

This study of the Alq n (n = 1 , 2, 3) complexes in vacuum shows that the system is stabilised by each addition of a q molecule on the Al atom. When three molecules are complexed to the Al atom, the Alq mer 3 isomer is found to be the most stable one by 120 meV, ‡ close to the value of 210 meV calculated by Lima et al. (M06-2X/6-31+G(d,p)/SDD level of theory). 3 High electrons transfers are observed from the Al atom to the molecules upon complexation with two oxidation states of the Al atom in the systems (Al I in Alq and Al III in Alq 2 and Alq 3 ). It leads to a strong bonding between the Al atom and the q molecules. As experimental works proposed that q complexes are present on aluminium surfaces during the corrosion inhibition process by 8Hq molecules [9][10][11][12][13] , the adsorption of complexes on the Al(111) surface was investigated and the results are presented in the following. First, the deposition on Al(111) of an Alq 3 complex previously formed in vacuum was studied.

The geometries of the fac and mer isomers were reoptimised employing an implicit solvation model (implemented in the VASP code) to simulate the effects of the aqueous solution. Both isomers were found to be isoenergetic in solution, with a fac/mer populations ratio of 50/50. The adsorption of both the fac and mer isomers was thus considered in the following.

Table 1 Formation energies of the Alq n (n = 1 , 2, 3) complexes and van der Waals energy contribution, formation energies per molecule, interaction energies between the q molecules and the Al atom and van der Waals energy contribution, interaction energies per molecule, interaction energies between the q molecules and van der Waals energy contribution and total deformation energies of the n molecules (in eV) complexes were absorbed on an Al(111) surface. For both mer and fac isomers, initial geometries were chosen so that O and N atoms are directed towards the aluminium surface, in order to favour the bonding between the complex and the surface. When the complex is adsorbed on the Al(111) surface, it can access two configurations, i.e. an upward configuration, denoted by 'up' hereafter, and a downward configuration, designed by 'down' hereafter. These notations are the ones proposed by Yanagisawa et al. 24 who characterised the adsorption modes relatively to the molecular permanent dipole. The latter is directed to the vacuum side for the up configuration while it is directed to the substrate in the down configuration. The mer/up, mer/down, and the fac/up adsorption modes are show in Figure 4. The fac/down orientation was not investigated because it has been shown to be the least stable one.

Species

The coverage θ of the Al(111) surface by the adsorbed complexes is defined as the number of surface Al atoms covered by the complexes n covered The adsorption energies of the Alq 3 complexes on the Al(111) surface, the interaction energies, and the complex and Al(111) slab deformation energies are summarised in Table 4, together with the distances between the Al atom of the complex and the aluminium surface, the net charges on the Al atom of the complex and the charge variations due to adsorption.

The fac/up configuration is the most stable one, with the three O atoms of the Alq 3 complex bonded to the Al surface (top positions, d O-Al = 2.007 ± 0.002 Å). This is consistent with the previous works on the adsorption of Alq 3 on Al(111) 24 and on Co(111) 18,25 . The fac/up configuration is energetically followed by the mer/up and mer/down configurations that are higher in energy by 0.54 and 0.69 eV, respectively. In the mer/up configuration, the Alq 3 complex is bonded to the surface by two O atoms (top positions, d O-Al = 1.971 Å and 2.003 Å) whereas in the mer/down configuration, the Alq 3 complex bonds to the Al The complexes are strongly in interaction with the Al(111) surface with values ranging from -1.52 eV to -2.21 eV for the adsorption energies of the Alq 3 complexes on the Al(111) surface and -3.12 to -4.25 eV for the interaction energies. The Al(111) surface is deformed in the adsorption zone (E Al(111) de f = 0.57 ± 0.06 eV), and it is also the case for the complexes which deform in order to optimise their interaction with the surface. The deformation energies of the complexes are much lower for the up configurations (0.81 ± 0.07 eV) than for the mer/down one (2.09 eV) in which the bonding of the O, C, and N atoms with the surface induced a strong deformation. The mer/down configuration is thus the less stable adsorption mode (-1.52 eV) even if its complex-surface interaction energy is the strongest one (-4.25 eV).

The distance between the Al atom in the complex and the Al(111) surface is above 3.5 Å and the charge on the Al atom in the complex remains unchanged upon adsorption for the fac/up and mer/up configurations and is very slightly modified for the mer/down configuration (∆q Al = -0.09 e). Whatever the complex, an electron transfer occurs during adsorption, the electrons moving from the metal surface to the complex. The largest transfer takes place for the mer/down configuration (1.57 e), this significant charge transfer inducing, as explained above, strong complex-surface interactions and molecule deformation.

To investigate the stability of the adsorption of the Alq 3 complex on the Al(111) surface at ambient (300K) and high (500K) temperatures, molecular dynamics runs of 3 ps were performed for the fac/up configuration. The system stayed in a similar fac/up To summarise, we showed by depositing the Alq 3 complexes pre-formed in vacuum on the aluminium surface, that the most stable configuration is the fac/up one. Following the second scenario, which consists of the direct adsorption of three q molecules on an adatom of the aluminium surface, the formed complex could be more stable. In the next section we explore this path.

) q O (∆q O ) q N (∆q N ) q O (∆q O ) q N (∆q N ) q O (∆q O ) q N (∆q N ) Isolated q -0.

Direct formation of an Alq 3 complex on an Al-Al(111)

surface.

In the following paragraph, the adsorption of one q molecule on the Al-Al(111) surface is investigated, followed by the one of three q molecules. In all systems, the Al adatom is located on a threefold hcp site of the Al(111) surface as it was found to be 46 meV more stable than on a threefold fcc site. In that case, the adatom lies at 2.034 Å perpendicularly to the Al(111) surface plane and has a net charge of +0.12 e. For comparison, the DFT-D calculated Al(111) interlayer distance in the bulk is 2.314 Å. In the perpendicular structure, which lies 0.13 eV higher in energy (Table 5), the molecule is only bonded on top of the adatom (d O-Al = 1.813 Å, d N-Al = 1.998 Å). In both cases, the q molecule is thus strongly bonded to Al-Al(111) with adsorption energies of -4.10 ± 0.07 eV (in Table 5). For comparison, DFT-D adsorption energies of a single q molecule on a perfect Al(111) surface were computed and are -3.61 eV and -3.48 eV for a tilted and a parallel adsorption mode, respectively. The comparison of the adsorption energies thus indicates that the bonding of a q molecule on Al-Al( 111) is stronger than its bonding on a flat Al(111) surface. The complex formed on the surface in the parallel adsorption mode is stabilised by the very strong interaction between the molecule and the Al-Al( 111) surface (-6.14 eV for the interaction energy to be compared to -5.35 eV for the Alq complex in vacuum) whereas the q molecule is less bonded to the Al atoms for the perpendicular adsorption mode (-4.57 eV for the interaction energy). The q molecule is highly deformed in its parallel adsorption mode (1.42 eV for the deformation energy) as well as the Al-Al(111) surface which has a significant deformation energy (0.55 eV) due to a local deformation of the surface plane and the displacement of the adatom (+0.373 Å) which has occurred during the molecule adsorption. This is a first insight into the formation of an Alq species on the Al surface. For the perpendicular mode, the q molecule and the Al-Al(111) slab are much less deformed (0.40 eV and 0.14 eV) and the distance between the adatom and the Al(111) surface is only slightly decreased (-0.054 Å).

One q molecule adsorbed on

For both adsorption modes, the Al adatom is in its Al I oxidation state with a net charge of +0.82 e and +1.03 e, similarly to the case of the isolated Alq complex. These results show that the q-Al adatom part of the system behaves as the isolated Alq complex.

For the parallel adsorption mode, there is an electron transfer of 2.11 e from the metal to the q molecule. These electrons come from the atoms of the surface plane and from the Al adatom, which undergoes a charge variation of 0.70 e. For the perpendicular geometry, the q molecule gains 0.95 e that comes from the Al adatom. When a q molecule is adsorbed on a perfect Al(111) surface, the electron transfer from the Al(111) surface to the molecule is 1.92 e and 1.20 e for the parallel and tilted adsorption modes, respectively. Thus, for the parallel adsorption of the q molecule, the additional bonding to the adatom increases by 0. We demonstrated that the bonding of one q molecule on the Al-Al(111) surface is strong and that the parallel and perpendicular adsorption modes are close in energy, due to a balance Table 5 One q molecule adsorbed on Al-Al(111). Adsorption energies, interaction energies, molecule and Al-Al(111) surface deformation energies (in eV), distance (distance variation) between the Al adatom and the surface (in Å), net charge on the Al adatom (charge variation), charge variation on the q molecule upon adsorption (in e) between interaction and deformation energies. In the following paragraph, we investigate the adsorption of three q molecules on Al-Al(111).

Three q molecules adsorbed on Al-Al(111).

To explore the conformations space of a system containing 3 q molecules adsorbed on the Al-Al(111) surface, we performed two ab initio molecular dynamics runs (MD1 and MD2) differing by their starting geometry. In the initial geometry of MD1, the three q molecules were adsorbed parallel to the surface and bonded to the Al adatom, while for MD2 two q molecules lied parallel to the surface and bonded to the adatom and the third one was located on top of the Al adatom. Twenty geometries were extracted along each MD trajectory and optimised. Some selected distances are depicted in Figure 6 for the structures issued from MD1 and in Figure 7 In Figure 8, the total energy of the system is depicted relatively to the most stable one issued from MD1 (geometry 10). The values seem to be split into two groups: one corresponding to the lowest energies (geometries 1-2, 9-13, 15 and 17) and one corresponding to higher energies (geometries 3-8, 14, 16, 18-20). The energies of these two groups are different of at least 0.2 eV. For the most stable geometries, except for geometries 1 and 2, two Osurface Al atoms bonds are present. Moreover, one N-Al adatom and one or two N-surface Al bonds are formed. An exception is geometry 17, which has a low energy despite the absence of Nadatom bond. This latter seems to be compensated for by the presence of a bond between the Al surface and the C7 atom of molecule C. The structures extracted from the MD2 trajectory differ much less than those from MD1. Indeed, the evolution of the selected distances depicted in Figure 7 shows few variations of the position and orientation of the three molecules A, B and C, relatively to both the Al adatom and the surface. Two O atoms are always bonded to the surface (Figure 7(b)), two N atoms to the adatom (Figure 7(c)), and one N atom to the surface (Figure 7(d)). The only significant variations are related to the bonding between C atoms of the molecules and Al atoms of the surface (Figure 7(e) and (f)). In Figure 9, the energy is plotted relatively to the most stable form issued from MD2 (geometry 17). The energies are in a very small range (< 0.05 eV), except for geometries 2, 17 and 18. The geometry 2 is the less stable one. In that geometry, the C9 atom of the molecule B is in interaction with an Al atom of the surface. This destabilises the system, as the corresponding molecule is highly deformed (E molB de f = 1.82 eV). For the more stable geometries (17 and 18), the formation of one C7-surface bond (2.192 ± 0.001 Å) stabilises the system.

Most stable geometries for three q molecules adsorbed on Al-Al(111).

Geometry 10 from MD1 and geometry 17 from MD2 were reoptimised with more accurate computational parameters and Bader charge analyses were conducted on these structures. The geometry 17 from MD2 is the most stable one, the geometry 10 being 0.38 eV higher in energy. During the accurate optimisation of geometry 10, additional N-adatom and C7-Al bonds were formed. The structures of the systems are pre-sented in Figure 10. In both geometries, one molecule lies parallel to the metal surface, another is slightly tilted on the surface and the third one adopts an adsorption mode that is close to the perpendicular one. The interatomic distances between the oxygen and nitrogen atoms and the aluminium adatom are characteristic of covalent bonds. These latter are of the same order as the sum of the covalent radii for a simple bond (r O cov + r Al cov = 0.63 + 1.26 = 1.89 Å; r N cov + r Al cov = 0.71 + 1.26 = 1.97 Å) and largely smaller than the sum of the van der Waals radii (r O vdW + r Al vdW = 1.52 + 1.84 = 3.36 Å; r N vdW + r Al vdW = 1.55 + 1.84 = 3.39 Å). In the isolated Alq 3 complexes, the O-Al distance is around 1.9 Å and the N-Al distance is around 2.1 Å. Thus the adsorption of 3 q molecules on Al-Al(111) leads to the formation of O-adatom and N-adatom bonds shorter than the ones observed in the free complex. In addition, one of the nitrogen atoms is not bonded to the adatom but to an Al surface atom (molecule B).

The molecule A is, in both structures, the one that has a quasiperpendicular adsorption mode. It is in interaction with the Al adatom via its O and N atoms (O-Al and N-Al distances of 1.812/1.978 Å and 2.042/1.974 Å in geometries 10/17) and, in geometry 17, it is also in interaction with the surface via its O atom (1.947 Å).

The molecule B is the one that is positioned with a slight tilt relatively to the surface, the tilt being significantly larger in structure 10 than in structure 17. It is bonded to the adatom by its O atom (distance O-adatom: 1.887/1.882 Å for geometries The analysis of the electron population in geometries 10 and 17 is presented in Table 6 and gives information about the charge distribution within the systems and the bonding between the q molecules and the Al-Al(111) surface. In both geometries, the Al adatom is close to an Al III oxidation state, with a charge of respectively +2.56 e and +2.48 e in geometries 10 and 17, similarly to the net charge of the Al atom in the isolated or adsorbed Alq 3 complexes. A very strong electrons transfer occurs from the Al-Al(111) surface to the q molecules (respectively 4. The adsorption energy of the three q molecules on Al-Al(111) is -13.08 eV in the most stable geometry located during our conformational exploration (geometry 17 MD2). To locally characterise the strength of the interactions within the 3q/Al-Al(111) system, interaction energies between different subparts of the system were computed and the corresponding values are summarised in Table 7, together with the deformation energies. For this stable species on the aluminium surface, the interaction energy between the three q molecules and the Al-Al(111) surface is -17.20 eV, evidencing a very strong bonding of the molecules on the Al-Al(111) surface.

As described above, in geometry 17, the molecules A, B and C have different adsorption topologies. Even if its adsorption geometry is close to a perpendicular mode, the molecule A interacts strongly with the rest of the system, with an interaction energy of -6.40 eV (with a contribution of 14% of the van der Waals interactions), and is slightly deformed with a deformation energy of 0.49 eV during the formation of the 3 q/Al-Al(111) species. Similarly, the molecule C is slightly deformed upon ad- surface and is thus strongly deformed upon adsorption (deformation energy: 2.03 eV).

To investigate the strength of the interaction of the Al(111) surface with the Alq 3 -like species, the interaction energy between the three molecules (A + B + C) and the Al adatom, as in the isolated complex, and the interaction energy between the complexlike system (A + B + C + Al adatom) and the Al(111) surface are also presented in Table 7. The interaction energy between the three molecules (A + B + C) and the Al adatom has a value of -13.48 eV. It again demonstrates that the three molecules are strongly bonded to the Al atom, forming an Alq 3 complex with a geometry that is nearly as stable as the optimised complex when isolated from the Al(111) surface (E

Discussion

The direct adsorption of Alq 3 complexes on an Al(111) surface and the adsorption of three q molecules on an Al-Al(111) surface can only be compared in term of total energy and the relative total energies ∆E are presented in Table 6. The system formed by the adsorption of three q molecules on an Al-Al(111) surface is the most stable one. The Alq 3 complex adsorbed on Al(111) in the fac/up configuration is 0.25 eV higher in energy followed by the mer/up and mer/down configurations at 0.74 and 0.94 eV, respectively. In the most stable system, the Al adatom is at 2.738 Å from the Al(111) surface, i.e. about 1 Å smaller than the distances between the Al atom of the complex and the Al(111) surface in the mer/up, mer/down and fac/up geometries. The transition from the direct adsorption of Alq 3 complexes to the most stable config-uration (geometry 17) would therefore require a significant structural reorganisation. The charge on the central aluminium atom is very similar in the four structures, the Al atom being in strong interaction with the q molecules with an oxydation state close to Al III . With regard to the formation of Alq 3 complex layers on aluminium surfaces, it is conceivable that these two pathways may exist depending, inter alia, on the species present in solution (presence of free q or not), on the number of defects present on the aluminium surface, and on the activation barriers in play in the different pathways. The determination of these latter will be the purpose of forthcoming studies.

A very interesting point highlighted by this study is that the charge transferred from the metal surface to the metal-organic complexes is substantial, resulting in negatively charged complexes. This charge distribution raises the question whether such charge distribution is also observed at high coverage, i.e. in metal-organic layers covering aluminium surfaces. If this is the case and if the coverage of the metal by the organic layer is total, it could be envisaged that electrostatic repulsion could be part of the corrosion inhibition mechanism, as a globally negatively charged surface layer would repel the anions responsible for corrosion initiation.

Conclusions

The atomic scale structure of Alq 3 complexes adsorbed on Al(111) surfaces has been investigated following two scenarios: the direct adsorption of Alq 3 complexes pre-formed in vacuum on flat Al(111) surfaces, or the sequential adsorption of three molecules on a defective Al(111) surface, on which Al adatoms are present. During the Alq 3 formation in vacuum, each addition of a q molecule on the Al atom stabilises the system, the oxidation state of the Al atom evolving from Al I in Alq to Al III in Alq 2 and Alq 3 . Two conformations of the Alq 3 complex were investigated in vacuum and both were subsequently deposited on a flat Al(111) surface. The Alq 3 complex bonding to the surface is strong, with adsorption energies between -1.52 and -2.21 eV. This deposition lead to the formation of O-Al, N-Al and C-Al bonds between the q molecules of the complex and the Al(111) surface and to a significant electron transfer from the surface to the complexes ranging from 0.73 to 1.57 e. The formation of Alq 3 complex via the sequential adsorption of q molecules on an Al adatom present on the aluminium surface was investigated thanks to AIMD runs, which allowed the conformational space to be explored. The obtained adsorption conformations are more stable than the one obtained from direct adsorption on Al(111) of pre-formed Alq 3 complexes in vacuum. In the most stable conformation, two q molecules are almost parallel to the Al surface and the third one is almost perpendicular, the three molecules being bonded to the Al adatom but only two being bonded to the surface. The electron transfer from the Al-Al(111) surface is 4.40 e and is distributed unevenly on the three molecules with an enhanced transfer (2.18 e) to the molecule most bonded to the surface, via its O, N and C atoms.

It is conceivable that the two studied pathways leading to the formation of Alq 3 layers on aluminium surfaces could (co-)exist. Parameters such as the presence in solution of free q molecules, the number of defects on the aluminium surface, or the determination of the activation barriers at play in the different pathways would allow to go further in the comparison of these pathways. Based on structures, energetics and charge distributions data, one can conclude that the bonding between the q molecules and the Al atoms has an iono-covalent nature. The nature of this interaction will be studied in more detail in a forthcoming paper presenting a topological study of these bonds and the effect of the solvent on these adsorption modes will be evaluated. The significant charge transferred from the metal surface to the metal-organic complexes lead to the open question of a possible electrostatic contribution to corrosion inhibition mechanisms.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

Notes and references

Fig. 1

 1 Fig. 1 The 8Hq molecule.

Fig. 2

 2 Fig.2The q molecule. O atom: red ball; N atom: blue ball; C atoms: brown balls; H atoms: white balls.

Fig. 3

 3 Fig. 3 Isolated Alq n (n = 1 , 2, 3) complexes. The Alq 3 species exists in meridional (mer ) and facial (fac) isomers. A, B and C indicate the q molecules in the complex. Al atoms: big blue balls; O, N, C and H atoms: same colors as in Figure 2.

  and experimental2,3,55 studies. The shorter bond lengths in the Alq 2 species, which has a symmetric structure, show stronger O-Al and N-Al interactions than in the other Alq n species. The geometrical and energetic properties of the Alq 2 species are different from the ones of the other complexes with a stronger deformation of the molecules (0.52 eV per molecule for Alq 2 versus 0.36 eV and 0.29 ± 0.01 eV for Alq and Alq 3 ) and a strong repulsion between them (E inter-mol int = 0.24 eV).

3 and the Alq mer 3 complexes

 33 Aldivided by the total number of Al surface atoms in the supercell (n sur f ace Al = 64 atoms). A complete surface coverage (θ = 1 ML) corresponds to the case in which all surface Al atoms are covered by the complexes. As the Alq fac do not have the same geometrical form and can adsorb in an up or down position, their adsorption lead to the coverage of different numbers of surface Al atoms. Thus, θ = 16 ML, 0.19 ML and 0.20 ML for the fac/up, mer/up and mer/down conformations respectively.

  (111) surface by one O atom, one N atom, and one C atom (top positions, d O-Al = 1.957 Å, d N-Al = 2.059 Å and d C-Al = 2.156 Å).

  Al-Al(111). Two adsorption modes for a q molecule on Al-Al(111) were studied in which the molecule is either parallel or perpendicular to the surface. The geometry optimisations lead to two stable geometries, presented in Figure 5. In the most stable one, the q molecule is parallel to the Al(111) surface and is bonded both to the Al adatom (d O-Al = 1.927 Å, d N-Al = 1.914 Å) and to the Al(111) surface (d O-Al = 2.051 Å, d C-Al = 2.209 Å).

  2 e the electron transfer to the molecule. The comparison of the tilted adsorption on Al(111) and the perpendicular adsorption on Al-Al(111) is less obvious since the molecule is further from the surface when adsorbed on the adatom.

Fig. 4

 4 Fig. 4 Alq 3 complexes adsorbed on an Al(111) surface. Top: meridional up (mer/up). Middle: meridional down (mer/down) Bottom: facial up (fac/up); for practical purposes, only two layers of the Al(111) slab are shown.

  for the ones issued from MD2 and the associated energies are shown in Figures 8 and 9, together with the number of bonds existing between the q molecules and the Al-Al(111) surface. Three O-Al adatom bonds are present in all the structures,

Fig. 5

 5 Fig. 5 One q molecule adsorbed on Al-Al(111). Top: parallel adsorption mode. Bottom: perpendicular mode.

Fig. 6

 6 Fig. 6 Three q molecules adsorbed on Al-Al(111). Geometries optimised from 20 snapshots of MD1; (a) atom numbering; distances (in Å) (b) between the O atom and the Al adatom, (c) between the O atom and a surface Al atom (shortest distances), (d) between the N atom and the Al adatom, (e) between the N atom and a surface Al atom (shortest distances), (f) between the C7 atom and a surface Al atom (shortest distances). Data for molecules A in blue, B in red and C in green.

Fig. 7

 7 Fig. 7 Three q molecules adsorbed on Al-Al(111). Geometries optimised from 20 snapshots of MD2; for atom numbering, see inset; distances (in Å) (a) between the O atom and the Al adatom. (b) between the O atom and a surface Al atom (shortest distances), (c) between the N atom and the Al adatom, (d) between the N atom and a surface Al atom (shortest distances), (e) between the C7 atom and a surface Al atom (shortest distances), (f) between the C9 atom and a surface Al atom (shortest distances). Data for molecules A in blue, B in red and C in green.

10/ 17 )

 17 , and to the surface via its O and N atoms, with distances between O and N and the surface Al atoms of 1.956/1.926 Å and 2.067/1.930 Å for geometries 10/17. In addition, a C7-surface bond is present in geometry 17 (2.185 Å).Finally, the molecule C, which is parallel to the Al(111) surface, is bonded to the Al adatom both via its O and N atoms (distances O-adatom and N-adatom of 1.823/1.844 Å and 2.024/1.985 Å for geometries 10/17). In geometry 10, the N and C7 atoms are also bonded to surface Al atoms (distances N-surface Al atom and C7-surface Al atom of 2.183 Å and 2.153 Å).

  28 e and 4.40 e for geometries 10 and 17). The electron gain redistribution is different for the three molecules and depends particularly on their direct bonding with the Al(111) surface. The molecule A has a charge of -0.85 e for the geometry 10 and -0.97 e for the geometry 17 . Its perpendicular adsorption conformation makes a direct bonding with surface Al atoms difficult to form. However, one O-surface Al atom bond is observed in geometry 17 and the electron transfer is slightly higher. The q molecules bonded to the adatom and to Al(111) via their O (or N) and C7 atoms gains 2.18 e. It corresponds to C in geometry 10 and to B in geometry 17. The electron transfer to the third molecule, B in geometry 10 and C in geometry 17, is lower (1.25 e).

Fig. 8

 8 Fig. 8 Three q molecules adsorbed on Al-Al(111). Geometries optimised from 20 snapshots of MD1; relative total energies ∆E (in eV) in black (the reference is the energy of the geometry 10); number of bonds between O atoms and surface Al atoms in dark orange; number of bonds between N atoms and the Al adatom in green-blue; number of bonds between N atoms and surface Al atoms in light blue; number of bonds between C7 atoms and surface Al atoms in brown. Atoms a and atom b are considered bonded if d a-b < 2.3 Å.

Fig. 9

 9 Fig. 9 Three q molecules adsorbed on Al-Al(111). Geometries optimised from 20 snapshots of MD2; relative total energies ∆E (in eV) in black (the reference is the energy of the geometry 17); number of bonds between O atoms and surface Al atoms in dark orange; number of bonds between N atoms and the Al adatom in green-blue; number of bonds between N atoms and surface Al atoms in light blue; number of bonds between C7 atoms and surface Al atoms in brown; number of bonds between C9 atoms and surface Al atoms in beige. Atoms a and atom b are considered bonded if d a-b < 2.3 Å.

Fig. 10

 10 Fig. 10 Three q molecules adsorbed on Al-Al(111). Top: geometry 10 from MD1. Bottom: geometry 17 from MD2 (most stable form).

  3mol/Al int = -15.13 ± 0.07 eV in that case). The interaction energy between the Alq 3 -like species (A + B + C + Al adatom) and the Al(111) surface, with a value of -6.23 eV, evidences a strong adsorption of the complex on the Al(111) surface. This latter is significantly stronger that in the case of the direct adsorption on the Al(111) surface of the Alq 3 species pre-formed in vacuum (E cx/Al(111) int = -3.61 ± 0.49 eV). The electron density variation ∆ρ occuring during the adsorption of the 3 q molecules on the Al-Al(111) surface is presented in Figure 11 for geometry 17. A large redistribution of the electron density is observed and the electron accumulation between the molecules, the adatom and the aluminium surface confirms the formation of direct O-Al, N-Al and C-Al bonds. Thanks to this electron excess zone between the molecules, the Al adatom and the Al(111) surface, and the strong electrons transfer from the Al-Al(111) surface to the molecules highlighted above by Bader analysis, the bonding of the molecules with the metal atoms can be characterised as iono-covalent.

Fig. 11

 11 Fig. 11 Electron density variation upon adsorption for three q molecules adsorbed on Al-Al(111). Isosurface value: -0.0675/+0.0675 e• Å -3 ; yellow (blue): electron excess (deficit) regions.

  This work was granted access to the HPC resources of CALMIP supercomputing center under the allocation 2019-[p12174]. It was also supported by the French Ministry of Higher Education and Research [YB, PhD grant].

Table 2

 2 Bond lengths (in Å) in the Alq n (n = 1 , 2, 3) complexes (shown in Figure3)

		E cx f orm	E cx f orm (vdW )	E cx f orm /n	E int nmol/Al	E int nmol/Al	(vdW )	E	nmol/Al int	/n	E inter-mol int	E inter-mol int	(vdW )	E nmol de f
	Alq	-4.99	-0.07	-4.99	-5.35	-0.07		-5.35					0.36
	Alq 2 Alq mer 3 Alq fac 3	-9.51 -13.74 -13.62	-0.21 -0.64 -0.71	-4.75 -4.58 -4.54	-10.79 -15.20 -15.06	-0.14 -0.22 -0.22		-5.39 -5.07 -5.02		0.24 0.56 0.61	-0.09 -0.45 -0.52		1.04 0.90 0.83
	Species		mol A			mol B					mol C
			d O-Al	d N-Al		d O-Al		d N-Al		d O-Al		d N-Al
	Alq		1.848	2.187
	Alq 2 Alq mer 3 Alq fac 3		1.770 1.891 1.856	1.897 2.039 2.097		1.770 1.860 1.855		1.897 2.052 2.099		1.888 1.856		2.092 2.103
	2.2

Alq 3 complexes on an aluminium surface 2.2.1 Adsorption of Alq 3 complexes on the flat Al(111).

	The Alq fac 3	and the Alq mer 3

Table 3

 3 Net charges q atom and charge variations ∆q atom (in e) on the Al, O and N atoms during the formation of Alq n (n = 1 , 2, 3) complexes

	Species	Al	molA	molB	molC
		(q Al

Table 4

 4 Alq 3 complexes adsorbed on an Al(111) surface. Relative total energies (the reference is the total energy of the fac/up configuration), adsorption energies of the Alq 3 complexes on the Al(111) surface, interaction energies, complex and Al(111) slab deformation energies (in eV), distance between the Al atom of the complex and the Al(111) surface (in Å), net charge on the Al atom of the complex and charge variation on the complex upon adsorption (in e)

		∆E cx/Al(111)	E ads cx/Al(111)	E int cx/Al(111)	E cx de f	E	Al(111) de f	d z (Al)	q Al	∆Q cx
	mer/up	0.54	-1.68	-3.12	0.88	0.56	3.862	+2.44	-0.73
	mer/down	0.69	-1.52	-4.25	2.09	0.63	3.897	+2.36	-1.57
	fac/up	0	-2.21	-3.45	0.73	0.51	3.509	+2.45	-0.97
	geometry for all the duration of the runs without breaking of O-Al
	bonds.

Table 6

 6 Alq 3 species on Al(111). First two lines: 3 q molecules adsorbed on Al-Al(111) in geometries 10 and 17. Last three lines: Alq 3 complexes adsorbed on Al(111). Relative total energies (in eV, the reference is the energy of geometry 17); charge on the Al atom in the complex or the adatom; charge on molecules A, B and C (in e); distance from the Al atom in the complex or the Al adatom to the Al(111) surface (in Å)

	Species	∆E	q Al	Q molA Q molB Q molC dz Al
	Geometry 10 0.38 +2.56 -0.85	-1.25	-2.18	3.206
	Geometry 17 0	+2.48 -0.97	-2.18	-1.25	2.738
	mer/up	0.79 +2.44 -0.94	-1.26	-0.98	3.862
	mer/down	0.94 +2.36 -1.26	-1.49	-1.25	3.897
	fac/up	0.25 +2.45 -1.08	-1.19	-1.14	3.509

Table 7

 7 Most stable geometry for three q molecules adsorbed on Al-Al(111): interaction energies between subparts of the system, van der Waals energy contribution and deformation energies (in eV)

	Part 1	Part 2	E int part1/part2	E int part1/part2	(vdW )	E	part1 de f
	molA + molB + molC	Al-Al(111)	-17.20	-2.82		2.95
	molA	molB + molC + Al-Al(111)	-6.40	-0.91		0.49
	molB	molA + molC + Al-Al(111)	-6.91	-1.55		2.03
	molC	molA + molB + Al-Al(111)	-5.87	-1.30		0.43
	molA	Al-Al(111)	-4.56	-0.49		0.49
	molB	Al-Al(111)	-6.15	-1.23		2.03
	molC	Al-Al(111)	-4.78	-1.11		0.43
	molA + molB + molC	Al adatom	-13.48	-0.20		2.95
	molA + molB + molC + Al adatom	Al(111)	-6.23	-2.83

[1][2][3][4][5][6][7][8][9][10][11][12][13][14]

| 1-14

1-14 J o u r n a l N a me , [ y e a r ] , [ v o l . ] ,

J o u r n a l N a me , [ y e a r ] , [ v o l . ] ,


Pascale Domingo
email: domingo@coria.fr

Zacharias Nikolaou
email: z.nicolaou@cyi.ac.cy

Andréa Seltz
email: seltz@coria.fr.

Luc Vervisch
email: vervisch@coria.fr

From discrete

come

Introduction

Virtual prototyping has become an essential ingredient in the development and optimisation of combustion systems, as furnaces, boilers, internal combustion engines and gas turbines [START_REF] Bilger | Paradigms in turbulent combustion research[END_REF][START_REF] Gicquel | Large Eddy Simulations of gaseous flames in gas turbine combustion chambers[END_REF][START_REF] Poinsot | Prediction and control of combustion instabilities in real engines[END_REF][START_REF] Mastorakos | Forced ignition of turbulent spray flames[END_REF][START_REF] Locci | Selective non-catalytic reduction (SNCR) of nitrogen oxide emissions: A perspective from numerical modeling[END_REF]. Even relying on high performance computing, the spatial resolutions reached by the grids do not allow for fully resolving the smallest scales of the flow motions and scalar signals when simulating real combustion chambers. The introduction of modeling is therefore mandatory to deal with unresolved sub-grid scale (SGS) fluctuations affecting the transport of mass, momentum and energy. In addition, chemical reactions mechanisms are strongly non-linear and specific tools must also be developed to compute burning rates over grids which are coarse at the scale of the thin reaction zones observed in flames [START_REF] Bray | The challenge of turbulent combustion[END_REF][START_REF] Peters | Turbulent Combustion[END_REF][START_REF] Veynante | Turbulent combustion modeling[END_REF][START_REF] Pitsch | Large Eddy Simulation of turbulent combustion[END_REF].

In large-eddy simulation (LES), a spatial filtering operation is applied to the first principle equations of aerothermochemistry and any quantity φ(x, t) resolved by the grid reads

φ(x, t) = +∞ ∫ -∞ φ(x , t)G ∆ (x -x)dx , (1)
where G ∆ (xx) denotes a filtering operation of characteristic size ∆ and φ(x, t) is the fully resolved signal i.e., the solution that would be observed in a direct numerical simulation (DNS) resolving all the scales [START_REF] Lesieur | Large-Eddy Simulations of Turbulence[END_REF].

Because of unresolved SGS effects, the filtered value of any non-linear term T(φ) differs from T(φ). So far, closures were developed in the literature in the form of mathematical expressions, providing estimation of the sub-grid scale effects from the knowledge of the quantities resolved by the grid, then T(φ) = T (φ), where T (φ) denotes the physical model used to express T(φ) from the resolved field φ.

The application to turbulent flames of a different paradigm for modeling, which is based on a direct data analysis, is discussed in this chapter. Two main points are addressed:

• First, two methods of signal reconstruction based, respectively, on approximate inverse filtering and iterative deconvolution are discussed. In these two approaches, a tentative true signal φ (x, t) = L -1 ∆ φ(x, t) is approximated from φ(x, t), the resolved field, to compute the non-linear terms which are then explicitly filtered,

T(φ)(x, t) = T L -1 ∆ φ(x, t) . (2)
• Second, a neural network is trained to relate the filtered non-linear terms (i.e., T(φ)) to three-dimensional distributions of resolved quantities surrounding the considered grid point (i.e., T(φ) or any function of φ). A mapping function H is then obtained so that

T(φ)(x, t) = H T(φ(x 1 , t)), • • • , T(φ(x N , t) , (3)
with x j the N points selected around x to build an input image for a neural network previously trained.

Both approaches are evaluated a priori from different DNS databases, for the prediction of unresolved SGS convective and diffusive fluxes and for the modeling of filtered chemical sources.

Scalar signal deconvolution and machine learning for premixed turbulent flames 3

2 Flame filtering and signal reconstruction

Filtered scalar balance equation

The budget equations for scalars describing flames (chemical species mass fractions, energy, temperature) may be cast in the generic form,

∂ ρφ ∂t + ∇ • (ρuφ) = ∇ • (ρDφ) + ω , ( 4
)
where ρ is the density u the velocity vector, D a molecular diffusion coefficient and ω a chemical source. Applying the filtering operation (1) leads to

∂ ρ φ ∂t + ∇ • (ρ ũ φ) = ∇ • (ρD( φ)∇ φ) + ω + ∇ • τ , (5)
where φ = ρφ/ρ and

τ = ρD∇φ -ρD( φ) , ∇ φ SGS diffusion -ρuφ -ρ ũ φ SGS convection , (6)
are two unresolved SGS fluxes. The molecular diffusion coefficient usually depends non-linearly on the flow thermochemical properties and the notation D( φ) denotes D computed from the resolved fields.

Signal deconvolution

Two methodologies for signal reconstructions are now discussed after introducing a relation useful to determine the filter size which should actually be applied according to the mesh resolution and the number of points required within the flame.

Relation between flame thickness, filter size and mesh size

For a mesh of characteristic size h and a laminar premixed flame thickness δ L , the relation between the number of points within the filtered signal and the filter size ∆ may be approximated as [START_REF] Nikolaou | A priori assessment of an iterative deconvolution method for les sub-grid scale variance modelling[END_REF]:

h ∆ = 1 n π 6 + δ 2 L ∆ 2 , (7)
where a Gaussian filtering operation has been assumed Fig. 1 Temperature profiles across a one-dimensional premixed flame. Line: T (x); Square: T (x) from deconvolution (Eq. ( 10)); Dash: T (x); Circle: T (x). Solving balance equation for filtered 1D flame from deconvolution Triangle: T (x). ∆/δ f l = 8. Reprinted with permission [START_REF] Domingo | Large Eddy Simulation of premixed turbulent combustion using approximate deconvolution and explicit flame filtering[END_REF].

G(x) = (6/(π∆ 2 )) 3/2 exp(-6x • x/∆ 2 ) , (8) 0
with a flame front modeled with an error function. This relation shows that for a LES with ∆ = 5δ L , the mesh size should verify h = 0.75∆/n. Since n > 1, h is therefore expected smaller than the filter size applied to the scalar fields.

Approximate deconvolution and explicite filtering (ADEF)

Most filtering operations may be approximated as a diffusive process formulated implicitly φ(x, t) = φ(x, t) + d∇ 2 φ(x, t) ,

with d = ∆ 2 /24 for a Gaussian filter [START_REF] Katopodes | Explicit filtering and reconstruction turbulence modeling for large-eddy simulation of neutral boundary layer flow[END_REF]. From this relation, an approximate deconvolution operator is readily obtained in which all derivatives are resolved over the coarse mesh (h < ∆):

φ (x, t) = L -1 ∆ φ(x, t) = φ(x, t) -d∇ 2 φ(x, t) . (10)
To illustrate these filtering/deconvolution operations, Fig. 1 shows T(x) (line), the temperature signal through a one-dimensional unstrained stoichiometric methane/air premixed flame, along with T(x) (dashed line) and T(x) (circle). The deconvoluted signal T (x) (square) perfectly matches T(x). Moreover, the solving of the budget equations of the one-dimensional filtered flame with all the non-linear terms computed from approximate deconvolution and filtering, returns the expected T(x) (triangle). This was done here for ∆/δ L = 8 with a second order interpolation of the scalar signals inside computing cells (h = 0.5∆) [START_REF] Domingo | Large Eddy Simulation of premixed turbulent combustion using approximate deconvolution and explicit flame filtering[END_REF]. Flame deconvolution with an appropriate sub-grid interpolation or a regularisation procedure thus appears as a robust tool [START_REF] Wang | Regularized deconvolution method for turbulent combustion modeling[END_REF][START_REF] Mehl | Evaluation of deconvolution modelling applied to numerical combustion[END_REF]. Sometimes deconvolution is also associated to a scale similarity hypothesis [START_REF] Vreman | A similarity sub-grid model for premixed turbulent combustion[END_REF]. Approximate deconvolution is easily generalised to three-dimensional turbulent flames. Figure 2 shows snapshots of the simulation of a stoichiometric methane/air turbulent bunsen flame studied experimentally [START_REF] Chen | The detailed flame structure of highly stretched turbulent premixed methane-air flames[END_REF]. These simulations are performed with the flow solver SiTCom [START_REF] Bouheraoua | Large-eddy of a supersonic lifted jet flame: Analysis of the turbulent flame base[END_REF], considering the Navier-Stokes equations in their fully compressible form together with a balance equation for c, a filtered progress variable (c = 0 in fresh gases and c = 1 in burnt gases [START_REF] Bray | The challenge of turbulent combustion[END_REF]). Tabulated detailed chemistry [START_REF] Gicquel | Laminar premixed hydrogen / air counterflow flame simulations using flame prolongation of ILDM with differential diffusion[END_REF][START_REF] Van Oijen | Modeling of complex premixed burner systems by using flamelet-generated manifolds[END_REF] is coupled to deconvolution of this progress variable with Eq. ( 10) to simulate the premixed turbulent jet flame. The tabulation is based on the GRI-3.0 [21] detailed chemistry, with a progress variable defined from CO, CO 2 , H 2 O and NO x as in [START_REF] Godel | Tabulation of nox chemistry for large-eddy simulation of non-premixed turbulent flames[END_REF]. The thermochemical quantities are then read from a lookup table φ = φ T (c) and the filtered quantities may be written from deconvolution and explicit filtering,

ρ φ = ρ φ T (c ) = L -1 ∆ [ρ]φ T L -1 ∆ [ρ c] L -1 ∆ [ρ] . (11)
The convective terms are discretised with a fourth-order centered skew-symmetriclike scheme [START_REF] Ducros | High-order fluxes for conservative skew-symmetric-like schemes in structured meshes: application to compressible flows[END_REF] and the diffusive terms with a fourth-order centered scheme. Time is advanced explicitly with a third order Runge-Kutta method and NSCBC boundary conditions [START_REF] Lodato | Three-dimensional boundary conditions for direct and large-eddy simulation of compressible viscous flows[END_REF] are imposed at inlet and outlet, with the measured profiles with synthetic turbulence [START_REF] Klein | A digital filter based generation of inflow data for spatially developing direct numerical or large eddy simulations[END_REF] prescribed at inlet.

The time sequence (∆t = 1.8ms) of iso-burning rate (Fig. 2) illustrates the high level of resolved flame wrinkling in this jet flame simulation featuring u /S L = 5, with u the characteristic velocity fluctuation and S L the laminar flame speed. The formation of pockets at the flame tip is also observed. The mesh resolution1 is of the order of 170 µm. The approximate deconvolution and explicit filtering approach coupled with tabulated chemistry allows for simulating the flame and predicting species mass fractions without resorting to any adjustment of model parameters. Results are given in Fig. 3. The experimentalist reported a measurement error range between 8% and 15% for major species concentration and up to 25% for minor species [START_REF] Gicquel | Laminar premixed hydrogen / air counterflow flame simulations using flame prolongation of ILDM with differential diffusion[END_REF][START_REF] Van Oijen | Modeling of complex premixed burner systems by using flamelet-generated manifolds[END_REF]. The major species O 2 , CH 4 , H 2 O and CO 2 agree with measurements within this range (Fig. 3-right). The chemistry tabulation limited to a progress variable as input does not include strain rate as a parameter in these simulations, explaining the higher departure for the minor species CO, H 2 and OH.

Iterative deconvolution and explicit filtering (IDEF)

As an alternative to approximate deconvolution based on discretised inverse filtering, the Van Cittert algorithm is another well-established iterative deconvolution algorithm [START_REF] Van Cittert | Zum einfluss der spaltbreite auf die intensitätverteilung in spektralinien. II[END_REF][START_REF] Jansson | Deconvolution with applications in spectroscopy[END_REF]. The calculation of derivatives is avoided by expressing the deconvoluted scalar field from

φ n+1 = φ n + b(φ -φ n ) . ( 12
)
The initial condition of the iteration is φ 0 = φ and b is defined so that the deconvoluted variable stays within specific bounds [START_REF] Nikolaou | A priori assessment of an iterative deconvolution method for les sub-grid scale variance modelling[END_REF]. If φ l and φ h are the low and high bounds, on option consists in writting:

b =              φ n -φ l φ a -φ l if φ ≤ φ a φ n -φ h φ a -φ h otherwise (13)
with the average value φ a = (φ l + φ h )/2. The number of iterations is monitored through a volume averaged error e v = |φ -φ | , so that iterations can proceed as long as e v decreases, in practice for well resolved premixed flames (∆ < 10δ L ), the iterations never exceed a few decades.

Correlations between conditionally averaged progress variable obtained from (12) compared to DNS values are shown in Fig. 4. The data are from a DNS of a methane/air planar premixed flame interacting with homogeneous turbulence [START_REF] Nikolaou | Direct numerical simulation of complex fuel combustion with detailed chemistry: physical insight and mean reaction rate modeling[END_REF]. In this DNS, u /S L = 14, the turbulent Reynolds number is Re T = 230, the Damköhler number Da = 1.17, and the Karlovitz number Ka = 13.

These three-dimensional direct simulations have been performed using the SENGA2 code [START_REF] Cant | Senga2 user guide. cued/a?thermo/tr67[END_REF], which solves the compressible reacting flow Navier-Stokes equations using a 10th order finite difference scheme for interior points, and a 4th order Runge-Kutta scheme for the time-stepping. The DNS-filtered variables were sampled onto a LES (h/∆ = 0.25) mesh using high-order Lagrange polynomials, values which are then used for the deconvolution [START_REF] Nikolaou | A priori assessment of an iterative deconvolution method for les sub-grid scale variance modelling[END_REF]. The strong correlation between the prediction and the DNS reference (Fig. 4) confirms the potential of the iterative procedure for signal reconstruction. The time-averaged Pearson correlation coefficients are above 0.9850 for the three filter sizes (∆ = δ L , 2δ L and 3δ L ) considered. In order to better quantify the quality of the deconvolution process but also to circumvent any biases which may accompany correlation coefficients as defined above, mean percentage errors are also calculated using, e(x, y) = 1

N s x,y,z,t 100 • |x -y| |y| ( 22
)
where N s is the total number of sample points. The above error estimate includes contributions both from under-predicting and over-predicting a target variable y using a model variable x, and is therefore a stringent error measure. In addition, in order to avoid any On this basis, the iterative deconvolution was further applied to compute unresolved SGS fluxes, as the convective one [START_REF] Nikolaou | Scalar flux modelling in turbulent flames using iterative deconvolution[END_REF]

τ u = ρuφ -ρ ũ φ = (ρu) (ρφ) ρ -ρ ũ φ . ( 14
)
Evaluation against the same turbulent planar flame DNS configuration [START_REF] Nikolaou | Direct numerical simulation of complex fuel combustion with detailed chemistry: physical insight and mean reaction rate modeling[END_REF] of this SGS flux projected in the direction normal to the mean flame front is shown in Fig. 5.

In this figure, comparisons against the well-established Clark model [START_REF] Clark | Evaluation of sub-grid scalar models using an accurately simulated turbulent flow[END_REF] are also added, this physical model reads

f 1x = τ u • e 1 = ρ ∆ 2 12 ∂ ũ1 ∂ x j ∂ φ ∂ x j . ( 15
)
As expected the transition from counter gradient SGS turbulent transport to the usual gradient turbulent flux is found in the data by varying the turbulence intensity [START_REF] Veynante | Gradient and counter-gradient scalar transport in turbulent premixed flames[END_REF]. Indeed, f 1x > 0 (counter gradient transport) is seen in Fig. 5(a) for the progress variable larger than 0.3, and, f 1x < 0 (gradient transport) is observed in Fig. 5(b). The transport through the flame which was mainly driven by gas expansion for u /S L = 3.18 is rather governed by turbulent mixing for u /S L = 9. Interestingly, the deconvolution procedure perfectly reproduces this transition, while the physical model given by the relation [START_REF] Mehl | Evaluation of deconvolution modelling applied to numerical combustion[END_REF] misses this very basic, but essential, property of SGS turbulent transport in premixed turbulent flame. Following the same line, iterative deconvolution was also found effective to model the SGS stress tensor in premixed turbulent V-flame [START_REF] Nikolaou | Unresolved stress tensor modeling in turbulent premixed v-flames using iterative deconvolution: An a priori assessment[END_REF].

Machine learning for turbulent combustion modeling

DNS embedded in LES database

An additional DNS database was obtained in the configuration of the premixed turbulent bunsen flame discussed above in § 2.2.2 [START_REF] Chen | The detailed flame structure of highly stretched turbulent premixed methane-air flames[END_REF]. This DNS is located downstream of the well-resolved LES of the first sections of the piloted premixed stoichiometric fuel-air round jet (Fig. 6), LES that provides turbulent flame conditions for the DNS inlet plane [START_REF] Domingo | Dns and approximate deconvolution as a tool to analyse onedimensional filtered flame sub-grid scale modeling[END_REF][START_REF] Cifuentes | Local volumetric dilatation rate and scalar geometries in a premixed methane-air turbulent jet flame[END_REF][START_REF] Cifuentes | Effects of the local flow topologies upon the structure of a premixed methane-air turbulent jet flame[END_REF]. In this simulation, both LES and DNS are run simultaneously.

This dual simulation was achieved by embedding, inside the LES mesh, a zone where the resolution is sufficiently high to fully resolve the thin reaction zones and the smallest turbulent scales. The nozzle diameter is D = 12 mm, the jet Reynolds number is 24,000 (bulk nozzle velocity of 30 m • s -1 and turbulent kinetic energy of 3.82 m 2 • s -2 ). The pilot is composed of fully burnt gases at T b = 2200 K. The LES mesh contains 171 million nodes for a domain 16D×8D×8D, leading to a resolution of the order of 150 µm. The resolution in the DNS zone is uniformed at 50 µm, (a) the actual flux in the (inhomogeneous) direction of mean flow x, i.e., f 1 , and (b) the divergence of the flux, F = ∂f i /∂x i . The scalar flux in the mean direction of flow propagation, f 1 , is the dominant term since this is also the mean direction of flame propagation and there are large gradients in the progress variable c in this direction. The divergence of f , i.e., F , allows us to (a) examine the collective ability of the model in modeling the scalar-flux terms in the two other homogeneous directions y and z as well, and (b) determine whether the gradients of the scalar flux are accurately recovered in addition to its actual local magnitude. Note that in an actual LES the term (-F ) would be the one contributing to the right-hand side of the transport equation for the filtered progress variable [START_REF] Ducros | High-order fluxes for conservative skew-symmetric-like schemes in structured meshes: application to compressible flows[END_REF]. The data are averaged in the homogeneous directions y and z at every point in the direction of the mean incoming turbulent velocity field x. For any variable v we have

v x (i) = 1 N y N z j,k v(x i ,y j ,z k ), (20)
where N y and N z are the numbers of LES mesh points in the y and z directions, respectively, as given in Tables II and III for the different filter widths. A correlation coefficient is also used in order to quantify the local spatial agreement between the modeled variable and the actual variable as obtained from the DNS. For any two variables x and y this is defined as

R(x,y) = N i=1 (x i -xi )(y i -ȳi ) N i=1 (x i -xi ) 2 N i=1 (y i -ȳi ) 2 , ( 21
)
where the overbar in Eq. ( 21) denotes the sample mean and not filtering. The correlation coefficients for any two quantities are calculated at every time step and are time averaged over successive instances in order to increase the statistical accuracy of the results. Figures 3(a) and 3(b) show instantaneous values of the scalar-flux component f 1 as predicted using the IDEF method and Clark's model against the DNS values for the lowest turbulence level, case A, for + = 1. Both the DNS values and the models' predictions are normalized using the maximum instantaneous f 1 value as obtained from the DNS. It is clear that there is a strong correlation for both the IDEF-based model and the Clark model. Figures 4(a to secure a full resolution of the flow and flame scales for the chosen operating conditions [START_REF] Domingo | Dns and approximate deconvolution as a tool to analyse onedimensional filtered flame sub-grid scale modeling[END_REF].

Both LES and DNS rely on chemistry tabulation with a stoichiometric premixed flamelet with fresh gases at T o = 300 K (GRI-3.0 mechanism [21] and progress variable defined as in [START_REF] Godel | Tabulation of nox chemistry for large-eddy simulation of non-premixed turbulent flames[END_REF]), without any SGS modeling in the DNS part. The flame thermal thickness based on the progress variable field is of the order of δ L ≈ 400 µm. The mesh in the DNS zone is composed of 28.58 million nodes (243×343×343), over a physical domain of 12 mm×18 mm×18mm. The DNS starts at 4.5D downstream of the nozzle, and at 5.5D the mesh is coarsened again to progressively evolve toward LES resolution (Fig. 6). Modeling based on a progress variable presumed pdf approach is used for the chemical sources and species in the LES zones [START_REF] Domingo | Large-Eddy Simulation of a lifted methane-air jet flame in a vitiated coflow[END_REF] and the SGS momentum fluxes are approximated with the Vreman model [38].

As for the jet flame LES of § 2.2.2, these simulations are performed using the flow solver SiTCom [START_REF] Bouheraoua | Large-eddy of a supersonic lifted jet flame: Analysis of the turbulent flame base[END_REF]. A Gaussian filtering operation (Eq. ( 8)), with filter size ∆ = 0.3 mm = 0.75 δ L , ∆ = 0.6 mm =1.50 δ L and ∆ = 0.9 mm = 2.25 δ L , is applied to the DNS variables in order to generate a priori LES filtered quantities. Thus the resolution of the a priori fields varies from well-resolved to coarse LES from the reaction zone point of view (∆ = 0.9 mm is 18 times larger than the DNS grid resolution). The Q-criterion is defined as

2 Q = ( ij ij -S ij S ij ) = -(∂ u j /∂ x i )(∂ u i /∂ x j )
, where ij and S ij are the anti-symmetrical and the symmetrical parts of the velocity gradient tensor, respectively. Q > 0 indicates zones of flow rotation [74] . Figure 5 shows a zoom of the wrinkled flame surface interacting with the vortical structures. The progress variable iso-surface c = 0 . 8 , corresponding to the peak reaction zone, is seen in orange. In grey, an iso-Q is displayed for the locations where c > 0.8 (burnt gases). These iso-Q and iso-c surfaces show that the flame interacts, at the same time, with large flow scales generated in the fuel-jet shear layer and with turbulence inside the fresh gases. The flame surface in the DNS sub-zone is displayed in Fig. 6 , with the wrinkling of the flame over a quite large range of scales. A pocket of fresh gases detaching from the main flame surface appears, a pocket that will lead to an isolated island of fresh gases burning as they are convected downstream in the burnt products. At some other locations, the protruding parts of the iso-surface can merge and a pocket filled with burnt gases is formed at this particular instant in time ( Fig. 7 ).

The one-dimensional flame assumption is supposedly valid when the local reaction zone is convected, and propagates, with the flow as a whole, thus keeping its inner structure unchanged. There has been much discussion in the literature concerning the exact impact of velocity fluctuations on local flame properties according to the ratio between u ′ , the amplitude of the velocity fluctuations, and S L , the flame speed and this particular point is out of the scope of the present work. However, some very basic properties of the three-dimensional scalar field may be examined and related to the behaviour of the error brought by the one-dimensional flame deconvolution.

surface leading to burnt gases pocket

Pocket burning

Large scale wrinkling

Small scale wrinkling Fig. 6. Iso-c at 0.8 in the DNS zone.

In the DNS subzone, filtering and deconvolution operators discussed above are applied in an a priori manner and conditional statistics are collected, averaging over the DNS computational domain. Figure 8 shows snapshots of the progress variable source term filtered at different levels from = 0 . 75 δL = 0 . 3 mm up to = 3 . 00 δL = 1 . 2 mm. As expected, the filtering decreases the peak values and makes the reaction zones thicker. For large filter sizes, this is even followed by a change in flame topology, as observed in the merging of an isolated pocket, seen by comparing Fig. 8 (e) with other snapshots for smaller filter sizes down to DNS.

The DNS field c ( x , t ) is filtered at a length scale , applying the relation [START_REF] Domingo | Large Eddy Simulation of premixed turbulent combustion using approximate deconvolution and explicit flame filtering[END_REF] , to obtain an a priori c ( x , t ) field. Then, the one-dimensional flame deconvolution L -1

, L ( c ( x , t )) ( Eq. ( 17) ) is applied and the progress variable distribution is compared to the original DNS field c ( x , t ) ( Fig. 1 ). The result of this exercise is visible in Fig. 9 , displaying a transverse cut of the progress variable in the jet flame. Overall the flame structure is recovered, with however a tendency of L -1

, L ( c ( x , t )) to provide iso-contours of progress variable which stay mostly parallel, as expected with a one-dimensional flame model. In zones where the flame is highly wrinkled and where the local flame thickness features strong variations, this implies that the one-dimensional flame deconvolution misses the local flame structure.

A direct measure of the error introduced by the onedimensional flame hypothesis is defined as:

e ( x , t ) = c( x , t ) -L -1 , L [ c ( x , t )] c( x , t ) , (19)
and statistical conditional mean of this error are computed within the DNS zone,

E(φ ) = ⟨ e ( x , t ) | φ⟩ , (20)
where φ denotes the variable used to condition the statistical mean. φ = c ( x , t ) leads to E( c ) , the conditional mean of the error introduced by the one-dimensional flame deconvolution versus the filtered progress variable, which is shown in Fig. 10 . It is seen that the error peaks in the preheat zone, with a maximum that can reach up to 60% for = 3 δL . The snapshot of e ( x , t ) illustrates in Fig. 15 (b) the development of the error on the fresh gas side, but also the fact that e ( x , t ) is larger where the curvature is more pronounced, with an amplitude that depends on the sign of the curvature, which is a quantity that has been the subject of numerous studies in premixed turbulent combustion [75][76][77][78][79][80][81] . In the DNS zone, the flame surface is highly convoluted (Fig. 6-right), with flame wrinkling covering a large spectrum of length scales, including complex flame topological properties such as pockets of fresh or burnt gases.

Figure 7 shows ω + | c , the statistical mean over the DNS domain of the normalised filtered progress variable source, conditioned on values of c. The subscript ' + ' denotes source terms normalised by their maximum value in the reference freelypropagating laminar premixed flame. The result obtained using a 1D laminar flame, ω + ( c), is also shown as a solid line. The maximum of ω + | c decreases while increasing the filter size and thus while increasing unresolved fluctuations [START_REF] Bray | The challenge of turbulent combustion[END_REF]. Following the thickening of the filtered flame front, the response of this conditional filtered source term also spreads in progress variable space for increasing filter sizes, up to ∆ = 0.9 mm. Figure 7 clearly demonstrate that ω + ω + ( c) and not only the amplitude, but also the shape of the burning rate versus the progress variable, vary with the filter size. A simple re-scaling by a factor would therefore not be sufficient to properly model ω from ω( c). Similarly, Fig. 8 shows the conditional mean of the sum of the SGS convective and diffusive fluxes (Eq. ( 6)), which also requires modeling.

Machine learning for turbulent combustion modeling

The underlying idea discussed in this section consists of finding a mapping function G, so that the filtered chemical source ω(x, t) is constructed from the three dimensional distribution surrounding x of the chemical sources computed from the resolved field, i.e.,

ω(x, t) = G ω( c(x 1 , t)), • • • , ω( c(x N , t)) , ( 16
)
where N is the number of LES mesh cells involved in the mapping. Similarly, a second mapping function F is sought to express ∇ • τ, the divergence of the SGS fluxes (Eq. 6), from the three-dimensional distribution of resolved diffusive budgets ∇ • (ρD c ( c)∇ c) surrounding the x location, P. Domingo, Z. Nikolaou, A. Seltz, L. Vervisch During the training phase, a series of convolution/sampling operations are performed iteratively, in which the neural weights are adjusted until a satisfying minimal error is found between the value of the labels and the neural network prediction. A series of convolution/sampling operations are applied to the database, to extract its features using a number of different kernels [START_REF] Schmidhuber | Deep learning in neural networks: An overview[END_REF][START_REF] Lecun | Deep learning[END_REF]:

• Each image is first convoluted with 32 different filter kernels obtained with random values from a truncated normal distribution. Meaningful values of the returned features are extracted with a max pooling non-linear function to avoid excessive computational costs. • The process is repeated with 64 filters, decomposing the image into several most meaningful features, needed for seeking out the inner properties of the fluxes and sources. • Two fully connected layers are built to process the 64 obtained features, and to classify the image by giving the probability it belongs to each label.

The training of the network was conducted using the TensorFlow (www.tensorflow.org) library and breakdown of the network structure is given in Fig. 10. The training error function is based on cross entropy [START_REF] De Boer | A tutorial on the cross-entropy method[END_REF] and this training is controlled by the Adam optimizer [START_REF] Kingma | ADAM: A method for stochastic optimization[END_REF] for stochastic gradient descent, with a learning rate of 10 -4 . In using the networks, the N = 27 values of the chemical sources and of the divergence of the fluxes computed from the resolved progress variable field in the test box surrounding the LES cell (Fig. 9), are the input. For each filter size, 1000 filtered DNS fields are used for a priori tests (the noised images of the training phase do not enter these tests).

Figure 11 shows the comparison between the conditional means of the predicted divergence of the fluxes and their DNS counterpart. For the trained filter size, the CNN returns the expected value, for the untrained one some departure is observed, but compared to what can be observed in some case using a physical model, results are more than encouraging. (Compare for instance these results against those of Fig. 5 with the Clark model.)

The filtered chemical source are also very well predicted both the shape and the amplitudes are recovered by the CNN (Fig. 12). Interestingly, even the data corresponding to the filter sizes that were not used for training (0.4 mm and 0.45 mm) are quite well reproduced.

Conclusion

Large datasets from direct numerical simulation have been mostly used so far to explore reacting flow physics and to test well-established or novel modeling ideas. With rapidly growing computing facilities, followed by the development of well-resolved large eddy simulations for real combustion systems, the possibility of relying on modeling based on a direct analysis of the signals is now emerging.

Instead of expressing unclosed terms from algebraic relations involving quantities resolved by the mesh, deconvolution may be used to directly reconstruct the nonlinear terms to close the scalar balance equations, then some form of sub-grid scale interpolation is required. Both discretised and iterative deconvolutions have been found valuable along these lines.

In addition, the capabilities of machine learning tools to estimate unclosed terms from the resolved ones have been demonstrated, which opens many perspectives such as the coupling between advanced computational flame dynamics tools and artificial intelligence.

The deconvolution of scalar signals and the direct mapping form LES resolved scales to unknown filtered terms from CNN do not operate at the same level (signal reconstruction or direct prediction) and could thus be combined as complementary tools.

Fig. 2

 2 Fig.2Time sequence of iso-filtered source of progress variable, ∆t = 1.8 ms. Reprinted with permission[START_REF] Domingo | Large Eddy Simulation of premixed turbulent combustion using approximate deconvolution and explicit flame filtering[END_REF].

Fig. 3

 3 Fig. 3 Mass fractions in %. Symbol measurements [17]. Line: LES. In each subfigure: Left: Red: CO×10. Black: H 2 ×100. Blue: OH×75. Right: Blue: O 2 . Green: CH 4 . Red: H 2 O. Black: CO 2 . Reprinted with permission [26].

Fig. 7

 7 Fig. 7 Conditionally averaged progress variable as obtained from IDEF based on the actual value as obtained from the DNS for a case A, and b case B

Fig. 4

 4 Fig.[START_REF] Mastorakos | Forced ignition of turbulent spray flames[END_REF] Conditionally averaged progress variable as obtained from IDEF based on the actual value as obtained from the DNS. ∆ + = ∆/δ L Reprinted with permission[START_REF] Nikolaou | A priori assessment of an iterative deconvolution method for les sub-grid scale variance modelling[END_REF].

FIG. 7 .FIG. 8 .

 78 FIG. 7. Averaged (in homogeneous directions y and z) flux term f 1x for case A: (a) + = 1, (b) + = 2, and (c) + = 3. 043201-11

FIG. 9 .Fig. 5

 95 FIG. 9. Averaged (in homogeneous directions y and z) flux term f 1x for case C: (a) + = 1, (b) + = 2, and (c) + = 3. 043201-12

Fig. 5 .

 5 Fig. 5. Zoom iso-Q at 10,0 0 0 s -2 . Orange: iso-c at 0.8. are computed with a fourth-order centred skew-symmetric-like scheme [70] . The diffusive terms are approximated with a fourthorder centred scheme. Time is advanced explicitly with a third order Runge-Kutta method and the 3D-NSCBC boundary conditions [71] are imposed at inlet and outlet. To resolve eventual stiff scalar gradients, in the LES zones an artificial dissipation scheme was added introducing fourth order dissipative terms [72] (turned off in the DNS part). The measured profiles are prescribed at the inlet of the jet, with an addition of synthetic turbulence [73] .

Fig. 6 Fig. 7

 67 Fig.6LES-DNS snapshot of the jet-flame simulation[START_REF] Domingo | Dns and approximate deconvolution as a tool to analyse onedimensional filtered flame sub-grid scale modeling[END_REF]. Mesh and iso-progress variable c = 0.8. h: resolution. Right: zoom of iso-c = 0.8 in the DNS zone (different angle view). Reprinted with permission[START_REF] Domingo | Dns and approximate deconvolution as a tool to analyse onedimensional filtered flame sub-grid scale modeling[END_REF].

Fig. 8 ∇

 8 Fig. 8 ∇ • τ | c vs c. Filter size ×: 0.3 mm, : 0.6 mm, •: 0.9 mm. Reprinted with permission [39].

Fig. 11 ∇

 11 Fig. 11 ∇ • τ | c vs c. Symbols: DNS reference. Line: CNN prediction. Reprinted with permission [40].

  Fig. 12 ω + | c vs c. Symbols: DNS reference. Line: CNN prediction.

Table 4

 4 Global time-averaged Pearson correlation coefficients

A resolution of 50 µm would be necessary to fully resolve the flame (i.e., DNS) with tabulated chemistry and between 10 µm and 80 µm to resolve the Kolmogorov scale.

Filter DNS

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t >

!(e c(x 1 , t))

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t >

!(e c(x N , t))

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t >

!(x, t)

< l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t >

These two mapping function are built from convolutional neural networks (CNN), which conveniently for interpolating over a large number of reference values. Typically, if a set of N L reference values of filtered chemical sources ω[i] and divergence of the SGS fluxes ∇ • τ[i] are known and associated through the mappings to well defined properties of the resolved field (here sources and diffusive budgets computed form the resolved scalar fields), well trained neural networks will return Ω i (x, t) T i (x, t) as the weights to be affected to the i-th values contributing to the filtered quantities, i.e., ω(x, t) =

The ω[i], ∇ • τ[i], Ω i and T i are calibrated by training the neural network from the DNS database. For a given filter size ∆, a three-dimensional test-box of size (2h) 3 is built around every of the M = 28.58 million DNS nodes, where h is defined from Eq. [START_REF] Peters | Turbulent Combustion[END_REF]. Figure 9 illustrates this process. The test box is centered at x and contains N = 27 points which hold the three-dimensional distributions of ω( c(x j , t)) and ∇ • (ρD c ( c)∇ c) (x j , t), for j = 1, • • • , N. These data are stored and constitute the 'images' that will be processed by the CNN. From the DNS, the 'labels' of each i-th image are ω

Two networks of similar structures (same number of layers, convolution kernels, etc.) are trained, one for the chemical source and one for the SGS fluxes. For each value of ∆, the following procedure is applied: Fig. 10 Typical structure of a convolutional network for scalar deconvolution. More details on the TensorFlow functions used may be found in [START_REF] Seltz | Direct mapping from LES resolved scales to filtered-flame generated manifolds using convolutional neural networks[END_REF]. Reprinted with permission [START_REF] Seltz | Direct mapping from LES resolved scales to filtered-flame generated manifolds using convolutional neural networks[END_REF].

• First, 1000 images with their associated ω[i] and ∇ • τ[i] values are drawn. 20 values of c * uniformly distributed between 0 and 1 (∆c * = 0.05) are defined.

For each value of c * , 50 images are selected randomly such that c(x, t) ∈ [c * -∆c * /2; c * + ∆c * /2] (x is at the center of the test box, see Fig. 9). • To avoid overfitting, uncorrelated random perturbations are added to the images, ω( c(x j , t)) and ∇ • (ρD c ( c)∇ c) (x j , t) for j = 1, • • • , N = 27, as 10% of their maximum in the test box, thus a second image is available for each label. N L = 2000 images then exist for 1000 labels.

Using two filter sizes, the database to train the network contains 4000 images and 2000 labels and for every quantity learned, a set of 27×4000 =108000 data ( ω( c(x j , t)) and ∇ • (ρD c ( c)∇ c) (x j , t)) is involved, data which are associated to the 2000 reference labels ( ω[i] and ∇ • τ[i] for i = 1, • • • , N L ).