Adrian Ungureanu

Alexandru Chirieac

Carmen Ciotonea

Irina Mazilu

Cezar Catrinescu

Sabine Petit

Eric Marceau

Sébastien Royer

Emil Dumitriu

Enhancement of the dispersion and catalytic performances of copper in the hydrogenation of cinnamaldehyde by incorporation of aluminium into mesoporous SBA-15 silica

Keywords: copper nanoparticles, mesoporous materials, aluminium doping, hydrogenation, cinnamaldehyde

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

Introduction

Over the last two decades, copper-based nanostructured catalysts have attracted significant attention for hydrogenation reactions. 1 However, the control over the size of copper nanoparticles (NPs) is often difficult, especially on mesostructured supports, when standard methods such as impregnation are used to introduce the metal. Populations of large copper NPs, with low dispersion, poor catalytic activity and poor stability, are usually reported in the final catalysts. [2][3][4][5][6][7][8][9][10] Several strategies have been tested to stabilize copper nanoparticles (NPs) at sizes below 10 nm in mesostructured silicas: thermal treatments carried out at low pressures, 11,12 modification of silicas by APTES grafting, 2,13-15 impregnation using organic solvents 16 or organometallic complexes, 17 procedures based on self-induced assembly 18 or deposition-precipitation, 2,[19][20][21][22] and decomposition of inorganic precursors by autocombustion. 23 When standard impregnation procedures are used, it has been shown that the addition of second metals to copper (e.g., Ni, Cr), [24][25][26][27] or the functionalization 28,29 or pore occlusion of the mesostructured silica support by surfactants, 27,30 promote copper dispersion, with a resulting increase of catalytic activity in the hydrogenation of cinnamaldehyde. The hydrogenation of cinnamaldehyde (CNA) is known to be challenging as the unsaturated alcohol, cinnamyl alcohol (CNOL), should be selectively produced via C=O bond hydrogenation. However, the preferred path on most catalysts, including silica-supported Cu systems, is the hydrogenation of the C=C bond owing to favorable thermodynamics and lower activation barriers, [30][31][32][33][34] While the literature on noble metals is abundant, fewer attempts have been made so far to develop a selective copper-based catalytic system. 26,29,30,[33][34][35][36][37][38][39][40] It was reported in particular that Cu + or Cr + sites adjacent to the active metallic Cu 0 centers can function as electrophilic centers that polarize the C=O bond via the lone electron pair in oxygen, thus improving the hydrogenation of the carbonyl group in CNA to CNOL, yet at low reaction rates. 26,29 Previously, Marchi et al. reported that ternary Cu-Zn-Al and quaternary Cu-Ni(Co)-Zn-Al catalysts were one order of magnitude more active than Cu/SiO2 in CNA hydrogenation, while selective towards CNOL. 34 In these systems, Cu 0 particles appeared as highly dispersed over a super-stoichiometric zinc aluminate spinel, and in close interaction with M 2+ cations. On the other hand, it has been demonstrated that the introduction of Al atoms at the support surface has a positive effect on both the activity and chemoselectivity of Pt-based catalysts in the hydrogenation of CNA. 41 The presence of Lewis sites onto which the CNA molecule adsorbs via its C=O bond was key to improve the selectivity toward the unsaturated alcohol, without the intervention of a second transition metal. It can thus be tempting to solve simultaneously the problem of catalytic activity on the one hand, by promoting copper dispersion thanks to strong metal-support interactions, and the problem of selectivity toward cinnamyl alcohol on the other hand, by providing surface Lewis sites with enhanced affinity for the C=O bond. To this end, ordered mesoporous silica (OMS) supports can be modified by addition of aluminium (ordered mesoporous silica-aluminas, OMSA). A facile way to incorporate high amounts of aluminum into the structure of SBA-15 mesoporous silica is to combine the direct synthesis of silica in aqueous medium with Al introduction via a two-step pH adjustment method. 42 This approach was successfully used to prepare Al-SBA-15 with a wide range of Si/Al ratios (130-5.2). 43 Compared with Herein, we show that even if introduced by standard impregnation, monometallic copper can be stabilized as highly dispersed nanoparticles on water-stable Al-SBA-15 OMSA. We compare these catalysts with materials prepared on SBA-15 OMS and ordered mesoporous alumina (OMA) supports, to highlight the positive effect of aluminum incorporation on the catalyst performances, and we demonstrate that the combination of highly dispersed copper nanoparticles and presence of Al results in catalysts displaying a high activity and selectivity to cinnamyl alcohol.

Results and Discussion

Structural, textural and surface characterization Al-SBA-15 supports (OMSA10, OMSA20) were prepared with Al2O3 contents of 10 and 20 wt.% (Si/Al atomic ratio = 7.7 and 3.4, respectively), while SBA-15 silica (OMS) and organized mesoporous alumina (OMA) were used as reference supports (experimental details are supplied in the final section of the article). The low-angle XRD patterns of OMSA supports (ESI, Figure S1) exhibit one well-resolved peak and two weak peaks assigned to (100), ( 110) and (200) reflections associated with hexagonal p6mm symmetry of well-ordered mesostructured SBA-15. In line with our previous results, 43 it can be seen that upon Al introduction via the pH adjustment method, the short-and long-range mesopore ordering and pore size of Al-SBA-15 materials are well preserved at high Al loadings, because the concentration of Al species is high enough to protect the mesoporous structure against hydrolytic degradation. 5 wt.% copper was introduced by incipient wetness impregnation with an aqueous solution of copper nitrate on OMS and OMSA, or by melt infiltration of hydrated copper nitrate in the case of OMA. The choice of the latter procedure was made necessary by the high sensitivity of mesoporous aluminas to exposure to water [START_REF] Bejenaru | [END_REF]48 ICP-OES analysis confirmed that Al2O3 and Cu contents were close to the targeted values (Table 1). Besides, the Cu/OMS pattern presents intense and narrow diffraction peaks in the range 2 = 30 -80°, characteristic of monoclinic CuO (ICDD 048-1548) (Figure 1a). These CuO crystallites are large (mean size of 31.5 nm, estimated by the Scherrer equation), and thus located on the external surface of the support given the pore size of SBA-15 (8.4 nm). With the introduction of aluminium in the support (Cu/OMSA10), the diffraction peaks assigned to large CuO particles become much less intense and broader (Figure 1b), indicating that the average size of these crystallites has decreased (17 nm, see Table 1).

For the sample with the highest amount of aluminium (Cu/OMSA20), no diffraction peaks can be observed. Copper phases are thus highly dispersed, at a size below the detection limit in XRD (~ 3 nm). These results might be considered as first evidence of the positive role of the Al2O3 coating in improving the dispersion of copper species upon impregnation, due to an enhanced interaction between the OMSA support and copper species. 44 Very weak diffraction peaks are equally observed for Cu/OMA. Textural properties of the supports and of copper-loaded materials were investigated by N2 physisorption at -196 °C.

Figure 2 shows the isotherms and pore size distributions for copper catalysts (characteristics of the bare supports are shown in ESI, Figure S2), whereas Table 1 summarizes the textural properties of all materials. For Cu/OMS, sharp adsorption/desorption branches specific to highly-ordered materials are observed (isotherm of type IV with H1 hysteresis loop 45,[START_REF] Rouquerol | Adsorption by Powders and Porous Solids: Principles, Methodology and Applications[END_REF] ), and lead to a narrow pore size distribution centered at 8.2 nm, characteristic of bare mesostructured SBA-15 (Figure 2a). The decrease of specific surface area and pore volume after metal deposition is explained in part by the mixing of SBA-15 with non-porous CuO, that furthermore is present as bulky particles which do not clog the porosity. This interpretation is confirmed by TEM-EDXS analysis (ESI, Figure S3, showing large aggregates of extraporous CuO particles). A very low copper content is measured in the silica grains (< 0.5 wt.% by EDXS).

The hysteresis loops of the two Cu/OMSA materials are also well defined. Compared with Cu/OMS, OMSA-derived samples show a decreased specific surface area and pore volume, butworthy of emphasis -even at high loadings of alumina, the narrow pore size distribution (centered at 9.1 and 8.4 nm for Cu/OMSA10 and Cu/OMSA20, respectively) is maintained (Table 1, Figure 2c andd).

On the other hand, Cu/OMA displays a physisorption isotherm of type IV with a hysteresis loop of type H2, indicating an ordered mesoporous material, though less structured than OMS. [START_REF] Bejenaru | [END_REF] Surface area and pore volume were preserved after copper introduction by the melt infiltration method, that avoids the use of water as solvent.

Nature and reducibility of the copper phases

The chemical and electronic characteristics of the copper phases were investigated by XPS for the Cu/OMSA20 sample, after calcination and after reduction at 500°C. Cu 2p core level spectra are presented in Figure 4. The calcined sample shows the Cu 2p3/2 and 2p1/2 doublet peaks centered at 933.3 and 952.8 eV, respectively, with a spin-orbit splitting (SOS) of 19.6 eV. The absence of a peak located at 932.4 eV indicates that there are no bulk-like copper oxide NPs. [START_REF] Moulder | Handbook of X Ray Photoelectron Spectroscopy: A Reference Book of Standard Spectra for Identification and Interpretation of Xps Data[END_REF] A shoulder at 935.4 eV on the 2p3/2 peak suggests a charge transfer from metal ions toward the support, and the formation of species in which copper is in strong interaction with the support. 44,[52][53][54] A small satellite peak, arising from the 3d-4s transition during the relaxation process, 55 is observed at 944.5 eV, (11.2 eV from the main peak). XPS thus indicates the existence of Cu 2+ species in two different chemical environments: the contribution at BE = 933.3 eV is attributed to highly dispersed CuO, while the value at BE = 935.4 eV is assigned to copper in a CuAl2O4 spinel phase. The low ratio between the intensities of the main peak and satellite peak (Isat/Ipp = 0.07) may also be assigned to highly dispersed CuOx species in octahedral coordination 56 (in contrast with a value of 0.55 for bulk CuO 52,53 ). A correlation of this low Isat/Ipp value with the existence of Cu + species would not be relevant because no signal specific to Cu + /Cu 0 species at ~ 932 eV can be detected for the calcined sample. [57][58][59] After reduction at 500 °C, the binding energy is shifted to a lower value as expected for metallic Cu (932.8 eV) while

A B C

the satellite peak disappears, revealing that all the copper species are in a reduced form. 44 Thermo-programmed reduction (TPR) was used to determine the reducibility of the copper species and confirm their characteristics in terms of dispersion and interaction with the support (Figure 5). The TPR profile for Cu/OMS exhibited a single hydrogen consumption peak with a maximum at 286°C, corresponding to the reduction of bulk CuO. 60 This interpretation is in line with the XRD and TEM results obtained for this material (large CuO particles outside the SBA-15 porosity).

When the aluminum loading increases in the sample, the reduction profile becomes more complex and reduction temperatures shift to higher values. These changes can be related to the CuO particle size, 61 or to interactions with the Al2O3-functionalized support surface. Due to the progressive evolution of the chemical nature of the support, the reduction behavior of copper species supported on -alumina was thus considered as a reference, with copper aluminate reducing at the highest temperature. 62 The most complex profile belongs to Cu/OMSA10, with multiple reduction steps. Three maxima of reduction could be defined, at 220 °C, 308°C and 410 °C, respectively assigned to the small amount of external bulk-like oxide particles, highly dispersed copper oxide particles in the porosity, and copper species in strong interaction with the surface (CuAl2O4). 60,61,[63][64][65][66] The profile of Cu/OMSA20 resembles to that of Cu/OMA sample, suggesting similar characteristics of Cu species in the two materials. The main reduction peaks are observed at 310 and 390 °C (Cu/OMSA20), and 340 and 410 °C (Cu/OMA). Most copper is thus in strong interaction with the support, as dispersed CuO or copper aluminate, expected to predominate at low Cu loading and for calcination temperatures up to 500 °C. 65 For Cu/OMSA20, a shoulder located at 220 °C could originate from the reduction of residual CuO in weak interaction with the support.

Dispersion of copper active phase

The dispersion of metallic copper after reduction at 500°C was determined by N2O chemisorption. [66][67][68][69] According to TPR data, this temperature is consistent with the reducibility of the copper species and will be chosen for reduction of copper prior to the catalytic reaction. Results are expressed as copper dispersion (DCu), average particle size (dCu), and metallic copper accessible surface area (SCu) (Table 2). The dispersion of copper and the associated metallic surface area both increase as the aluminum content in the support increases. In particular, the values found for Cu/OMSA20 are very similar to those obtained for Cu/OMA, confirming that on OMSA20, copper particles present characteristics close to those on mesoporous alumina.

Catalytic properties

The catalytic properties of the above-discussed copperbased solids were evaluated in the liquid phase hydrogenation of cinnamaldehyde (CNA). As shown in Scheme 1, this process involves parallel and consecutive reactions, with the possible formation of partially hydrogenated molecules (hydrocinnamaldehyde -HCNA; cinnamyl alcohol -CNOL) and of a fully hydrogenated product (hydrocinnamyl alcohol -HCNOL).

Cu/OMS shows a very low activity (~9 %) after 6 h of reaction (Figure 6 A). It is related to the poor dispersion of copper and consequently to a low exposed active surface. 26,34 Activity increases with the increase of Cu exposed surface area (Table 2), and thus with the aluminum content in the support. Cu/OMSA10 and Cu/OMSA20 leads to CNA conversions of ~48 and ~80 % after 6 h of reaction, respectively. Interestingly, the activity obtained for the Cu/OMA catalyst is essentially similar to that of Cu/OMSA20, despite the fact that Cu/OMA presents a slightly higher Cu NPs dispersion. The dominant product at low conversion is that resulting from the hydrogenation of the C=O bond, namely CNOL (selectivity of 57 % for Cu/OMSA10, 62% for Cu/OMSA20, 78 % for Cu/OMA). Selectivities to CNOL are however observed to decrease with increasing conversion (41, 47 and 51 % measured for Cu/OMSA10, Cu/OMSA20 and Cu/OMA, respectively, at XCNA = 40 %), while the selectivity to HCNA, resulting from the hydrogenation of the C=C bond, increases. Once conversion has reached 40 %, the selectivity to HCNOL, the total hydrogenation product, increases while the CNOL selectivity further decreases, down to 30 % at XCNA  80 %. At any conversion below 40%, the contribution of HCNOL to the sum of the products lies in the range 8-20% with no influence of the amount of aluminium, while the selectivity to CNOL increases and that to HCNA decreases when more aluminum is incorporated.

Scheme 1. Reaction pathways for the hydrogenation of cinnamaldehyde

The production of CNOL is usually related to the existence of large metal particles. The OMS-supported catalyst (large Cu particles, 31.5 nm) exhibits similar selectivites to CNOL and HCNA, at a low level of conversion (9%). Chambers et al. reported that 18 nm Cu particles enhance CNOL selectivity. 33 Rudolf et al. showed that the size of small nanoparticles of metallic copper (2.3 -6 nm) supported on a mesoporous silica did not influence the adsorption of CAN, and that a low selectivity to CNOL (< 5 % for a conversion of 40 %) was obtained for all catalysts. 28 In our case, it is the presence of small Cu particles and of an alumina-rich support surface that both ensures a high activity and favours the hydrogenation of the C=O bond with respect to that of the C=C bond (selectivity to CNOL between 60 and 80% at conversion < 10%, between 50 and 60% at conversion = 30%). Gutierrez et al. reported selectivities to CNOL on alumina-supported Cu catalysts very close to the values obtained here on Cu/OMA and Cu/OMSA20 (51% at a conversion of 30%), but on larger copper particles (6-19 nm, compared to the values calculated here, 3-4 nm). 37 The TOF value calculated from our results on Cu/OMSA20 (6 10 -3 s -1 at 110°C) compares well with the TOF values given for their most active Cu catalysts (5-7 10 -3 s -1 at 100°C in isopropanol). Selectivities to CNOL mentioned by Marchi et al. 34 on ternary or quaternary Al-containing catalysts, on which Cu was highly dispersed, are also similar to ours (selectivity to CNOL = 47-53% at a conversion of 30%, 36-38% at a conversion of 60%). Surface Lewis acid sites brought by the introduction of aluminum species and in the close vicinity of the metallic sites are most likely involved in the adsorption of CNA molecules via polar C=O bond, 41 as also observed for other systems, like Rh supported on oxides such as TiOx, VOx, ZrOx etc. 70 , Pt on supports manifesting SMSI effects (TiO2), 71 or bimetallic catalysts containing one of the metallic elements in a cationic state. 72

Conclusions

Al-SBA-15 mesoporous silica-alumina materials with high Al content were used as supports for copper NPs. Al was introduced during the synthesis of the SBA-15 support following a two-step pH-adjustment method, and Cu was introduced by a standard incipient wetness impregnation procedure using an aqueous solution of copper nitrate. Reference samples were prepared by using ordered mesoporous silica and alumina supports. The results evidenced that a significant gain in copper NPs dispersion was achieved over the aluminum containing supports, with a dispersion enhancement from < 5 % (when copper is dispersed over silica support) up to > 25 % (when copper is dispersed over an aluminum-rich support). It has been shown that a water-stable SBA-15 support containing 20 wt% of alumina displays a dispersing effect similar to a pure mesoporous alumina support, on which impregnation with aqueous solutions cannot be carried out owing to a lack of stability when exposed to water. The stabilization of copper over aluminum-rich supports is suggested to stem from CuAl2O4 species generated during thermal treatment.

In line with the high exposed metallic area, copper-based catalysts on aluminum-rich supports (20 wt% alumina and pure alumina) showed a high activity in the hydrogenation of cinnamaldehyde in liquid phase. These catalysts exhibited a high selectivity to cinnamyl alcohol, showing a preferential hydrogenation of the C=O bond up to a conversion of around 50%. Hence, our results show that the introduction of aluminum into the SBA-15 silica has multiple consequences: stabilization of copper as highly dispersed NPs on the internal surface of the mesopores, and improvement of both the activity and selectivity of the copper catalysts by promoting dispersion and changing the adsorption mode of CNA molecules onto the catalyst surface.

Experimental Section

Materials

All chemicals required to prepare the supports and the catalysts were used as purchased: tetraethylorthosilicate (Si(OC2H5)4, TEOS, 98%, Sigma-Aldrich), nonionic triblock copolymer Pluronic P123 (poly(ethyleneoxide)-block-poly-(propyleneoxide)-blockpoly(ethyleneoxide)-block, EO20PO70EO20, molecular weight = 5800, BASF Corp.), distilled water, hydrochloric acid (HCl, 37%, Sigma-Aldrich), nitric acid (HNO3, 67%, Sigma-Aldrich), aluminum iso-propoxide (C9H21AlO3, 98%, Sigma-Aldrich), ethanol (C2H6O, 100%, Sigma-Aldrich), aluminum nitrate (Al(NO3)2•9H2O, 98.5%, Sigma-Aldrich), ammonium hydroxide solution (28% in water, Fluka), and copper nitrate (Cu(NO3)2•3H2O, 98%, Sigma-Aldrich).

For the hydrogenation of trans-cinnamaldehyde, the chemicals were also used as purchased: trans-cinnamaldehyde (C6H5C3H3O, 98%, Merck) as reagent and iso-propanol (C3H8O, 99.9%, Sigma-Aldrich) as solvent.

Synthesis of samples Supports

OMS (ordered mesoporous silica) support was synthesized according to the procedure proposed by Zhao et al. 73 4 g of Pluronic P123 was dissolved in a 1.6 M solution of HCl at 40 C. The appropriate amount of TEOS (8.5 g) was then added dropwise to the above solution, which was subjected to magnetic stirring for 24 h. The resulting gel was submitted to hydrothermal treatment for 48 h at 100 C. After recovering by filtration, washing, and drying, the structure directing agent was removed from the pores by calcination under stagnant air at 550 C for 6 h (in a muffle furnace, using a heating ramp of 1.5 C min -1 ). OMSA (ordered mesoporous silica-alumina) supports were prepared by a two-step pH-adjustment method, according to a procedure previously established. 43 The synthesis was carried out as follows: (step 1) was performed as for OMS support, with aluminum source added after 4 h of stirring, during the aging step of 24 h. The amounts of Al(NO3)3•9H2O were calculated to obtain an alumina content of 10 and 20 wt.%, in the final materials. The aged gel was subjected to a first hydrothermal treatment at 100 °C for 72 h; for (step 2), the suspension resulting from (step 1) was cooled down to ambient temperature and then, the pH value of the mother liquor was adjusted to 7.5 by using a 4 M ammonia solution. The pH was continuously monitored with a calibrated pH-meter. After pH adjustment, the mixture was subjected to a second hydrothermal treatment at 100 °C for 48 h. The final solids were recovered by filtration, washing with distilled water, and drying at 60 °C for 12 h. The open porous structures of OMSA were obtained by calcination (see OMS preparation section). The calcined aluminosilicate samples were designated as OMSAX, where X stands for the aluminum content in the material expressed as wt.% Al2O3 (X = 10 and 20). OMA (ordered mesoporous alumina) support was prepared using the EISA method proposed by Yuan et al. 74 1 g of Pluronic P123 was dissolved at room temperature in 40 mL of anhydrous ethanol. 3 mL of HNO3 (67%) and 2.04 g of aluminum iso-propoxide were added into the mixture under vigorous stirring at room temperature for 5 h. Solvent evaporation was performed at 60 C. When the total evaporation was achieved, the yellow solid was recovered and calcined in stagnant air at 500C (1 C min -1 ).

Catalysts

Catalysts were denoted by Cu/MS where MS stands for the mesoporous support (MS = OMS, OMSAX). All nominal loadings were 5 wt.% metallic Cu. Metal oxide-loaded OMS and OMSA samples were prepared by incipient wetness impregnation using aqueous copper nitrate solutions, followed by mild drying (25 ± 1°C under air for 5 days). 25 The powders were submitted to calcination under stagnant air at 500 °C for 6 h (heating ramp of 1.5 °C min -1 ) to obtain the oxidic form of Cu. After calcination, the solids were stored under ambient conditions in sealed flasks. Cu/OMA catalyst was prepared by a solvent-free, melt infiltration (MI) procedure. 27,75,76 This method was chosen because, when contacted with water (the solvent usually used for impregnation), OMA is subjected to a significant loss of pore structure ordering. [START_REF] Bejenaru | [END_REF] The MI method involved mixing the appropriate amounts of Cu(NO3)2•3H2O and OMA, followed by grinding for 45 min. The solid mixture was heat-treated at 120 o C in sealed conditions. The obtained infiltrate was calcined under the same conditions as the other catalysts. It can be noted that MI failed to introduce copper nitrate into the porosity of siliceous SBA-15, resulting in the formation of large metallic particles outside the pores.

Physico-chemical characterization

Inductively coupled plasma optical emission spectrometry (ICP-OES) was performed on a Perkin sequential scanning spectrometer to determine the amount of copper in the catalysts. Before analysis, a known amount of calcined sample was introduced in a diluted HF-HCl solution and then digested under microwave. Powder X-ray diffraction (XRD) was performed on a Bruker AXS D5005 X-ray diffractometer, using Cu Kα radiation (λ = 1.54184 Å).

For high-angle analysis, the data were collected in the 2θ range 10 -80° with a step of 0.05° (step time of 8 s). Phase identification was made by comparison with the ICDD database. The average crystallite size of CuO was calculated according to the Scherrer equation for the most intense diffraction lines at 35.5 and 38.7°.

Nitrogen physisorption was carried out on an Autosorb MP-1 automated gas sorption system from Quantachrome Instruments. The textural data were calculated from the corresponding isotherms using conventional algorithms (B.E.T. for surface area measurement, and NL-DFT for pore size distribution). X-ray photoelectron spectroscopy (XPS) was performed on an ESCALAB 250 (Thermo Electron, Thermo Fisher Scientific, WI) photoelectron spectrometer equipped with a monochromatic Al-K source (1486.6 eV, powered at 20 mA and 10 kV). Binding energies were determined relatively to the Si 2p XPS peak set at 103.2 eV.

Transmission electronic microscopy (TEM) coupled with energydispersive X-ray spectroscopy (EDXS) was used to characterize the pore structure of the supports and the distribution of Cu. The micrographs were obtained on a JEOL 2100 instrument (operated at 200 kV with a LaB6 source and equipped with a Gatan Ultra scan camera). EDXS was carried out with a Hypernine (Premium) detector (active area: 30 mm 2 ) using the software SM-JED 2300T for data acquisition and treatment. Before analysis, the sample was included in a resin and crosscuts of ∼100 nm were realized by ultramicrotomy.

Temperature programmed reduction (TPR) and N2O chemisorption were performed on a Chembet Pulsar TPR/TPD instrument from Quantachrome. About 30 mg of calcined samples were inserted in a U-shape microreactor. Before each TPR run, the catalyst was activated at 500 °C for 1 h under a flow of simulated air (40 mL min - 1 ). After cooling down to 50 °C, the H2 containing flow was stabilized (40 mL min -1 , 5 vol.% H2 in Ar) and the TPR was performed from 50 to 700 °C, with a temperature ramp of 5 °C min -1 . The dispersion of metallic copper (DCu), average copper particle size (dCu) and active surface area (SCu) were determined using nitrous oxide chemisorption with the Quantachrome device. 58-60 30 mg of samples were introduced in the microreactor. The first reduction was carried out under the same conditions as used for TPR, up to 500 °C. The gas was switched to Ar. The reactor was cooled down to 50 °C, kept at this temperature for 30 min, and then the gas was switched to N2O (40 mL min -1 ) for 30 min in order to oxidize the zero-valent Cu surface atoms to Cu2O. The sample was again flushed with pure Ar for 30 min, and then a second TPR run was performed.

Catalytic testing

Prior to the catalytic tests, the calcined materials were crushed and sieved to select a granulometric fraction smaller than 0.126 mm. Catalysts were reduced under hydrogen flow (1 L h -1 ) at 500 °C for 10 h (heating rate 6 °C min -1 ) in order to obtain the metallic catalysts. Cinnamaldehyde hydrogenation was performed in liquid phase in a 4593 Parr reactor. The experimental conditions were: reaction temperature 110 °C, 1 mL of trans-cinnamaldehyde, 40 mL of iso-propanol, 250 mg of catalyst, H2 pressure 10 bar, stirring rate of 750 rpm. Aliquots of reaction mixture were periodically taken for analysis. Analysis was performed by GC (HP 5890 equipped with a DB-5 capillary column and a FID detector). Identification of the reactants and products was achieved from the retention times of pure compounds, and occasionally by GC-MS (Agilent 6890N system equipped with an Agilent 5973 MSD detector and a DB-5ms column). For the quantitative analyses, the FID response factors for each compound were taken into consideration.

Entry for the

Figure 1 .

 1 Figure 1. XRD patterns for Cu/OMS (a), Cu/OMSA10 (b), Cu/OMSA20 (c), and Cu/OMA (d)

[

  a] DXRD = CuO crystallite size evaluated by Scherrer equation; [b] SBET = specific surface area obtained using the BET equation; [c] Vp = total pore volume at P/P0 = 0.97; [d] Dp = mean diameter of pores evaluated by NL-DFT for cylindrical pores/equilibrium model. Wide-angle XRD (WA-XRD) patterns of the calcined Cucontaining samples are presented in Figure 1. All the samples based on OMS and OMSA display a broad peak at 2 ~24° characteristic of amorphous silica. No diffraction peaks belonging to alumina can be observed in the diffractograms of OMSA-based systems, indicating a dispersed or amorphous state for aluminic species. The diffractogram of Cu/OMA exhibits broad features assigned to poorly crystallized alumina.

Figure 2 .

 2 Figure 2. N2 physisorption isotherms (left) and NL-DFT pore size distribution (right) for Cu/OMS (a), Cu/OMSA10 (b), Cu/OMSA20 (c), and Cu/OMA (d)

Figure 3 .

 3 Figure 3. TEM images for calcined materials: (A) Cu/OMSA10, (B) Cu/OMSA20, and (C) Cu/OMA.

Figure 4 .

 4 Figure 4. Cu 2p XPS spectra of calcined (a) and reduced (b) Cu/OMSA20

Figure 5 .

 5 Figure 5. TPR profiles for Cu/OMS (a), Cu/OMSA10 (b), Cu/OMSA20 (c), and Cu/OMA (d)

Figure 6 .

 6 Figure 6. Catalytic activity (A) and selectivity to CNOL (B), HCNA (C) and HCNOL (D) for Cu/OMS (a), Cu/OMSA10 (b), Cu/OMSA20 (c), and Cu/OMA (d). Reaction conditions: Treduction = 500°C; Treaction = 110°C; 250 mg of catalyst; 1 mL of CNA, 40 mL of iso-propanol; PH2 = 10 bar; stirring rate = 750 rpm.

Table 1 .

 1 Physico-chemical properties for supports and copper containing materials

	Sample	ICP		XRD		N2 physisorption
		Al2O3, wt.%	Cu, wt.%	DXRD [a] , nm	SBET [b] , m 2 g -1	Vp [c] , cm	3 g -1	Dp [d] , nm
	OMS	-	-		810	1.12		8.4
	Cu/OMS	-	4.6	31.5	725	1.07		8.2
	OMSA10	9.1	-		351	0.86		9.2
	Cu/OMSA10		4.4	17.0	323	0.79		9.1
	OMSA20	20.5	-		360	0.76		8.4
	Cu/OMSA20		4.4	n.d.	313	0.63		8.4
	OMA	100	-		195	0.37		7.5
	Cu/OMA	-	4.9	n.d.	197	0.37		8.2

Table 2 .

 2 N2O chemisorption data for the copper-based catalysts

	Sample	DCu [a] , %	dCu [b] , nm	SCu	[c] , m	2 Cugcat -1	SCu [c] , m	CugCu 2 -1
	Cu/OMS	4.6	21.6		1.4	31.1
	Cu/OMSA10	18.8	5.3		5.6	127.2
	Cu/OMSA20	25.6	3.9		7.6	173.2
	Cu/OMA	28.0	3.6		9.3	189.4

[a] Cu dispersion; [b] Average particle size of Cu particles; [c] Accessible surface of metallic Cu

-Page No. Enhancement of copper dispersion and catalytic performances in the hydrogenation of cinnamaldehyde by incorporation of aluminium into mesoporous SBA-15

  Table of Contents (Please choose one layout)

	Layout 1:
	FULL PAPER
	Thanks to Al incorporation, Cu is
	stabilized within the pores of Al-SBA-
	15 ordered mesoporous system as
	highly dispersed nanoparticles,
	showing high catalytic activity and
	chemoselectivity towards cinnamyl
	alcohol (CNOL) in the liquid-phase
	hydrogenation of cinnamaldehyde
	(CNA)

A. Ungureanu,* A, Chirieac, C. Ciotonea, I. Mazilu, C. Catrinescu, S. Petit, E. Marceau,* S. Royer, E. Dumitriu Page No.

Acknowledgements

This work was supported by three grants of the Romanian National Authority for Scientific Research, CNCS-UEFISCDI (Projects PN-II-RU-TE-2012-3-0403, PN-II-ID-PCE-2011-3-0868, and PN-II-CT-RO-FR-2012-1-0052 -Bilateral program Hubert Curien-Brancusi). C Ciotonea (784016L) and I. Mazilu (812773E) acknowledge the Ministère des affaires étrangères et du développement international for Eiffel Excellence fellowships. C. Ciotonea acknowledges the Fondation de l'Université de Poitiers for financial support. C. Ciotonea, S. Royer and E. Marceau acknowledge Chevreul Institute (FR 2638), Ministère de l'Enseignement Supérieur et de la Recherche, Région Nord -Pas de Calais and FEDER for funding.


Xin He

Yan Wen

Chenhui Zhang

Peng Li

Dongxing Zheng

Aitian Chen

Aurelien Manchon

Xixiang Zhang

Spin transport in multilayer graphene away from the charge neutrality point

Keywords: Multilayer graphene, Ionic liquid gating, Spin transport, Elliott-Yafet mechanism

Graphene is considered as a promising material in spintronics due to its long spin relaxation time and long spin relaxation length. However, its spin transport properties have been studied at low carrier density only, beyond which much is still unknown. In this study, we explore the spin transport and spin precession properties in multilayer graphene at high carrier density using ionic liquid gating. We find that the spin relaxation time is directly proportional to the momentum relaxation time, indicating that the Elliott-Yafet mechanism still dominates the spin relaxation in multilayer graphene away from the charge neutrality point.

Introduction

Due to the high charge-carrier mobility [START_REF] Dean | Boron nitride substrates for highquality graphene electronics[END_REF], weak spin-orbit coupling [START_REF] Gmitra | Band-structure topologies of graphene: Spin-orbit coupling effects from first principles[END_REF][START_REF] Min | Intrinsic and Rashba spin-orbit interactions in graphene sheets[END_REF], and efficiently gate-tunable transport properties [START_REF] Schwierz | Graphene transistors[END_REF], graphene has become a prospective material for spintronics [START_REF] Kamalakar | Long distance spin communication in chemical vapour deposited graphene[END_REF][START_REF] Roche | Graphene spintronics: the European Flagship perspective[END_REF][START_REF] Khokhriakov | Two-dimensional spintronic circuit architectures on large scale graphene[END_REF]. When the study of spintronics in graphene just rose, the spin-relaxation time and spin-relaxation length were measured to be 170 ps and 2 µm [START_REF] Tombros | Electronic spin transport and spin precession in single graphene layers at room temperature[END_REF], respectively, which are much higher than those in the traditional materials [START_REF] Sinova | Spin Hall effects[END_REF], yet much lower than the theoretical predictions [START_REF] Huertas-Hernando | Spin-orbit-mediated spin relaxation in graphene[END_REF][START_REF] Pesin | Spintronics and pseudospintronics in graphene and topological insulators[END_REF]. The poor quality of the tunneling layer [START_REF] Han | Tunneling spin injection into single layer graphene[END_REF], the polymer residues during the fabrication process [START_REF] Avsar | Electronic spin transport in dual-gated bilayer graphene[END_REF], and the presence of scattering centers in the substrate [START_REF] Guimaraes | Spin transport in high-quality suspended graphene devices[END_REF][START_REF] Zomer | Long-distance spin transport in highmobility graphene on hexagonal boron nitride[END_REF] are considered to be the main factors that degrade the spin transport properties of graphene.

Correspondingly, several measures have been proposed to address these problems, such as the optimized insulating layer which facilitates the tunneling injection [START_REF] Han | Tunneling spin injection into single layer graphene[END_REF][START_REF] Kamalakar | Enhanced tunnel spin injection into graphene using chemical vapor deposited hexagonal boron nitride[END_REF][START_REF] Singh | Nanosecond spin relaxation times in single layer graphene spin valves with hexagonal boron nitride tunnel barriers[END_REF][START_REF] Kamalakar | Inversion of Spin Signal and Spin Filtering in Ferromagnet vertical bar Hexagonal Boron Nitride-Graphene van der Waals Heterostructures[END_REF], the hexagonal boron nitride (hBN) layer which serves as a substrate or encapsulating layer to provide a smooth interface with reduced trapped charges [START_REF] Gurram | Electrical spin injection, transport, and detection in graphenehexagonal boron nitride van der Waals heterostructures: progress and perspectives[END_REF], and the suspended sample structure to avoid impurity scattering [START_REF] Guimaraes | Spin transport in high-quality suspended graphene devices[END_REF][START_REF] Drogeler | Spin lifetimes exceeding 12 ns in graphene nonlocal spin valve devices[END_REF][START_REF] Drogeler | Nanosecond spin lifetimes in single-and few-layer graphene-hBN heterostructures at room temperature[END_REF]. Thanks to these efforts, the spin-relaxation time and spinrelaxation length of graphene have been improved up to 12.6 ns and 30.5 µm [START_REF] Drogeler | Spin lifetimes exceeding 12 ns in graphene nonlocal spin valve devices[END_REF], respectively, approaching the theoretical prediction gradually [START_REF] Huertas-Hernando | Spin-orbit-mediated spin relaxation in graphene[END_REF][START_REF] Pesin | Spintronics and pseudospintronics in graphene and topological insulators[END_REF]. Alternatively, multilayer graphene was also proposed as a promising strategy to improve the spin transport performance since the outer layers of multilayer graphene can screen the electrical potential of impurities [START_REF] Li | Spin diffusion and non-local spin-valve effect in an exfoliated multilayer graphene with a Co electrode[END_REF][START_REF] Maassen | Comparison between charge and spin transport in few-layer graphene[END_REF]. However, this original study was limited to relatively low carrier density (close to 1.25×10 !" 𝑐𝑚 !! ) because of the low efficiency of the traditional solid gate [START_REF] Maassen | Comparison between charge and spin transport in few-layer graphene[END_REF]. Hence, the conclusion that the spin relaxation is mainly due to the Elliott-Yafet mechanism was drawn only in the vicinity of the charge neutrality point [START_REF] Maassen | Comparison between charge and spin transport in few-layer graphene[END_REF]. When the Fermi level is shifted away from the charge neutrality point, the energy band structure and thus the electron properties will be quite different.

Furthermore, the electron-electron interaction will be weakened due to the high-carrier density induced charge screening, which will influence the momentum relaxation time. Therefore, how spin transport evolves and which mechanism dominates the spin relaxation in the high-carrierdensity range of multilayer graphene is yet to be explored. Ionic liquid (IL) gating has been widely used in 2D materials due to its high efficiency in tuning carrier density [START_REF] Du | Electric double-layer transistors: a review of recent progress[END_REF][START_REF] Fujimoto | Electric-double-layer field-effect transistors with ionic liquids[END_REF]. It has been central to the exploration of the electronic and optical properties [START_REF] Ye | Accessing the transport properties of graphene and its multilayers at high carrier density[END_REF][START_REF] Braga | Quantitative determination of the band gap of WS2 with ambipolar ionic liquid-gated transistors[END_REF][START_REF] He | Formation of p-n-p junction with ionic liquid gate in graphene[END_REF][START_REF] He | Exploration of exciton behavior in atomically thin WS2 layers by ionic gating[END_REF], for the tuning and realization of voltage-controlled superconductivity [START_REF] Saito | Gate-induced superconductivity in two-dimensional atomic crystals[END_REF], and manipulation of the magnetic properties [START_REF] Deng | Gate-tunable room-temperature ferromagnetism in two-dimensional Fe3GeTe2[END_REF]. These studies all prove that IL gating is a mature and reliable technique.

In the present study, we explore the spin transport in multilayer graphene by tuning its carrier density using IL gating at 200 K. In a heavily doped multilayer graphene sample, the electron and hole densities reach as high as 4.3×10 !" 𝑐𝑚 !! and 5.7×10 !" 𝑐𝑚 !! , respectively.

The spin relaxation time and the momentum relaxation time are found to follow a directly proportional relationship, which indicates that the Elliott-Yafet mechanism, rather than the D'yakonov-Perel' mechanism, plays an important role in the high-carrier density range.

Results and analysis

The schematic and the optical images of the device for the spin and charge transport measurements are shown in Fig. 1(a-c), while the image with IL on the device can be found in Supplementary Fig. S1. All the electrodes are made of Co/Au (60 nm/ 20 nm). As presented in Fig. 1(a), the outer four electrodes [START_REF] Dean | Boron nitride substrates for highquality graphene electronics[END_REF][START_REF] Gmitra | Band-structure topologies of graphene: Spin-orbit coupling effects from first principles[END_REF][START_REF] Kamalakar | Long distance spin communication in chemical vapour deposited graphene[END_REF][START_REF] Roche | Graphene spintronics: the European Flagship perspective[END_REF] are designed for charge transport measurement using a four-probe configuration, while the inner four electrodes [START_REF] Gmitra | Band-structure topologies of graphene: Spin-orbit coupling effects from first principles[END_REF][START_REF] Min | Intrinsic and Rashba spin-orbit interactions in graphene sheets[END_REF][START_REF] Schwierz | Graphene transistors[END_REF][START_REF] Kamalakar | Long distance spin communication in chemical vapour deposited graphene[END_REF] are designed for spin transport measurement; the widths of the central two electrodes [START_REF] Min | Intrinsic and Rashba spin-orbit interactions in graphene sheets[END_REF][START_REF] Schwierz | Graphene transistors[END_REF] are fabricated to be 270 and 660 nm, respectively, forming injector and detector electrodes with different coercive fields.

The thickness of the multilayer graphene flake was measured to be 4 nm by atomic force microscope. Considering the spacing of adjacent graphene layers is 0.335 nm [START_REF] Maassen | Comparison between charge and spin transport in few-layer graphene[END_REF], the layer number of our sample can be estimated to be 12, with an error of about 1 layer. Note that all the measurements in this report were performed at 200 K unless otherwise specified. After dropping ionic liquid to cover the whole sample and setting the IL gate floating, the measurements for both spin transport and spin precession were first conducted at 300 K. Subsequently, the system was cooled down quickly to 200 K, the freezing point of ionic liquid [START_REF] Perera | Improved carrier mobility in few-layer MoS2 field-effect transistors with ionic-liquid gating[END_REF], so as to avoid chemical reaction between graphene and ionic liquid. The same experiments were also performed at 200 K. As shown in Fig. 1(a), under an in-plane magnetic field (B y ), a current I is injected from electrode 3 (injector) to 2, and then spins accumulate underneath electrode 3 and diffuse along the graphene channel to both sides, inducing a voltage drop V NL between electrode 4 (detector) and 5 and defining the nonlocal resistance R NL = V NL /I (spin transport) [START_REF] Tombros | Electronic spin transport and spin precession in single graphene layers at room temperature[END_REF]. The detector probes the chemical potential of the spins impinging from the injector. Therefore, switching the magnetization configuration of the injector and detector from parallel (P) to antiparallel (AP) by sweeping the in-plane magnetic field leads to a sudden decrease of both the nonlocal voltage (ΔV NL ) and the nonlocal resistance (ΔR NL = ΔV NL /I), as shown in Fig. 2(a) [START_REF] Tombros | Electronic spin transport and spin precession in single graphene layers at room temperature[END_REF][START_REF] Jedema | Electrical spin injection and accumulation at room temperature in an all-metal mesoscopic spin valve[END_REF]. The non-local spin valve signal obtained at 200 K is enhanced in comparison to that measured at 300K (Fig. 2(a)), due to the decreased Coulomb scattering in graphene with relatively low carrier mobility, which is consistent with previous reports [START_REF] Maassen | Long spin relaxation times in wafer scale epitaxial graphene on SiC[END_REF][START_REF] Yamaguchi | Spin relaxation in weak localization regime in multilayer graphene spin valves[END_REF].

Applying an out-of-plane magnetic field (B z ), a current I is also injected from electrode 3 to 2, and spins precess as they diffuse from electrode 3 to 4, inducing a voltage drop V NL between electrode 4 and 5 and defining the nonlocal resistance R NL = V NL /I as well (spin precession). The non-local Hanle spin precession signals in both 300 K and 200 K are plotted in Supplementary Fig. S2(a) and (b), respectively. At low magnetic fields, the signal weakens because of the reduction of the spin polarization upon increasing the field from zero. When the signal corresponding to the P (AP) configuration decreases (increases) to zero, it indicates that the average angle of the precession reaches 90° [START_REF] Jedema | Electrical detection of spin precession in a metallic mesoscopic spin valve[END_REF]. Upon increasing the field further, R NL changes sign and reaches a minimum (maximum) gradually, indicating the reversal of spins (the average precession angle reaches 180°) [START_REF] Jedema | Electrical detection of spin precession in a metallic mesoscopic spin valve[END_REF]. However, even though the vertical magnetic field was applied up to 6 T, we could not observe the saturation of R NL that indicates the magnetization of the electrodes has completely rotated out of the plane. This is because the large-resistance background blurs the spin precession signal [START_REF] Guimaraes | Controlling spin relaxation in hexagonal BN-encapsulated graphene with a transverse electric field[END_REF]. Thus, we can only estimate that R NL saturates at around 1.5 T, at which the slopes of the curves decrease suddenly (Supplementary Fig. S2). To verify the validity of IL gating at 200 K, we performed electronic measurements and then obtained the standard transfer curve of graphene, with the charge neutrality point located at V CNP = 0.59 V (Fig. S3). This means that IL gating works even at 200 K, suggesting that the molecules of IL are still mobile at 200 K. Moreover, from Fig. S3, the mobilities (µ) of both electron and hole are extracted to be 242 cm ! /𝑉 • 𝑠 and 197 cm ! /𝑉 • 𝑠, respectively, by using the

equation, µ = Δσ/(C !" ΔV ! ) = (Δ( ! !" ! ! ) • ! ! )/(C !" ΔV ! )
, where σ is the conductivity, I ds is the drain-source current, V d is the drain-source voltage, V g is the gate voltage, C IL is the capacitance per unit area of IL, L and W are the channel length and width, respectively [START_REF] Kim | Large-scale graphene micropatterns via self-assembly-mediated process for flexible device application[END_REF]. C IL is obtained to be 2.86 µF/cm 2 and 2.55 µF/cm 2 for the electron and hole, respectively (Supplementary Note 1), much lower than the values reported previously [START_REF] Ye | Accessing the transport properties of graphene and its multilayers at high carrier density[END_REF][START_REF] Zhang | Ambipolar MoS2 thin flake transistors[END_REF]. There are two factors leading to the small C IL : the extra TiO X layer and the screening effect in graphene. First, the inserted TiO X layer (about 1 nm thick) increases the effective dielectric thickness of the geometrical capacitance C g , which is defined by C ! = ε !" ε ! /d, where ε IL is the relative permittivity of ionic liquid, ε 0 is the vacuum permittivity, and d is the dielectric thickness [START_REF] Uesugi | Electric double-layer capacitance between an ionic liquid and few-layer graphene[END_REF]. Second, the charges distributed in the first 3 or 4 layers screen the electric field and quench its penetration, and then C g saturates to a constant as the layer number is larger than 4 [START_REF] Uesugi | Electric double-layer capacitance between an ionic liquid and few-layer graphene[END_REF]. In Fig. 2(b), the IL gate voltage is transformed to the carrier concentration N by using the formula, N = C !"

(! ! !! !"# ) ! ,
where e<0 is the elementary charge [START_REF] Perera | Improved carrier mobility in few-layer MoS2 field-effect transistors with ionic-liquid gating[END_REF][START_REF] Misra | Electric field gating with ionic liquids[END_REF]. Then the highest electron (hole) density can be calculated to be 4.3×10 !" 𝑐𝑚 !! (5.7×10 !" 𝑐𝑚 !! ). Also in Fig. 2(b), the changes of the non-local resistances were extracted by performing the non-local spin transport at a series of carrier densities, as a comparison to the carrier density dependence of the conductivity. The relationship between ΔR NL and σ in Fig. 2(b) follows neither ∆𝑅 !" ~𝜎 nor ∆𝑅 !" ~! ! , implying that the electrode contact could be a pinhole type [START_REF] Han | Tunneling spin injection into single layer graphene[END_REF]. As discussed in the previous study, the contact interface mainly affects the polarization of the injected current [START_REF] Maassen | Comparison between charge and spin transport in few-layer graphene[END_REF], and the pinhole-type contacts could also create long spin relaxation time and spin relaxation length [START_REF] Drogeler | Spin lifetimes exceeding 12 ns in graphene nonlocal spin valve devices[END_REF], thus the contacts barely influence the spin diffusion and spin relaxation in our study. The spin precession at different carrier densities is studied applying a perpendicular magnetic field at 200 K, as shown in Fig. 3. As an example, Fig. 3(a) shows the nonlocal Hanle spin precession signals with both P and AP magnetization of the central electrodes at V g = 2 V, corresponding to the electron density of 2.5×10 !" 𝑐𝑚 !! . Since the magnetic field is applied perpendicular to the plane, the spins with in-plane orientation begin to precess, following the Bloch equation for the spin accumulation µ S , [START_REF] Maassen | Comparison between charge and spin transport in few-layer graphene[END_REF]. The first, second and third term describes the spin diffusion, spin relaxation, and spin precession, respectively, with D S is the spin diffusion coefficient, τ S is the spin relaxation time, and ω ! = !! ! 𝐇 ħ is the Larmor frequency (g = 2 is the effective Lande factor, µ B is the Bohr magneton, and H is the outof-plane magnetic field). Moreover, D S and τ S can be extracted by fitting the Hanle precession signal at each carrier density with the solution of the above Bloch

D ! ∇ ! µ ! - ! ! ! ! + ω ! ×µ ! = 0
equation,R !" !(!") = (±)S ! ! !"! ! ! exp (- ! ! !! ! ! )cos (ω ! t)exp (- ! ! !
)dt [START_REF] Pi | Manipulation of spin transport in graphene by surface chemical doping[END_REF] (where 𝑅 !" !(!") is the nonlocal magnetoresistance, +(-) sign is for the P (AP) state, S 0 is the amplitude, and L is the channel length between the central two electrodes), while the spin diffusion length λ S can be calculated by λ ! = D ! τ ! [START_REF] Maassen | Comparison between charge and spin transport in few-layer graphene[END_REF]. Note that before fitting the result, the P curve of the Hanle effect was subtracted from the AP curve (e.g., Fig. 3(a)) to remove the background that is not related to the spin signal [START_REF] Maassen | Comparison between charge and spin transport in few-layer graphene[END_REF]. Therefore, the relationships between σ and n, D S and n, τ S and n, and λ S and n are presented in Fig. 3(b-e), respectively. D S , τ S, and λ S all diminish with decreasing the electron or the hole density, and reach the minimum values in the vicinity of the charge neutrality point, resembling the relationship between σ and n. Compared with the previous data obtained using the traditional solid gate in multilayer graphene [START_REF] Maassen | Comparison between charge and spin transport in few-layer graphene[END_REF], these parameters show more significant changes because of the stronger carrier density modulation induced by the IL gating. The asymmetry of D S for electrons and holes in Fig. 3(c) may be due to the nonuniformities in the carrier density [START_REF] Jozsa | Linear scaling between momentum and spin scattering in graphene[END_REF], while the asymmetry of τ S in Fig. 3(d) could be traced to the different impacts of the spin scatterings to the electrons and holes. In semiconductors, there are two main mechanisms that contribute to spin relaxation. The first one is the Elliott-Yafet mechanism, which stems from momentum transfer between spin and orbital angular momenta during momentum scattering events, resulting in τ S directly proportional to τ P (τ P is the momentum relaxation time) [START_REF] Jozsa | Linear scaling between momentum and spin scattering in graphene[END_REF][START_REF] Avsar | Toward wafer scale fabrication of graphene based spin valve devices[END_REF]. The second one is the D'yakonov-Perel' mechanism, which stems from scattering-induced destructive interference of spins coherently precessing around the spin-orbit field in the band structure. In this case, τ S is inversely proportional to τ P [START_REF] Jozsa | Linear scaling between momentum and spin scattering in graphene[END_REF][START_REF] Avsar | Toward wafer scale fabrication of graphene based spin valve devices[END_REF]. Thus, to identify which spin scattering mechanism dominates the spin relaxation, we need to determine the relationship between τ S and τ P . These relationships are consistent with the previous studies in both monolayer and multilayer graphene [START_REF] Maassen | Comparison between charge and spin transport in few-layer graphene[END_REF][START_REF] Jozsa | Linear scaling between momentum and spin scattering in graphene[END_REF][START_REF] Han | Spin relaxation in single-layer and bilayer graphene[END_REF]. Combing the facts that D ! ∝ D ! (D P is the momentum diffusion coefficient) [START_REF] Maassen | Comparison between charge and spin transport in few-layer graphene[END_REF] and D ! ∝ τ ! [START_REF] Pi | Manipulation of spin transport in graphene by surface chemical doping[END_REF][START_REF] Jozsa | Linear scaling between momentum and spin scattering in graphene[END_REF], we obtain τ ! ∝ τ ! . This indicates that the Elliott-Yafet mechanism still plays a more important role in the spin relaxation process of multilayer graphene at high carrier density [START_REF] Jozsa | Linear scaling between momentum and spin scattering in graphene[END_REF][START_REF] Avsar | Toward wafer scale fabrication of graphene based spin valve devices[END_REF]. Moreover, a large perpendicular electric field can lead to a considerable in-plane Rashba spin-orbit field, which will affect the spin transport properties, i.g., the spin relaxation time [START_REF] Guimaraes | Controlling spin relaxation in hexagonal BN-encapsulated graphene with a transverse electric field[END_REF]. In our experiment, the gate voltage up to 3 V is applied, and the electric field penetrates to the first 3 or 4 layers of graphene. Considering the ionic liquid molecule (diameter ̴ 1 nm) and the TiO X layer (thickness ̴ 1 nm), the electric field reaches around 1 V/nm in graphene theoretically. However, we do not see the deviation of the spin relaxation time in Fig. 3(d), which means that the induced spin-orbit field is quite weak, because of the strong screening of the charges to the applied electric field.

Outlook

Recently, to realize more functions in graphene-based spin devices, the heterostructures such as graphene/topological insulators [START_REF] Khokhriakov | Tailoring emergent spin phenomena in Dirac material heterostructures[END_REF], graphene/magnetic insulator [START_REF] Karpiak | Magnetic proximity in a van der Waals heterostructure of magnetic insulator and graphene[END_REF], and graphene/Transition metal dichalcogenides [START_REF] Dankert | Electrical gate control of spin current in van der Waals heterostructures at room temperature[END_REF] have been constructed. These heterostructures combine both the excellent spin properties of graphene and the exotic properties of other materials via proximity effect. Moreover, the nano-graphene based spin device is also drawing more attention [START_REF] Fan | Giant decreasing of spin current in a single molecular junction with twisted zigzag graphene nanoribbon electrodes[END_REF][START_REF] Fan | Redox control of magnetic transport properties of a single anthraquinone molecule with different contacted geometries[END_REF][START_REF] Fan | Symmetry-dependent spin transport properties of a single phenalenyl or pyrene molecular device[END_REF][START_REF] Enoki | Magnetism of Nanographene-Based Microporous Carbon and Its Applications: Interplay of Edge Geometry and Chemistry Details in the Edge State[END_REF]. Due to the weak tuning capability of the traditional solid gate, the gate voltage that is used to switch on/off the spin current is always as high as several tens of volts. We believe that the IL gate can be used to improve the gating effect in these devices and allow them to work practically at low voltages.

Conclusion

In conclusion, we provide more insights into the spin dynamics in multilayer graphene and verified that IL gating is a powerful technique investigating spintronics in 2D materials.

After applying IL gate, the electron and hole densities in multilayer graphene are increased up to 4.3×10 !" 𝑐𝑚 !! and 5.7×10 !" 𝑐𝑚 !! , respectively. The spin relaxation time and the momentum relaxation time are found to follow a directly proportional relationship, indicating that the Elliott-Yafet mechanism still dominates the spin relaxation in multilayer graphene at high carrier density.

Methods

Multilayer graphene samples were exfoliated from natural graphite onto SiO 2 /Si (300 nm / 500 µm) substrate. To avoid the extra etching process, the samples which were relatively long and narrow were located by optical microscope (Carl Zeiss Imager.A2 Vario with AxioCam

Fig. 1 .

 1 Fig. 1. (a) Schematic of a nonlocal spin-valve device under IL gating. (b) Optical image of a nonlocal spin-valve device (scale bar: 50 µm). A 12-layer graphene flake is used as the channel. (c) Zoom in on the image framed by the red open box in (b); the scale bar is 10 µm.

Fig. 2 .

 2 Fig. 2. (a) Non-local spin valve signal of the 12-layer graphene flack applying an in-plane magnetic field at 200 K and 300 K. The magnetization configurations of the electrodes are illustrated by the arrows. The amplitude of the nonlocal resistance is defined as the change of the nonlocal resistance, ΔR NL . (b) σ (left) and ΔR NL (right) as functions of the carrier density n, respectively.

Fig. 3 .

 3 Fig. 3. (a) Nonlocal Hanle spin precession signal of the 12-layer graphene flake at the electron density of 2.5×10 !" 𝑐𝑚 !! , P and AP magnetizations of the central two electrodes are illustrated in the insets. σ (b), D S (c), τ S (d) and λ S (e) as a function of the carrier density n, respectively.

Fig. 4 .

 4 Fig. 4. (a) Linear relationship between τ S and D S , extracted from Fig. 3(c) and (d). (b) Linear relationship between λ S and D S , extracted from Fig. 3(c) and (e).

Fig. 4 (

 4 a) shows the linear dependence between τ S and D S (τ ! ∝ D ! ), extracted from Fig. 3(c) and (d); while in Fig. 4(b), λ Sand D S are also in direct proportion, extracted from Fig.3(c) and (e). The larger slopes for the hole regimes indicates the averagely lower spin-flip probability for holes in each scattering[START_REF] Maassen | Comparison between charge and spin transport in few-layer graphene[END_REF].

Acknowledgments

This publication is based on research supported by the King Abdullah University of Science and Technology (KAUST), under Award Nos. OSR-2018-CRG7-3717 and OSR-2017-CRG6-3427.

HRc), and then their thicknesses were confirmed by atomic force microscope (Asylum Research MFP-3D) using tapping mode. After that, 0.8 nm Ti was evaporated onto the samples using ebeam evaporator, and then they were exposed in air for 0.5 h to form an oxidation layer (the TiO X layer not only plays a role of tunneling layer to conquer the conductance mismatch problem but also prevents the sample from reacting with ionic liquid) [START_REF] Song | Anomalous Hall effect in one monolayer cobalt with electrical manipulation[END_REF]. Subsequently, six electrodes were fabricated for each sample using e-beam lithography (EBL) and followed ebeam evaporation of Co/Au (60 nm / 20 nm). After the lift-off process, the samples were immersed in chloroform (60 , 4 hours) to remove the polymer residues to facilitate the usage of

Conflicts of interest

The authors declare no competing financial interests.

Appendix A. Supplementary data

Supplementary data to this article can be found online at XXXXXXXXX.


Patrick Cassam-Chenaï

Gilles Lebeau

Smeared Coulomb potential orbitals: I-Asymptotic expansion

We consider an 1-electron model Hamiltonian, whose potential energy corresponds to the Coulomb potential of an infinite wire with charge Z distributed according to a Gaussian function. The time independent Schrödinger equation for this Hamiltonian is solved perturbationally in the asymptotic limit of small amplitude vibration (Gaussian function width close to zero). We propose to use the naturally polarized functions soobtained, as orbital basis sets for quantum chemical calculations. In particular, they should be well suited to perform electron-nucleus mean field configuration interaction calculations. Since the free-parameters of the model have the remarkable property to factorize the perturbative corrections to the eigenfunctions, these corrective part in factor can be simply added as additional functions to standard basis sets, leaving it to the molecular orbital calculation to optimize the free parameters within molecular orbital coefficients.

Introduction

This paper is dedicated to Prof. Graham Chandler for his 80 th birthday. His famous "Mclean and Chandler basis sets" have proved extremely useful to the quantum chemistry community. Some twenty years ago, we optimized Gaussian basis sets for molecular fragments together with D. Jayatilaka and G. S. Chandler 1 . However, this endeavour was suspended due to technical difficulties raised by electronegative atoms. In the present article, we come back to, perhaps, a more original approach to basis functions, where the latter are not selected on the ground of their technical advantages, as was the case initially for Gaussian-type orbitals (GTO) 2 , but because they are eigenfunctions of a model Hamiltonian and therefore have some physical relevance. It is hoped that this property can be taken advantage of in basis set truncation.

The model Hamiltonian we will consider, is a generalization of the hydrogenoid atom Coulomb Hamiltonian. So, our new orbitals will be part of the exponential-type orbitals (ETO) 3 family, like the hydrogenoid atom eigenfunctions, which constitute an asymptotic limit. Slater-type orbitals (STO) 4,5 is another type of ETO related to hydrogenoid orbitals (HO): they can be seen as "uncontracted" HO, that is to say, as HO with the Laguerre polynomial prefactor replaced by a simple monomial one. In contrast to GTO, the difficulty of computing multicenter integrals with STO has led to the introduction of more ETO family members such as Bessel-type orbitals (BTO) 6 or Coulomb-Sturmian orbitals (CSO) 7 . The techniques developed for the latter 8-10 will be equally relevant for integrals involving our deformed hydrogenoid orbitals (DHO).

It is our take that DHO will be particularly useful for the recently developed electronnucleus mean field configuration interaction (EN-MFCI) method [START_REF] Cassam-Chenaï | ??A920125022015[END_REF] . The EN-MFCI method affords one to obtain in a single calculation, the electronic and vibrational energy levels of a molecule, without making the "Born-Oppenheimer" (BO) approximation [START_REF] Cassam-Chenaï | ??A920125022015[END_REF]12 . In contrast, the traditional methods of Quantum Chemistry are set in the frame of this approximation.

They describe electronic clouds of fixed nuclear configurations and make use of orbital basis sets centered on nuclear positions. The latter basis sets are not appropriate for EN-MFCI calculations and their discrepancies have been bypassed so far, only by adding off-centered orbitals. However, the addition of such functions introduces linear dependencies within the orbital set and spans virtual molecular orbitals of little relevance for the description of low energy wave functions of the molecule. So, it appears important to develop new orbital basis functions for the EN-MFCI method, able to describe the smeared electronic cloud of oscillating nuclei in a molecular system. This article is the first part of a series aiming at deriving appropriate basis functions for EN-MFCI calculations. We propose to use the eigenfunctions of a one-electron model

Hamiltonian corresponding to a Coulomb potential convoluted with a Gaussian function.

The latter can be seen as a ground state vibrational eigenfunction for the nuclear internal motion in the harmonic approximation. In the next section, we show how such Hamiltonian can be expanded asymptotically near the infinitely small Gaussian width limit, where the hydrogenoid atom eigenfunctions are retrieved. Then, in section 3, we solve the eigenproblem for the first order corrected potential by means of the Rayleigh-Schrödinger perturbation theory. We tabulate the fist order corrected eigenfunctions and associated eigenvalues up to n = 7. In the last section, we provide details and perspectives on how these eigenfunctions can be employed in multi-electron, quantum chemical calculations. Numerical applications of EN-MFCI making use of these new functions is postponed to part two of this series.

Asymptotic expansion of a smeared Coulombic Hamiltonian

Let us consider the one-electron model Hamiltonian: H = -2µ + V ( r), with a potential of the form:

V ( r) = -Z a π +∞ -∞ exp[-az 2 0 ] r -r z 0 dz 0 , (1)
where, Z ∈ N * , a ∈ R * + , r = (x, y, z) and r z 0 = (0, 0, z 0 ) in Cartesian coordinates. When r is expressed in cylindrical coordinates, r = (ρ, φ, z), the potential depends only upon ρ and z,

V (ρ, z) = -Z a π +∞ -∞ exp[-az 2 0 ] ρ 2 + (z -z 0 ) 2 dz 0 . (2)
This potential corresponds to the Coulomb potential of an infinite wire with charge Z distributed according to a Gaussian function. In the limiting case of a Gaussian function sharply peaking at the origin (a → +∞, Dirac distribution limit), the system will tends towards a point-charge Z concentrated at the origin and the hydrogenoid atom eigenfunctions will be recovered.

However, by taking a Gaussian width parameter of the order of magnitude of a nucleus vibration amplitude, we will get basis functions corresponding to a Coulomb potential convoluted by a nuclear, vibrational, harmonic motion, that we may think particularly appropriate for EN-MFCI calculations.

Unfortunately, the Schrödinger equation for this potential is hard to solve because the whole z-axis is singular. So, we will restrict ourselves to the a → ∞ asymptotic limit, and expand the potential V (ρ, z) as

V (ρ, z) = -Z √ a π +∞ 0 dλ √ λ +∞ -∞ dz 0 exp[-az 2 0 -λ(ρ 2 + (z -z 0 ) 2 )] = -Z a π +∞ 0 dλ √ λ(a+λ) exp[-λ(ρ 2 + z 2 )]exp[ λ 2 a+λ z 2 ] = -Z √ π +∞ 0 dλ 1 √ λ + -λ 2 +z 2 λ 2 √ λa + λ 3/2 (3-12z 2 λ+4z 4 λ 2 ) 8a 2 + λ 5/2 (-15+90z 2 λ-60z 4 λ 2 +8z 6 λ 3 ) 48a 3 + o 1 a 7 2 exp[-λ(ρ 2 + z 2 )]
(3)

Swaping the limits and setting r =

ρ 2 + z 2 , ρ = r × sin(θ), z = r × cos(θ), with θ ∈ [0, π],
the potential becomes

V (r, θ) = -Z r + (1-3cos(θ) 2 )Z 4r 3 a - 3((3-30cos(θ) 2 +35cos(θ) 4 )Z) 32r 5 a 2 + 15(5-105cos(θ) 2 +315cos(θ) 4 -231cos(θ) 6 )Z 128r 7 a 3 + O 1 a 4 , (4)
where we recognize the hydrogenoid atom potential in the zero th -order term,

V (0) (r, θ) = - Z r . (5)
Now, at any order, the singularity is located at the single point r = 0. In the following, we will only consider the V (r, θ) potential truncated at first order,

V 1 (r, θ) := - Z r + (1 -3cos(θ) 2 ) Z 4r 3 a . ( 6
)
Since this potential is not bounded from below, we will not attempt to solve the eigenproblem exactly. Instead, we apply Rayleigh-Schrödinger perturbation theory, to obtain the first order corrections to the hydrogenoid atom eigenstates.

Perturbationally corrected eigenstates

So, we consider the following Hamiltonian in spherical coordinates and atomic units

H = -1 2µ 1 r 2 ∂ ∂r r 2 ∂ ∂r + 1 r 2 sin(θ) ∂ ∂θ sin(θ) ∂ ∂θ + 1 r 2 sin(θ) 2 ∂ 2 ∂φ 2 + V 1 (r, θ), (7)
defined on the Hilbert space of square integrable functions whose scalar product is expressed as,

ψ 1 |ψ 2 := +∞ 0 r 2 dr π 0 sin(θ)dθ 2π 0 dφ ψ * 1 (r, θ, φ)ψ 2 (r, θ, φ). (8)
A Galerkin-type approach similar to the one proposed in 13 , with spherical harmonics in place of Chebychev basis functions could be considered to solve its eigenvalue problem. However, it is more practical to approach the eigenstates perturbationally, starting from the well-known solutions of time-independent Schrödinger equation for the hydrogenoid atom,

ψ (0) n,l,m (r, θ, φ) = R n,l (r)Y l,m (θ, φ) (9)
with

R n,l (r) = 2µZ n 3 2 (n -l -1)! 2n[(n + l)!] exp - µZr n 2µZr n l L 2l+1 n-l-1 2µZr n , (10)
where L 2l+1 n-l-1 (x) denotes the generalized Laguerre polynomials, and Y l,m (θ, φ) the spherical harmonics. We note that the perturbation operator,

V (1) (r, θ) := 1 a (1 -3cos(θ) 2 ) Z 4r 3 . ( 11
)
is proportional to Y 2,0 (θ, φ),

V (1) (r, θ) = π 5 -Z a Y 2,0 (θ, φ) r 3 . ( 12
)
Given the following integral formula

π 0 sin(θ)dθ 2π 0 dφ Y l 1 ,m 1 (θ, φ)Y l 2 ,m 2 (θ, φ)Y l 3 ,m 3 (θ, φ) = (2l 1 + 1)(2l 2 + 1)(2l 3 + 1) 4π   l 1 l 2 l 3 0 0 0     l 1 l 2 l 3 m 1 m 2 m 3   , (13)
and the well-known relation for the conjugate of a spherical harmonic:

Y * l,m (θ, φ) = (-1) m Y l,-m (θ, φ), (14)
we deduce that, for a given (n, l, m)-triplet of quantum numbers, the state ψ

(0)
n,l,m can only be coupled at first order to states ψ Table 1: List of hydrogenoid eigenstates coupled by the perturbation operator of Eq.( 12) to a given hydrogenoid eigenstate (in terms of their associated quantum numbers).

(0) n ,l ,m 's such that (i) l ≥ |m|, (ii) for l ∈ {|l -2|, • • • , l + 2}, the 3-j symbol   l l 2 0 0 0   is
Zero-order states states possibly coupled by first-order perturbation n > 0, l = 0, m = 0

(n > 2, l = 2, m = 0) n > 1, l = 1, m ∈ {-1, 0, 1} (n > 1, l = 1, m = m); (n > 3, l = 3, m = m) n > l ≥ 2, m ∈ {-l, -l + 1, l -1, l} (n > l, l = l, m = m); (n > l + 2, l = l + 2, m = m) n > l ≥ 2, m ∈ {-l + 2, • • • , l -2} (n > l -2, l = l -2, m = m); (n > l, l = l, m = m); (n > l + 2, l = l + 2, m = m)
A priori, the perturbation operator needs to be diagonalized first in each degenerate n-subspace. The first order correction to the unperturbed energies, (that is the energies of the hydrogenoid atom, E

n,l,m = -µZ 2 2n 2 in hartree), are the eigenvalues of the matrix, ( ψ

n,l,m |V (1) |ψ

(0) n,l ,m ) (l,m),(l ,m ) = (δ m,m ψ (0)
n,l,m |V (1) |ψ

(0) n,l ,m ) (l,m),(l ,m ) . (15)
However, this matrix is already diagonal: For m = m and l = l a non-zero integration on angular variables implies l = l ± 2, as we have seen. Suppose without loss of generality, that l = l + 2, the problem being symmetrical in l and l . The generalized Laguerre polynomial L 2l+1 n-l-1 to the factor appearing in front of the generalized Laguerre polynomials in the radial integral R n,l | 1 r 3 |R n,l . So, for l = 0, there is no first order correction,

2µZr n = L 2l +5 n-l -3 2µZr
E n,0,0 = -µZ 2 2n 2 , (16)
and for l > 0, the correction is,

E n,l,m = E (0) n,l,m - π 5 Z a ψ (0) n,l,m | Y 2,0 (θ, φ) r 3 |ψ (0) n,l,m = - µZ 2 2n 2 - π 5 Z a +∞ 0 R 2 n,l (r) r dr × (-1) m (l + 1 2 ) 5 π   l l 2 0 0 0     l l 2 -m m 0   = - µZ 2 2n 2 - Z a × (-1) m (l + 1 2 )   l l 2 0 0 0     l l 2 -m m 0   +∞ 0 R 2 n,l (r) r dr. (17)
As expected, the spherical symmetry is broken: an (l, |m|)-dependency is introduced in the perturbed eigenvalues. In Tab.2, we provide the first order eigenvalues, which can be useful for basis set truncation purposes. We note that, at first order, degeneracy is not completely lifted, as for example, E 4,0,0 = E 4,3,±2 . If we factorize by the zero order energy, we see that the relative, first order corrections are all proportional to µ 2 Z 2 a .

Table 2: First-order corrected energies (up to n=4). For every pairs, (n, l), the sum over m ∈ {-l, -l + 1, • • • , l -1, l} of the first order corrections is zero. The spacing between the energies does not follow any scale, only the order between the levels is respected. A similar observation can be made for the first order corrected eigenfunctions,

l = 0 l = 1 l = 2 l = 3 E 4,1,±1 = -µZ 2 32 + µ 3 Z 4 1920a E 4,2,±2 = -µZ 2 32 + µ 3 Z 4 6760a E 4,3,±3 = -µZ 2 32 + µ 3 Z 4 16128a n = 4 E 4,0,0 = -µZ 2 32 E 4,3,±2 = -µZ 2 32 E 4,3,±1 = -µZ 2 32 -µ 3 Z 4 26880a E 4,3,0 = -µZ 2 32 -µ 3 Z 4 20160a E 4,2,±1 = -µZ 2 32 -µ 3 Z 4 13440a E 4,2,0 = -µZ 2 32 -µ 3 Z 4 6760a E 4,1,0 = -µZ 2 32 -µ 3 Z 4 960a ....
E 3,1,±1 = -µZ 2 18 + µ 3 Z 4 810a E 3,2,±2 = -µZ 2 18 + µ 3 Z 4 2835a n = 3 E 3,0,0 = -µZ 2 18 E 3,2,±1 = -µZ 2 18 -µ 3 Z 4 5670a E 3,2,0 = -µZ 2 18 -µ 3 Z 4 2835a E 3,1,0 = -
E 2,1,±1 = -µZ 2 8 + µ 3 Z 4 240a n = 2 E 2,0,0 = -µZ 2 8 E 2,1,0 = -µZ 2 8 -µ 3 Z 4 120a ....
ψ n,l,m = ψ (0) n,l,m + (n ,l ) =(n,l) -π 5 Z a ψ (0) n ,l ,m | Y 2,0 (θ,φ) r 3 |ψ (0) n,l,m E (0) n,l,m -E (0) n ,l ,m ψ (0) n ,l ,m = ψ (0) n,l,m + (n ,l ) =(n,l) -Z a +∞ 0 R n ,l (r)R n,l (r) r dr × (-1) m (l + 1 2 )(l + 1 2 )   l l 2 0 0 0     l l 2 -m m 0   -µZ 2 2n 2 + µZ 2 2n 2 ψ (0) n ,l ,m = ψ (0) n,l,m + (n ,l ) =(n,l) 2n 2 n 2 aµZ(n 2 -n 2 ) (-1) m (l + 1 2 )(l + 1 2 )   l l 2 0 0 0     l l 2 -m m 0   +∞ 0 R n ,l (r)R n,l (r) r dr ψ (0) n ,l ,m , (18)
the radial integral being proportional to µ 3 Z 3 (because of the 1 r 3 -term, as it appears when one makes the change of variables x = 2µZr), all the coupling coefficients in the expansion are proportional to µ 2 Z 2 a . Here also, the potentially divergent terms in the expansion, due to the degeneracy of the zero th order eigenvalues for a given n in the denominators, can be excluded since the corresponding numerators cancel out for the same reason that makes the first order matrix, Eq.( 15 Table 3: First-order corrected wave functions, ordered in increasing energy eigenvalue up to n=3 (in Appendix we provide a more comprehensive table up to n=7, "i-orbitals").

n = l = 0 m = 0 ψ 1,0,0 = ψ (0) 1,0,0 + √ 5 5aµZ n >2 1 1-1 n 2 +∞ 0 R n ,2 (r)R 1,0 (r) r dr ψ (0) n ,2,0 = ψ (0) 1,0,0 + µ 2 Z 2 a √ 6 480 ψ (0) 3,2,0 + 104 28125 ψ (0) 4,2,0 + 5 √ 14 6804 ψ (0) 5,2,0 + 5744 √ 21 12353145 ψ (0) 6,2,0 + 3299 √ 21 8847360 ψ (0) 7,2,0 + • • • n = l = 0 m = 0 ψ 2,0,0 = ψ (0) 2,0,0 + 4 √ 5 5aµZ n >2 1 1-4 n 2 +∞ 0 R n ,2 (r)R 2,0 (r) r dr ψ (0) n ,2,0 = ψ (0) 2,0,0 -µ 2 Z 2 a 32 √ 3 9375 ψ (0) 3,2,0 + 2 √ 2 1215 ψ (0) 4,2,0 + 3680 √ 7 7411887 ψ (0) 5,2,0 + 29 √ 42 215040 ψ (0) 6,2,0 + 63584 √ 42 645700815 ψ (0) 7,2,0 + • • • n = l = 1 m = 0 ψ 2,1,0 = ψ (0) 2,1,0 + 8 5aµZ n >2 1 1-4 n 2 +∞ 0 R n ,1 (r)R 2,1 (r) r dr ψ (0)
n ,1,0

+ 12 √ 21 35aµZ n >3 1 1-4 n 2 +∞ 0 R n ,3 (r)R 2,1 (r) r dr ψ (0) n ,3,0 = ψ (0) 2,1,0 + µ 2 Z 2 a 192 3125 ψ (0) 3,1,0 + 56 √ 10 6075 ψ (0) 4,1,0 + 8768 √ 5 1058841 ψ (0) 5,1,0 + 201 √ 35 89600 ψ (0) 6,1,0 + 2926784 √ 14 1076168025 ψ (0) 7,1,0 + • • • + µ 2 Z 2 a 8 √ 10 14175 ψ (0) 4,3,0 + 512 √ 5 823543 ψ (0) 5,3,0 + 123 √ 10 358400 ψ (0) 6,3,0 + 541184 √ 3 1076168025 ψ (0) 7,3,0 + • • • n = l = 1 m = ±1 ψ 2,1,±1 = ψ (0) 2,1,±1 - 4 5aµZ n >2 1 1-4 n 2 +∞ 0 R n ,1 (r)R 2,1 (r) r dr ψ (0) n ,1,±1 + 12 √ 14 35aµZ n >3 1 1-4 n 2 +∞ 0 R n ,3 (r)R 2,1 (r) r dr ψ (0) n ,3,±1 = ψ (0) 2,1,±1 -µ 2 Z 2 a 96 3125 ψ (0) 3,1,±1 + 28 √ 10 6075 ψ (0) 4,1,±1 + 4384 √ 5 1058841 ψ (0) 5,1,±1 + 201 √ 35 179200 ψ (0) 6,1,±1 + 1463392 √ 14 1076168025 ψ (0) 7,1,±1 + • • • + µ 2 Z 2 a 16 √ 15 42525 ψ (0) 4,3,±1 + 512 √ 30 2470629 ψ (0) 5,3,±1 + 41 √ 15 179200 ψ (0) 6,3,±1 + 541184 √ 2 1076168025 ψ (0) 7,3,±1 + • • • n = l = 0 m = 0 ψ 3,0,0 = ψ (0) 3,0,0 + 27 √ 5 35aµZ n >3 1 1-9 n 2 +∞ 0 R n ,2 (r)R 3,0 (r) r dr ψ (0) n ,2,0 = ψ (0) 3,0,0 -µ 2 Z 2 a 792 √ 3 588245 ψ (0) 4,2,0 + 45 √ 42 229376 ψ (0) 5,2,0 + 656 √ 7 2066715 ψ (0) 6,2,0 + 1809 √ 7 7812500 ψ (0) 7,2,0 + • • • n = l = 1 m = 0 ψ 3,1,0 = ψ (0) 3,1,0 + 18 5aµZ n >1 n =3 1 1-9 n 2 +∞ 0 R n ,1 (r)R 3,1 (r) r dr ψ (0) n ,1,0 + 27 √ 21 35aµZ n >3 1 1-9 n 2 +∞ 0 R n ,3 (r)R 3,1 (r) r dr ψ (0) n ,3,0 = ψ (0) 3,1,0 + µ 2 Z 2 a -192 3125 ψ (0) 2,1,0 + 60288 √ 10 2941225 ψ (0) 4,1,0 + 231 √ 5 16384 ψ (0) 5,1,0 + 35648 √ 35 10333575 ψ (0) 6,1,0 + 310191 √ 14 78125000 ψ (0) 7,1,0 + • • • -µ 2 Z 2 a 10368 √ 10 20588575 ψ (0) 4,3,0 + 81 √ 5 458752 ψ (0) 5,3,0 + 128 √ 10 3444525 ψ (0) 6,3,0 + 567 √ 3 39062500 ψ (0) 7,3,0 + • • • n = l = 1 m = ±1 ψ 3,1,±1 = ψ (0) 3,1,±1 - 9 5aµZ n >3 1 1-9 n 2 +∞ 0 R n ,1 (r)R 3,1 (r) r dr ψ (0) n ,1,±1 + 27 √ 14 35aµZ n >3 1 1-9 n 2 +∞ 0 R n ,3 (r)R 3,1 (r) r dr ψ (0) n ,3,±1 = ψ (0) 3,1,±1 + µ 2 Z 2
√ 2 39062500 ψ (0) 7,3,±1 + • • • n = l = 2 m = 0 ψ 3,2,0 = ψ (0) 3,2,0 + 9 √ 5 5aµZ n >0 n =3 1 1-9 n 2 +∞ 0 R n ,0 (r)R 3,2 (r) r dr ψ (0)
n ,0,0

+ 18 7aµZ n >3 1 1-9 n 2 +∞ 0 R n ,2 (r)R 3,2 (r) r dr ψ (0)
n ,2,0 (in harmonic quantum level n vib ). Hence, it is hoped that these approximate eigenstates could be appropriate to describe the electron density of the effective electronic

+ 54 √ 5 35aµZ n >4 1 1-9 n 2 +∞ 0 R n ,4 (r)R 3,2 (r) r dr ψ (0) n ,4,0 = ψ (0) 3,2,0 + µ 2 Z 2 a - √ 6 480 ψ (0) 1,0,0 + 32 √ 3 9375 ψ (0) 2,0,0 + 128 √ 6 1764735 ψ (0) 4,0,0 + 5 √ 30 196608 ψ (0) 5,0,0 + 32 295245 ψ (0) 6,0,0 + 392 √ 42 29296875 ψ (0) 7,0,0 + • • • + µ 2 Z 2 a 13824 √ 6 4117715 ψ (0) 4,2,0 + 675 √ 21 802816 ψ (0) 5,2,0 + 9472 √ 14 14467005 ψ (0) 6,2,0 + 7371 √ 14 15625000 ψ (0) 7,2,0 + • • • + µ 2 Z 2 a 135 √ 105 1605632 ψ (0) 5,4,0 + 512 √ 42 4822335 ψ (0) 6,4,0 + 7749 √ 231 214843750 ψ (0) 7,4,0 + • • • n = l = 2 m = ±1 ψ 3,2,±1 = ψ (0) 3,2,±1 + 9 7aµZ n >3 1 1-9 n 2 +∞ 0 R n ,2 (r)R 3,2 (r) r dr ψ (0) n ,2,±1 + 9 √ 6 7aµZ n >4 1 1-9 n 2 +∞ 0 R n ,4 (r)R 3,2 (r) r dr ψ (0) n ,4,±1 = ψ (0) 3,2,±1 + µ 2 Z 2 a 6912 √ 6 4117715 ψ (0) 4,2,±1 + 675 √ 21 1605632 ψ (0) 5,2,±1 + 4736 √ 14 14467005 ψ (0) 6,2,±1 + 7371 √ 14 31250000 ψ (0) 7,2,±1 + • • • + µ 2 Z 2 a 675 √ 14 3211264 ψ (0) 5,4,±1 + 512 √ 35 4822335 ψ (0) 6,4,±1 + 7749 √ 770 429687500 ψ (0) 7,4,±1 + • • • n = l = 2 m = ±2 ψ 3,2,±2 = ψ (0) 3,2,±2 -18 √ 5 35aµZ n >3 1 1-9 n 2 +∞ 0 R n ,2 (r)R 3,2 (r) r dr ψ (0) n ,2,±2 + 9 √ 3 7aµZ n >4 1 1-9 n 2 +∞ 0 R n ,4 (r)R 3,2 (r) r dr ψ (0) n ,4,±2 = ψ (0) 3,2,±2 -µ 2 Z 2 a 13824 √ 6 4117715 ψ (0) 4,2,±2 + 675 √ 21 802816 ψ (0) 5,2,±2 + 9472 √ 14 14467005 ψ (0) 6,2,±2 + 7371 √ 14 15625000 ψ (0) 7,2,±2 + • • • + µ 2 Z 2 a 675 √ 7 3211264 ψ (
Hamiltonians solved in the EN-MFCI method [START_REF] Cassam-Chenaï | ??A920125022015[END_REF] .

Although the model Hamiltonian depends upon 3 parameters, µ, Z, and a, its first order approximate eigenstates are parametrized by only two independent ones. They can conveniently be chosen as, c = µZ and b = µZ a . Furthermore, since the product b × c factorizes the first order perturbative term in the approximate DHOs expansion, Eq.( 18), we propose to just add the functions made of these corrective terms to standard basis sets, and leave it to molecular orbital calculations to optimize their linear combination with the other orbitals.

More precisely, instead of using ψ n,l,m (b, c) of Eq.( 18) (making explicit the dependency upon the free parameters), we propose to use

ψ (1) n,l,m (c) b • c = (-1) m 16n 2 l + 1 2 (n ,l ) =(n,l) n 2 (n 2 -n 2 ) l + 1 2   l l 2 0 0 0     l l 2 -m m 0   +∞ 0 R n ,l (x)R n,l (x) x dx ψ (0) n ,l ,m (c), (19)
which only depends upon c through the ψ An important issue is to calculate efficiently multicenter integrals involving the χ i 's and

DHO

(1) j 's orbitals, the former being usually GTOs and the latter ETOs. We will adopt a strategy similar to that of the SMILES module 14 , optionally included in the MOLPRO quantum chemistry package 15 to deal with STOs. That is to say, the DHO (1) 's will be transformed into an GTO-expansion, as also in Ref. [START_REF] Shaw | [END_REF] . However, different Cartesian prefactors will be mixed in the GTO-expansion of a given DHO. The BDF code [17][18][19] has been modified to accept contracted GTO having such mixed Cartesian prefactors. Full computational details will be provided in the second paper of this series 20 . Note that, a recently proposed, alternative strategy could be considered in the future 21 .

A hyperbolic cosine factor can be associated to ETO to provide a "double zeta" character to a minimal basis set 22 . This can be considered for DHO, as well. However, it seems more natural within our framework to combine different sets of DHO corresponding to atoms in different ionization states, to obtain multi-zeta basis sets. This last option preserves the desirable property for the DHOs to have associated energy eigenvalues, which can be used for truncation purposes.

Beside their application to the EN-MFCI method, which will be investigated in the next paper of the series, DHOs could also be useful for "clamped nuclei" quantum chemistry calculations, as they can be regarded as naturally "sp-hybridized" (unlike spherical harmonics STO). Actually, one could develop other model Hamiltonian distorted from spherical symmetry in more than one direction and obtain their first order approximate eigenfunctions by following the approach presented in this paper. For example, generalization of our con- Appendix: First order corrected eigenstates up to n=7

n = 1 E 0 1 = -µZ 2 2 l = 0 E 1 1,0,0 = 0 n = 2 E 0 2 = -µZ 2 8 l = 0 E 1 2,0,0 = 0 l = 1 E 1 2,1,0 = -µ 3 Z 4 120a E 1 2,1,±1 = µ 3 Z 4 240a n = 3 E 0 3 = -µZ 2 18 l = 0 E 1 3,0,0 = 0 l = 1 E 1 3,1,0 = -µ 3 Z 4 405a E 1 3,1,±1 = µ 3 Z 4 810a l = 2 E 1 3,2,0 = -µ 3 Z 4 2835a E 1 3,2,±1 = -µ 3 Z 4 5670a E 1 3,2,±2 = µ 3 Z 4 2835a n = 4 E 0 4 = -µZ 2 32 l = 0 E 1 4,0,0 = 0 l = 1 E 1 4,1,0 = -µ 3 Z 4 960a E 1 4,1,±1 = µ 3 Z 4 1920a l = 2 E 1 4,2,0 = -µ Z 6720a E 1 4,2,±1 = -µ Z 13440a E 1 4,2,±2 = µ Z 6720a l = 3 E 1 4,3,0 = -µ 3 Z 4 20160a E 1 4,3,±1 = -µ 3 Z 4 26880a E 1 4,3,±2 = 0 E 1 4,3,±3 = µ 3 Z 4 16128a n = 5 E 0 5 = -µZ 2 50 l = 0 E 1 5,0,0 = 0 l = 1 E 1 5,1,0 = -µ 3 Z 4 1875a E 1 5,1,±1 = µ 3 Z 4 3750a l = 2 E 1 5,2,0 = -µ 3 Z 4 13125a E 1 5,2,±1 = -µ 3 Z 4 26250a E 1 5,2,±2 = µ 3 Z 4 13125a l = 3 E 1 5,3,0 = -µ 3 Z 4 39375a E 1 5,3,±1 = -µ 3 Z 4 52500a E 1 5,3,±2 = 0 E 1 5,3,±3 = µ 3 Z 4 31500a l = 4 E 1 5,4,0 = -µ 3 Z 4 86625a E 1 5,4,±1 = -17µ 3 Z 4 1732500a E 1 5,4,±2 = -2µ 3 Z 4 433125a E 1 5,4,±3 = µ 3 Z 4 247500a E 1 5,4,±4 = µ 3 Z 4 61875a n = 6 E 0 6 = -µZ 2 72 l = 0 E 1 6,0,0 = 0 l = 1 E 1 6,1,0 = -µ 3 Z 4 3240a E 1 6,1,±1 = µ 3 Z 4 6480a l = 2 E 1 6,2,0 = -µ 3 Z 4 22680a E 1 6,2,±1 = -µ 3 Z 4 45360a E 1 6,2,±2 = µ 3 Z 4 22680a l = 3 E 1 6,3,0 = -µ 3 Z 4 68040a E 1 6,3,±1 = -µ 3 Z 4 90720a E 1 6,3,±2 = 0 E 1 6,3,±3 = µ 3 Z 4 54432a l = 4 E 1 6,4,0 = -µ 3 Z 4 149688a E 1 6,4,±1 = -17µ 3 Z 4 2993760a E 1 6,4,±2 = -µ 3 Z 4 374220a E 1 6,4,±3 = µ 3 Z 4 427680a E 1 6,4,±4 = µ 3 Z 4 106920a l = 5 E 1 6,5,0 = -µ 3 Z 4 277992a E 1 6,5,±1 = -µ 3 Z 4 308880a E 1 6,5,±2 = -µ 3 Z 4 463320a E 1 6,5,±3 = -µ 3 Z 4 2779920a E 1 6,5,±4 = µ 3 Z 4 463320a E 1 6,5,±5 = µ 3 Z 4 185328a n = 7 E 0 7 = -µZ 2 98 l = 0 E 1 7,0,0 = 0 l = 1 E 1 7,1,0 = -µ 3 Z 4 5145a E 1 7,1,±1 = µ 3 Z 4 10290a l = 2 E 1 7,2,0 = -µ 3 Z 4 36015a E 1 7,2,±1 = -µ 3 Z 4 72030a E 1 7,2,±2 = µ 3 Z 4 36015a l = 3 E 1 7,3,0 = -µ 3 Z 4 108045a E 1 7,3,±1 = -µ 3 Z 4 144060a E 1 7,3,±2 = 0 E 1 7,3,±3 = µ 3 Z 4 86436a l = 4 E 1 7,4,0 = -µ 3 Z 4 237699a E 1 7,4,±1 = -17µ 3 Z 4 4753980a E 1 7,4,±2 = -2µ 3 Z 4 1188495a E 1 7,4,±3 = µ 3 Z 4 679140a E 1 7,4,±4 = µ 3 Z 4 169785a l = 5 E 1 7,5,0 = -µ 3 Z 4 441441a E 1 7,5,±1 = -µ 3 Z 4 490490a E 1 7,5,±2 = -µ 3 Z 4 735735a E 1 7,5,±3 = -µ 3 Z 4 4414410a E 1 7,5,±4 = µ 3 Z 4 735735a E 1 7,5,±5 = µ 3 Z 4 294294a l = 6 E 1 7,6,0 = -µ Z 735735a E 1 7,6,±1 = -µ Z 792330a E 1 7,6,±2 = -µ Z 1030029a E 1 7,6,±3 = -µ Z 2060058a E 1 7,6,±4 = µ Z 5150145a E 1 7,6,±5 = µ Z 936390a E 1 7,6,±6 = µ Z 468195a
Table 4: Zero order eigenvalues and first order corrections up to n=7

n = 1 < 100|320 >=
n,l,m =

(n ,l ) =(n,l)

< n l m|nlm > µ 2 Only positive m-values are tabulated, their negative counterparts giving the same coupling coefficients. Note that, all numbers being real, < n l m|nlm >=< nlm|n l m >

  non zero, and, (iii) n > l . The allowed quantum number values are summed up in Tab.1.

  n can be expanded as a linear combination of generalized Laguerre polynomial L 2l +1 i 2µZr n with i ≤ n -l -3. All these polynomials are orthogonal to L 2l +1 n-l -1 2µZr n for the mesure 2µZr n (2l +1) exp -2µZr n which is proportional

  ), diagonal. The first terms in the expansion of the lowest eigenfunctions are given in Tab. 3. Isodensity surfaces of the density functions corresponding to these eigenfunctions are displayed in Figs. 1 to 10. They shows isosurfaces belonging sometimes to different homology groups for different density values, which give an idea of the repartition of the electronic charge as it becomes more and more concentrated. The parameter c := µ × Z, which controls how diffuse the unperturbed hydrogenoid orbitals are, is set to 1. The parameter b := µ×Z a , which tunes the intensity of the distorsion with respect to spherical symmetry, takes two values: b = 1 close to the no distorsion case (b = 0), and b = 32 for which the effect of the perturbed potential becomes clearly apparent. The perturbed 1s orbital, Fig. 1 for example, becomes elongated along the Coulomb potential spreading axis, as one expects. The variations for other orbitals are sometimes less intuitive.

  ,±2 + • • • 4 Application to quantum chemical calculations In the previous section, we have obtained first order perturbative corrections to the eigenfunctions of a smeared Coulomb potential along an arbitrary axis. Such a potential can represent the average potential felt by an electron bounded to a vibrating nucleus of effective charge Z, effective reduced mass µ and effective classical vibrating amplitude equal to 2n vib +1 2a

  coefficients in the occupied molecular orbitals, obtained for example in a HF calculation, implicitly determine the corresponding coefficient products, b • c, and consequently, if c is given, the free parameters b. For instance, in H 2 , there is only one occupied HF orbital of the form, ψ 1 = i u i χ i + j v j DHO (1) j where the χ i 's are standard basis functions, and the DHO (1) j are various DHOs on different centers and/or different quantum numbers, indexed by j. The u i , v j 's are scalar coefficients. One can avoid the explicit optimization of the b j • c j values for each DHO j (c j ), and just use the v j 's values.

volutedFigure 1 :Figure 2 :Figure 3 :Figure 4 :Figure 5 :Figure 7 :Figure 8 :Figure 10 :

 123457810 Figure 1: Comparison of hydrogen smeared Coulomb potential 1s orbital squared density levels.

  ................................................................................................................

  µZ 2 18 -µ 3 Z 4 405a ....................................................................................................................

  ................................................................................................................

	n = 1	E 1,0,0 = -µZ 2 2

Table 5 :

 5 Non-zero, first order coupling coefficients up to n=7 in µ 2

	units,
	Z 2	a

Acknowledgment

The authors were partially supported by ERC grant SCAPDE, 320845, and by the Döblin fédération de Recherches of CNRS.