Hongxia Wang

Meiyan Cui

Gaoliang Fu

Jiaye Zhang

Xingyu Ding

Irene Azaceta

Matthieu Bugnet

Demie M Kepaptsoglou

Vlado K Lazarov

Victor A De

La Pena O'shea

Kelvin H L Zhang
email: kelvinzhang@xmu.edu.cn

Shea Vad , O'

Freddy E Oropeza
email: freddy.oropeza@imdea.org

Victor De La Peña O'shea

Vertically aligned Ni/NiO nanocomposites with abundant oxygen deficient hetero-interfaces for enhanced overall water splitting

Keywords: Water splitting, Hydrogen evolution reaction, Electronic structure, Heterostructure, Vertically aligned nanocomposites Speciality: Inorganic Chemistry Water splitting, Hydrogen evolution reaction, Electronic structure, Heterostructure, Vertically aligned nanocomposites

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Electrolysis of water powered by renewable electricity to produce hydrogen is widely considered as a promising pathway for global clean and sustainable energy supply in the future. [1][2][3][4] The key to enable this technology is the development of low-cost and efficient electrocatalysts with enhanced activity and stability to accelerate the water splitting reactions, i.e., hydrogen evolution reaction (HER) and oxygen evolution reaction (OER). [5] At present, the most effective HER and OER electrocatalysts are Pt/Pd and RuO 2 /IrO 2 based materials, respectively. [6][7][8][9][10][11] Although these noble metal based catalysts exhibit high performance, their large-scale applications are limited by high costs and scarcity. Therefore, great efforts have been devoted to investigate earth-abundant transition metals (Co, Ni, Fe) based materials. [12][13][14] However, most of them generally perform well only for the half-cell reaction (either HER or OER) in either alkaline or acidic media. [15,16] Therefore, developing highly efficient bifunctional electrocatalysts which can simultaneously catalyze HER and OER process in the same electrolyte, is technologically crucial to increase the overall efficiency of water electrolysis.

Up to now, there are two mature water electrolysis systems: proton exchange membrane water electrolyser (PEMWE) performed in acidic condition, and alkaline water electrolyser (AWE) worked in alkaline media. [17] PEMWE has the advantages of high current densities, high efficiency and high gas purity. However, it faces problem with the scaling-up due to the high cost of the components like Pt/Pd and RuO 2 /IrO 2 based electrocatalyst materials. [18] To this regard, AWE electrolyser, which utilizes low-cost transition metal-based catalysts, is more appealing for largescale industrialization. Additionally, in alkaline media, most OER catalysts show faster reaction kinetics and better stability compared to those in acidic solution. However, so far the alkaline HER still suffers from sluggish reaction kinetics due to the additional water dissociation step, which decreases the efficiency of the overall process. Even for Pt, the HER performance in alkaline condition is limited by slow water dissociation due to the lack of optimal active sites for H-OH bond cleaving. [19,20] Motivated by this, designing active bifunctional catalysts capable of appropriate adsorption and transformation of hydrogen and oxygen involved intermediates for both OER and HER in alkaline media is highly desirable. Encouragingly, some studies reveal that multifunctional heterojunctions could accelerate the water dissociation step and achieve high intrinsic activity in alkaline media. For example, Subbaraman et al proposed Pt/Ni(OH) 2 heterostructured catalyst with outstanding HER performance, in which the Ni(OH) 2 clusters promote water dissociation, and Pt serves as the active site to generate molecular hydrogen. [21] Subsequently, considerable research efforts have been devoted to investigating heterostructure catalysts, which show superior multifunctional catalytic activities for HER, OER or both. [22,23] For instance, Xie and co-workers reported a Fe 3 O 4 -blended Ni/NiO heterostructure catalyst, exhibiting high performance toward both OER and HER, due to the synergetic effects between Ni/NiO compound and Fe 3 O 4 . [24] Although recent advances in nanostructured materials have led to the development of high-performance catalysts for HER and OER, the real active sites and reaction mechanism are still ambiguous. [25,26] Most of the previously reported heterostructured catalysts are based on nanostructures with different morphologies, exposed crystal facets, and defects, which would complicate the elucidation of possible active sites and the reaction mechanism. [27,28] Recently, vertically aligned nanocomposites (VAN), [29,30] which possess two or more epitaxially grown phases on substrates, open up new possibilities for the design of heterostructured catalysts. [31] VAN materials provide abundant well-defined hetero-interface on the surface with precise control over the crystal

F o r R e v i e w O n l y
orientation, roughness and thickness, therefore avoiding structural and morphological ambiguity inherent to the preparation of free standing powder catalysts. [32,33] Such VAN catalysts can provide ideal platforms to characterize the structural and electronic properties of multifunctional heterostructures, and to study the reaction mechanism. In addition, VAN can also significantly enhance the electrical conductivity to boost overall water splitting performance owing to the strong interfacial interaction between the components of the composite and the substrate, which minimizes issues of phase segregations and mechanical peeling of active components. Although VAN have been used as model catalysts for several chemical reactions, they are rarely applied for electrochemical water splitting.

In this work, we show that Ni/NiO VAN are efficient catalyst for both the alkaline HER and OER, owing to synergetic effect at the Ni/NiO interface that significantly reduce the energy barrier for water dissociation, and favors the formation of reactive key H* intermediates on the Ni side of the interface, and OH ads on the NiO side of the interface. For our study, vertically aligned Ni/NiO nanocomposite catalysts on the conductive (001) oriented Nb-doped SrTiO 3 (Nb:STO) substrates were prepared using pulsed laser deposition (PLD). Different from previously reported Ni-based heterostructure catalysts that are prepared by hydrothermal, [34] electrodeposition [35] or thermal annealing [36] NiO and Ni/NiO VAN thin films were epitaxially grown on Nb:STO(001) substrates by PLD with a 248 nm KrF excimer from NiO target at the oxygen partial pressure of 10 and 0.01 mTorr, respectively. The NiO target was made from NiO powder (99.998%, Alfa Aesar) by pelletizing and sintering in air at 1100 C for 12 hours. During the growth of NiO, the laser repetition rate and energy density, and the substrate temperature were set as 5 Hz, 1.2 J cm -2 and 600 °C, respectively.

Experimental Section

Ni thin films on Nb:STO(001) substrates were synthesized via magnetron sputtering with a thickness of ∼50 nm. Sputtering was performed at a rf power of 200 W and the Ar pressure of 3.8 mTorr for direct deposition from a Ni target.

Characterizations.

The crystal structure of NiO, Ni/NiO VAN and Ni was examined by high-resolution XRD using a PANalytical four-circle diffractometer in θ-2θ scans and reciprocal space mapping modes. Atomic force microscopy (AFM) with tapping-mode was applied to characterize surface morphologies. The detailed structure of Ni and Ni/NiO VAN were characterized by scanning transmission electron microscopy (STEM) and electron energy-loss spectroscopy (EELS) measurements, performed in a Nion UltraSTEM100 equipped with a Gatan Enfina spectrometer and a Nion UltraSTEM100MC monochromated instrument equipped with a Gatan Enfinium spectrometer. Both instruments were operated at 100kV, with electron optics set to form a 0.9 Å probe. Cross-sectional transmission electron microscopy (TEM) specimen preparation was carried out by focused ion beam (FIB) milling using an FEI Helios NanoLab 600 instrument. The x-ray absorption spectroscopy (XAS) spectra were conducted under the total electron yield (TEY) mode with an effective resolution of 0.2 eV at the I09 Surface and Interface Structural Analysis beamline of the Diamond Light Source (DLS). And X-ray adsorption near edge structure (XANES) measurements at the Ni K edges were performed at the XAFCA beamline of Singapore Synchrotron Light Source (SSLS). X-ray photoemission spectroscopy (XPS) were characterized by a monochromatic Al Kα X-ray (hʋ= 1486.6 eV) source at normal emission (electron take-off angle=90º relative to the surface plane) with a SPECS PHOIBOS 150 electron energy analyzer. Near ambient pressure XPS spectra were obtained through a lab-based spectrometer with a monochromated Al Kα source. A SPECS

F o r R e v i e w O n l y 6
PHOIBOS 150 NAP analyser was used. The entrance to the analyser is a nozzle with a 300 mm diameter orifice. The total energy resolution of the measurements was about 0.50 eV.

Electrochemical analysis.

The NiO, Ni/NiO VAN and Ni thin film electrodes were prepared by the method previously reported. Nb:STO substrates are conductive so that the electrical contact attached to a narrow copper stripe was made to the back of the NiO, Ni/NiO VAN and Ni thin films (instead of the front of thin films) using GaIn eutectic and Ag paste. The sides of thin films and copper stripe were covered with epoxy. The catalyst surface is calculated according to the actual exposed surface area. The CHI 660E electrochemical workstation was used for the electrochemical measurements for OER in 0.1M KOH using a typical three-electrode system. A Hg/HgO (0.1 M KOH) electrode and carbon rob were worked as the reference electrode and counter electrode, respectively. The electrolyte was prepared by high-purify KOH (Adamas, 99.999%) and Milli-Q water (18.2 MΩ-cm). In order to detach the impurity on the surface of electrode, cyclic voltammograms (CV) were firstly performed under the potential range of 0 to 0.60 V (vs RHE). Then, linear sweep voltammetry (LSV) polarization curves for OER were carried out with the potential range of -0.1 to 0.8 V (vs Hg/HgO) at a scan rate 10 mV s-1. The electrochemical impedance spectroscopy (EIS) measurement was conducted from 100 kHz to 1 Hz at a potential of 0.7 V (vs Hg/HgO). All potential values were converted to reversible hydrogen electrode (RHE) according to the Nernst equation: E RHE =E Hg/HgO + 0.059 pH + 0.164 V.

Results and Discussion

Synthesis of vertically aligned Ni/NiO heterointerfaces

To realize the growth of two phases, the PLD growth was conducted at a low oxygen partial pressure of 0.1 mTorr. Similar to the growth process of other self-assembled nanocomposites, [37] the growth of Ni/NiO VAN may go through a three-stage self-decomposition mechanism. [38] In and NiO (-111), respectively. The difference in lattice parameter between Ni and NiO results in a 3% lattice mismatch at the Ni/NiO interface, which may lead to the generation of oxygen vacancies at the Ni/NiO interface. [41] The chemical compositions of NiO, Ni/NiO VAN and Ni were studied by XAS and XPS. Figure 2a shows the XANES)at the Ni-K edge of each sample. The pre-edge peak at ~8330 eV is assigned to the electron transition of 1s to 3d state, [42] which is observed in NiO, but absent in metallic In order to investigate the chemical and electronic information of the Ni/NiO VAN interface, high energy-resolution EELS was conducted in the three regions (as shown in Figure 2c), including the Ni region (blue box), the region of Ni/NiO interface (pink box) and the NiO region (gray box). The Ni-L 2 and -L 3 edges, due to the electronic transition from 2p 1/2 and 2p 3/2 to unoccupied 3d state, are highly sensitive to changes in the occupancy of d states. [46] The energy position, sub-splitting of the L 3 and L 2 peaks, as well as their relative intensity ratio can be related to the valence state and O coordination of Ni. [47,48] Comparing the fine structures of Ni-L edges in the three regions marked in Figure 2d

Electrocatalytic Performance

The HER and OER activities of NiO, Ni/NiO VAN and Ni thin films were evaluated in a standard three-electrode setup in 0.1 M KOH solution, which have been commonly used for model thin film electrocatalysts. [54][55][56] Figure 3a presents the first linear sweep voltammetry (LSV) curves for the HER for each sample. The Ni/NiO VAN delivers the best HER catalytic performance. As shown in Figure 3b, the Ni/NiO VAN requires a low overpotential of 0.196 V to reach 10 A cm -2 , which is significantly lower than the overpotentials of metallic Ni (0.302 V) and NiO (0.420 V). To understand the kinetics of HER, the corresponding Tafel plots are calculated from the linear region of the polarization curves. It is widely accepted that Tafel slopes values around 120, 40 and 30 mV dec -1 are observed when the alkaline HER rate is determined by the Volmer, Heyrovsky and Tafel steps, respectively. As shown in Figure 3c, the Tafel slopes of the HER for all electrode samples have values above 100 mV dec -1 , suggesting that Volmer step (the water dissociation step) is the rate-determining step for the reaction driven by all studied catalyst. However, the Tafel slope of Ni/NiO VAN has a significantly smaller value (108 mV dec -1 ) than that for the HER driven by metallic Ni (156 mV dec -1 ) and NiO (238 mV dec -1 ), which means a faster Volmer step and a higher hydrogen generation rate for Ni/NiO VAN. Furthermore, the electrochemical impedance spectroscopy (EIS) was conducted to investigate the electron-transfer HER kinetics. As shown in To investigate the electrochemically active surface areas (ECSA) of the samples, CV measurement was performed to obtain the double layer capacitance (C dl ). As shown in Figure S8, the Ni/NiO VAN shows the highest C dl value (46.8 μF cm -2 ) compared with Ni (17.3 μF cm -2 ) and NiO (25.6 μF cm -2 ), thereby suggesting the highest catalytically relevant surface area. However, in the HER process, Ni shows better performance than NiO, which is not consistent with the trend of C dl values. Thus, this approach used for determining the ECSA is limited. Moreover, the C dl can be affected by several factors, such as the conductivity, the surface coordination of ions, and capacitance originated from the trap states and residual charge-transfer processes in non-Faradaic region. [57,58] Thus, the ECSA obtained by C dl may not reflect the real active area of electrocatalysts and can serve as an additional parameter for the performance indication.

Mechanism Study: Origin of Enhanced OER and HER Performance

As revealed by the electrocatalytic performance tests, Ni/NiO VAN interface seems to favor the Volmer step, i.e., the water dissociation, during the HER process. To better understand the adsorption and dissociation of water on the Ni/NiO heterostructure, we investigated the reactivity of water vapor on Ni/NiO VAN and NiO surfaces using NAP-XPS. In the experiment, we monitored the O 1s region of the spectra of the sample as water vapor was dosed into the chamber at pressures ranging from 7×10 -7 mbar to 2.5 mbar. Figure 4a shows photoelectron spectra of Ni/NiO VAN and NiO in the O 1s region under vacuum and 0.6 mbar water vapor pressure (see supplementary Figure S10 for other water pressure). The spectra can be fitted with four curves that represent contributions from four oxygenated species, including lattice oxygen at ~530 eV, adsorbed OH groups at ~531.8 eV, chemisorbed H 2 O (H 2 O ads ) at ~533.2 eV, and gaseous H 2 O (H 2 O gas ) at ~535 eV, [59,60] and their relative amounts can be obtained from peak fit area (Figure 4b). In particular, the absorbed OH* originates from the dissociative H 2 O, and therefore can be used as a measure of the degree of water dissociation on catalyst surface, i.e., the higher the value, the better the ability to dissociate water. Figure 5b shows plots of the percentage of adsorbed hydroxyl as a function of water vapor pressure. The percentage of OH ads content remains within 15-25 % for NiO under different water pressure, but OH ads content increases dramatically up to ~50 % for Ni/NiO VAN. The considerable increase of the hydroxylation on the Ni/NiO surface suggests that the dissociative adsorption of water is more favorable on Ni/NiO surface, and thereby accelerates the Volmer step which is the rate determining step for the HER process in alkaline conditions. At the surface with oxygen vacancies, the undercoordinated metal centers near oxygen vacancies are more active towards the adsorption of H 2 O and OH to achieve a saturated coordination than that at the surface without oxygen vacancies. [5] Therefore, we propose that the oxygen vacancies existed at Ni/NiO interface may promote the observed dissociative adsorption of H 2 O to form chemisorbed OH ads and H ads (as shown in Figure 4c), also can promote the adsorption of OH ads on the catalyst surface, which contributes to the highly enhanced electrocatalytic activity for the both OER and HER. Strikingly, the water dissociation barrier at the Ni/NiO VAN interface is reduced to 1.53 eV on the NiO sites (Figure 5a), indicating that oxygen vacancies at the Ni/NiO interface can promote water dissociation, which is consistent with the NAP-XPS results. On the one hand, the generated OH * will heal the oxygen vacancies and proceed the next reaction to produce oxygen. According to the previous reports, [61] the rate determining step of NiO for OER is the adsorption of OH * , and we found that the presence of oxygen vacancies at the interface for Ni/NiO significantly enhance the adsorption of OH * . Therefore, we conclude the underlying reason for the highest OER performance of Ni/NiO VAN, in comparison with individual NiO and Ni catalysts.

Previous reported works have proved that the optimal value of Gibbs free energy of H* (ΔG H* ) is 0, where the adsorption and desorption of H* at the electrocatalyst can be well balanced for H 2 generation. [62] Therefore, the ΔG H* value of a catalyst could be reasonably utilized to judge the catalytic performance for HER (the closer to the ΔG H* value of 0, the better the HER activity). Based on that, we computed values ofΔG H* in order to analyze the optimum adsorption and recombination sites to form H 2 . As shown in Figure 5b

Conclusion

In summary, we show that the Ni/NiO nanocomposite is an efficient catalyst for both alkaline Spin polarized DFT calculations were performed using the Vienna ab initio simulation package (VASP). [1,2] The generalized gradient approximation proposed by Perdew, Burke, and Ernzerhof (GGA-PBE) is selected for the exchange-correlation potential. [3] The pseudo-potential was described by the projector-augmented-wave (PAW) method. [4] The geometry optimization is performed until the Hellmann-Feynman force on each atom is smaller than 0.03 eV•Å -1 . The energy criterion is set to 10 -6 eV in iterative solution of the Kohn-Sham equation. Ni(001) and NiO(001) slab models were used for adsorption studies. For the study of NiO/Ni heterojunction, NiO(-100) and Ni(100) surface were built together and optimized, and the (001) surface of the heterojunction was cleaved with 12 Å vacuum above the model. The bottom two layers of atoms were fixed for all optimizations. Nudged elastic band method (NEB) [5,6] was used to search for transition state structures.

  photoelectron spectroscopy (NAP-XPS) shows that the enhanced catalytic performances are correlated with a substantial increase of the water dissociation capability of Ni/NiO VAN. Based on density-functional theory (DFT) calculations, we propose that the existence of oxygen vacancies at the interface of Ni/NiO VAN promotes the water dissociation and adsorption of OH ads . This work reveals insights into the mechanism of Ni/NiO heterostructure for overall water splitting, and provide an effective strategy for the design of novel heterostructure catalysts with well-defined interface.

5 2. 1 .

 51 Thin films growth.

  the first stage, NiO adatoms vaporized from NiO target nucleate on the Nb:STO substrate. During the deposition process, a low oxygen pressure triggers the phase segregation of Ni from NiO to form ultra-fine Ni clusters. Next, thermal energy facilitates the coalescence of Ni clusters, forming Ni nanodomains. In the third stage, the phase segregation is stabilized, leading to a vertically aligned addition, NiO and metallic Ni thin films were also prepared as reference samples (see details in the Experimental Section).

Figure

  Figure 1a presents XRD -2 scans for NiO, Ni/NiO VAN and Ni thin films. Sharp Bragg peaks at 43.2° and 51.8° corresponding to NiO (002) and metallic Ni (002), respectively, are observed in the scan of Ni/NiO VAN, indicating the formation of Ni/NiO heterostructures. To gain additional structural information, reciprocal space maps (RSMs) of NiO and Ni/NiO VAN near the (113) reflection of Nb:STO were measured (as shown in Figure 1b and c), from which the in-plane and out-of-plane lattice parameters were extracted. According to the lattice parameters, the unit cell volume can be calculated. As displayed in Figure S1, the unit cell volume of NiO in Ni/NiO VAN is larger than that of pure NiO. The expansion of unit cell volume of NiO in Ni/NiO VAN may be caused by the generation of oxygen vacancies at the Ni/NiO interface.[39, 40] Cross-sectional STEM images were collected to further investigate microscopic structure. Figure 1 d and e show the large-scale cross-sectional STEM images of NiO and Ni/NiO VAN. Different from the continuous flat layer observed in the NiO film, the image of Ni/NiO VAN has brighter regions of metallic Ni due to the high Ni concentration.[38] Metallic Ni and NiO phases coexist with a fairly uniform distribution of Ni clusters of 20~30 nm embedded in NiO matrix. The HAADF-STEM image and corresponding energy-dispersive X-ray spectroscopy (EDX) elemental maps on the sample region are shown in Figure 1f. We can see that the Ni-NiO coupling leads to a single crystalline stacking and distinctive interface between Ni and NiO. The high-resolution HAADF-STEM image (Figure 2g and h) near the Nb:STO substrate viewed along the Nb:STO [-1-10] direction, shows the epitaxial growth of the two crystalline phases. The lattice plane spacing in the Ni domain near the interface, 2.08 Å and 2.35 Å, correspond to Ni(002) and Ni(-111), respectively. In the NiO domain, lattice plane spacing of 2.14 Å and 2.35 Å correspond to NiO (002)

  Ni.[43] Metallic Ni exhibits a pre-edge peak at ~8332 eV. The enhancement of this pre-edge peak can be regarded as an indication of the existence of metallic Ni. For Ni/NiO VAN, the pre-edge peak at ~8332 eV can also be observed, indicating that both NiO and metallic Ni presented in Ni/NiO film. This could be further supported by the line-shape of Ni 2p XPS of Ni/NiO VAN, shown in Figure2b, which exhibits spectral feature associated with both NiO in the range of 855~858 eV, and metallic Ni at around 853.0 eV (see Ni 2p XPS of pure Ni and NiO in FigureS2).

Figure

  Figure S3 shows the O 1s spectra of Ni/NiO VAN, NiO and Ni. The O 1s spectra of Ni/NiO VAN and NiO are deconvoluted into three peaks that represent contributions from three oxygenated species, including lattice oxygen (O Latt ) at ~529.6 eV, defective oxygen at ~531.1 eV (V o ), adsorbed OH groups (OH ads ) at ~532.3 eV, while the O 1s spectrum of Ni only shows a weak peak assigned to the adsorbed OH groups.[44, 45] The percentage of V o shows a high value (35.4 %) for the Ni/NiO VAN, proving the existence of abundant oxygen vacancies in the Ni/NiO VAN.

  , a clear change of the near edge fine structure can be observed as we scan from the NiO to the Ni/NiO and Ni regions: (1) the sub-splitting of the L 3 (high energy shoulder) and L 2 (doublet) peaks gradually decrease; (2) the splitting is accompanied by a gradual shift of the edge onset towards lower energies; (3) there is a decrease of the L 3 /L 2 intensity ratio. These changes suggest the presence of unsaturated coordinated Ni atoms at the Ni/NiO interface.[49,50] The results obtained from the Ni-L 2,3 fine structures are further confirmed by the analysis of O-K edge spectra.As is shown in Figure2e, peak A at ∼530 eV arises from the electronic transition from O 1s to unoccupied state with partial O 2p character hybridized with empty Ni 3d states.[51] Peak C at ~540 eV and peak D at ~550 eV correspond to the electronic transition from O 1s to the hybridization orbital of O 2p with Ni 4s and 4p, respectively.[52] The origin of peak B can be ascribed to the configuration interaction.[53] Therefore, the density of O-K edge is related to the oxygen coordination and the Ni oxidation state. It can be seen that the O-K edge of the NiO region exhibits sharp A-D peaks, while at the Ni/NiO interface, the relative intensity of peaks B and C become weaker, while a slight redshift of the edge onset can also be observed. These results can be assigned to the breaking of long range order and further indicate the presence of oxygen vacancies at the Ni/NiO interface.[51] Oxygen vacancies at the NiO side of the interface (Figure2f) result in the unsaturated coordinated coordination Ni atoms and a higher d band occupancy of Ni due to the nominal electron density donation from oxygen vacancies to the Ni cations.

Figure S4 ,

 S4 Figure S4, the charge transfer resistance (R ct ) for Ni/NiO VAN (290.4 Ω) is smaller than that of Ni (516.1 Ω) and NiO (3708.8 Ω), indicating the faster Faradaic process and HER kinetics of Ni/NiOVAN. The stability of samples was also measured by chronopotentiometry. As displayed in FigureS5, the potentials of all samples maintained a constant current density of -20 µA cm -2 don't exhibit

11 that

 11 Inspired by the good HER and OER performance of Ni/NiO VAN, we use Ni/NiO VAN as both the anode and cathode to assemble an electrolyzer for overall water splitting in alkaline condition.For comparison, NiO and Ni catalysts were also used to assemble the alkaline electrolyzer. The LSV curves of Ni/NiO VAN‖Ni/NiO VAN, Ni‖Ni and NiO‖NiO are presented in FigureS9. It can be seen Ni/NiO VAN only requires a low cell voltage of 1.59 V to deliver a current density of 10 μA cm -2 , which is smaller than that of Ni (1.68 V) and NiO (1.87 V). The better water splitting performance of Ni/NiO VAN is ascribed to the synergetic effect at Ni/NiO interfaces.

  To clarify the role of Ni/NiO hetero-interface on the adsorption of reaction intermediates and the enhanced catalytic activity, DFT calculations were employed to evaluate the energy barrier for breaking the H-OH bond in H 2 O and the chemisorption free energies of hydrogen (ΔG H* ) on the (001) surface of pure Ni, pure NiO, Ni/NiO VAN, and Ni/NiO without oxygen vacancy (denoted as w/o-Ni/NiO). The structure models are shown in Figure S11. As shown in Figure S12, the energy barriers of pure Ni, pure NiO and w/o-Ni/NiO (001) surface for breaking the OH-H bond are 3.27 eV, 2.27 eV, 1.81 eV, respectively. Such high energy barriers hinder the dissociation of water.

  , the ΔG H* values of NiO sites near the interface and pure Ni are 0.66 eV and -0.58 eV, respectively, indicating that the adsorption of hydrogen on NiO sites near the interface is too weak, and the adsorption of hydrogen on pure Ni sites is too strong. However, the metallic Ni sites near the interface possessed the optimal ΔG H* value of -0.16 eV, suggesting that the metallic Ni sites near the interface show the best HER performance. To further explore the mechanism of boosted alkaline HER activity of the Ni/NiO electrocatalysts, the Bader charge analysis on Ni cations of Ni/NiO was conducted together with the average valence charge of Ni cations along the direction of the Ni/NiO heterointerface (Figure 5c). Combining with differential the Ni/NiO electrocatalysts (the inset of Figure 5c), we can confirm the occurrence of charge transport nearby the Ni/NiO heterointerface, providing 0.33 e -of charge transfer from metallic Ni phase to NiO phase. The electrons transfer between Ni and NiO may be driven by the difference in work function of NiO (5.2~5.6 eV) and that of metallic Ni (5.0~5.2 eV),[63, 64] i.e., the Fermi level of NiO is lower than that of Ni. When Ni and NiO are brought into contacted, electrons in metallic Ni will spontaneously transfer to NiO. According to the d-band theory,[65] when the Ni atoms loss electrons and d band experiences a downward shift, more antibonding states of hydrogen are lower than the Fermi level, and hence weakens the affinity of Ni towards hydrogen. Therefore, the Ni sites near the interface exhibit more appropriate binding strength with H atoms, which improve the HER activity compared with pure metallic Ni. In view of the above results, we propose a mechanism for the enhanced HER in alkaline media on the Ni/NiO VAN electrocatalyst. As shown in Figure 5d, H 2 O dissociation is assumed to occur on the NiO sites near the Ni/NiO VAN interface, where the oxygen vacancies reduce the dissociative energy barrier and promote to form chemisorbed OH ads and H ads . Then the produced H* will migrate to the metallic Ni sites for the further combination to form H 2 .

Figure 1

 1 Figure 1 (a) XRD θ-2θ scans around (002) Bragg peaks for NiO, Ni/NiO VAN and Ni films grown on Nb:STO substrates. RSMs for NiO (b) and Ni/NiO VAN (c) around Nb:STO (113) reflection. Large-scale cross-sectional STEM images of Ni (d) and Ni/NiO VAN (e). (f) The HAADF-STEM image and corresponding EDX elemental maps of Ni/NiO VAN. (g) and (h) Cross-sectional HAADF-STEM micrograph at the Ni/NiO interface.

Figure 2 .

 2 Figure 2. (a) Ni-K edge XAS of NiO, Ni/NiO VAN and Ni. (b) High resolution XPS spectra of Ni 2p 3/2 for Ni/NiO VAN. (c) STEM images of Ni/NiO VAN, indicating the selected region for EELS measurement. Ni-L edge (d) and O-K edge (e) EELS of Ni/NiO VAN. (f) Schematic diagram of Ni/NiO heterostructure with oxygen vacancies at the interface.

19 Figure 3 .

 193 Figure 3. (a) LSV curves of NiO, Ni/NiO VAN and Ni electrodes for HER. (b) Overpotentials to reach a current density of 10 A cm -2 . (c) The Tafel slopes for HER. (d) LSV curves of NiO, Ni/NiO VAN and Ni electrodes for OER. (e) Overpotentials to reach a current density of 10 A cm -2 . (f) The Tafel slopes for OER.

Figure 4 . 21 Figure 5 .

 4215 Figure 4. (a) O 1s NAP-XPS spectra of NiO, Ni/NiO VAN under different water vapor pressures. (b) Plots of the percentage of adsorbed hydroxyl as a function of water vapor pressure. (c) Schematic illustration of the enhanced water dissociation mechanism by oxygen vacancies at the Ni/NiO interface.

Figure S1 .

 S1 Figure S1. The out-of-plane lattice parameters, in-plane lattice parameters and unit cell volumes of the NiO and Ni/NiO VAN extracted from RSMs.

Figure S2 .

 S2 Figure S2. High resolution XPS spectra of Ni 2p 3/2 for (a) Ni and (b) NiO.

Figure S3 .

 S3 Figure S3. High resolution XPS spectra of O 1s for (a) Ni/NiO VAN, (b) NiO, and (c) Ni.

Figure S4 .

 S4 Figure S4. Nyquist plots of the electrochemical impedance spectroscopy at the potential of -0.4 V vs RHE. (b) The equivalent circuit image, R s , R ct and C dl mean the electrolyte resistance, charge transfer resistance and electrical double-layer capacitance, respectively.

Figure S5 .

 S5 Figure S5. Stability curves measured at constant current density of -20 μA cm -2 in 0.1 M KOH.

Figure S6 .

 S6 Figure S6. Nyquist plots of the electrochemical impedance spectroscopy at the potential of 1.6 V vs RHE. The inset is the equivalent circuit.

Figure S7 .

 S7 Figure S7. Stability curves measured at constant current density of 20 μA cm -2 in 0.1 M KOH.

Figure S8 .

 S8 Figure S8. CV curves of (a) Ni/NiO VAN, (b) Ni, and (c) NiO at different scan rates from 10 to 60 mV s -1 . (d) Liner fitting of the capacitive currents of the catalysts versus scan rate.

Figure S9 .

 S9 Figure S9. LSV curves of water electrolysis with a two-electrode configuration in 0.1 M KOH.

Figure S10 .

 S10 Figure S10. O 1s NAP-XPS spectra of NiO, Ni/NiO VAN under different water vapor pressures.

Figure S11 .

 S11 Figure S11. Atomic structure of H 2 O adsorbed on: (a) Ni (001), (b) NiO (001), (c) metallic Ni sides near the interface in Ni/NiO (001) and (d) NiO sites near the interface in Ni/NiO (001), and (e) NiO sites near the interface in w/o-Ni/NiO (001).

Figure S12 .

 S12 Figure S12. Calculated energy barriers of water dissociation kinetics of pure Ni (001), pure NiO (001), and NiO sites near the interface in w/o-Ni/NiO (001).

http://chem.scichina.com/english

https://engine.scichina.com/doi/10.1007/s11426-022-1326-2

Acknowledgements

Kelvin H.L. Zhang is grateful for funding support by the National Natural Science Foundation of China (Grant No. 21872116). F. E. Oropeza and V. A. de la Peña O'Shea are grateful for the funding supported. by the EU (ERC CoG HyMAP 648319) and Spanish AEI (NyMPhA PID2019-106315RB-I00). Also, this work has been funded by the regional government of Comunidad de Madrid and European Structural Funds through their financial support to FotoArt-CM project (S2018/NMT-4367). Besides, Fundación Ramon Areces funded this work though ArtLeaf project. Kelvin H.L. Zhang also acknowledge the Sino-German Mobility Program (Grant No. M-0377). SuperSTEM is the National Research Facility for Advanced Electron Microscopy, funded from the Engineering and Physics Research Council (EPSRC). M.B. is grateful to the SuperSTEM Laboratory for microscope access, and to the School of Chemical and Process Engineering at the University of Leeds for a visiting associate professorship and financial support.

Science China Chemistry

A c c e p t e d https://engine.scichina.com/doi/10.1007/s11426-022-1326-2 Science China Chemistry A c c e p t e d https://engine.scichina.com/doi/10.1007/s11426-022-1326-2 Science China Chemistry A c c e p t e d https://engine.scichina.com/doi/10.1007/s11426-022-1326-2 Science China Chemistry A c c e p t e d https://engine.scichina.com/doi/10.1007/s11426-022-1326-2

Supporting information Vertically Aligned Ni/NiO Nanocomposites with Abundant Oxygen Deficient

Hetero-interfaces for Enhanced Overall Water Splitting

Hongxia Wang 1# , Meiyan Cui 2# , Gaoliang Fu

#Equal Contribution

Email: freddy.oropeza@imdea.org; kelvinzhang@xmu.edu.cn


Laurentiu Stoleriu

Masamichi Nishino

Seiji Miyashita

Alexandru Stancu

Roman Bertoni

Eric Collet

Maciej Lorenc

Cristian Enachescu

Multiscale out-of-equilibrium dynamics driven by pulsed laser excitation in spin-crossover materials: A combined thermoelastic and mechanoelastic study

In this paper, we present an elastic model coupled with a heat propagation process in order to reproduce the out-of-equilibrium dynamics of spin crossover materials driven by femtosecond laser excitation: a multiscale out-of-equilibrium dynamics driven by pulsed laser excitation in spin-crossover materials (thermoelastic step), the thermal switching (thermal step), and the subsequent relaxation to the initial state on cooling. The simulations were performed for open boundaries two and three-dimensional samples, composed of individual molecules linked by springs, which stand for elastic interactions. This building-up of the samples allows the propagation of elastic waves, which leads to accumulation of high spin molecules towards edges at the maximum of the thermoelastic step. We first show that a simple model with a single "temperature" reproduces the thermoelastic, the thermal step and the relaxation to the original equilibrium state. However, the too slow thermalization of the lattice obtained in this model does not correspond to the experimental data. Therefore, to overcome this drawback, we consider either an inhomogeneous photoexcitation or different "temperatures" for the lattice and the spin states. The effect of the sample size, which prevents the existence of a thermal step in the case of nanoparticles is also discussed, as well as the three-dimensional model.

I. INTRODUCTION

Spin crossover systems are metal complexes showing the unique ability to have two metastable spin states, high spin (HS) and low spin (LS), characterized by specific optical and magnetic properties, that can be triggered by using external stimuli, such as temperature, pressure, light [START_REF] Gütlich | Spin Crossover in Transition Metal Compounds[END_REF]. The smaller cohesion energy in the HS state determines longer coordination bond lengths and larger molecular volume. A first consequence of the change of molecular volume is the existence of elastic interactions between the metal ions, which if strong enough, leads to cooperative switching accompanied by a first order phase transition (hysteresis) with temperature, light or pressure as input parameters. This property is important for practical applications, as the information storage.

The thermal conversion from LS to HS state occurs at the critical temperature T C . The LS state, stable at low temperatures, can be switched towards the long-lifetime metastable HS state by optical excitation of some specific charge transfer or metal-centred bands [START_REF] Zerdane | [END_REF]3], with a light irradiation, by the way of the so called light induced excited spin state trapping (LIESST) effect [4,5]. Pulsed laser irradiation allows the control of this bistability in an ultrafast way, by induc-* cristian.enachescu@uaic.ro ing out-of-equilibrium conditions [6][7][8][9][10][11][12]. In solids, the initial 49 femtosecond photoswitching is localized at molecular level 50 for few absorbing molecules towards HS state and triggers 51 a subsequent increase of the fraction of HS molecules at 52 longer time, driven by lattice expansion. This phenomenon is 53 known as thermoelastic step. Recent development of elastic 54 models for spin-crossover materials, which treat the molecular 55 volume variation accompanied by the spin state change and 56 elastic interactions between molecules, has clarified important 57 cooperative properties in spin crossover phenomena [13][14][15][16][17][18][19][20]. 58 It was demonstrated that this cooperative switching occurs due 59 to the propagation of the elastic waves, taking place on the 60 acoustic time scale, which corresponds to the ratio between 61 relevant spatial scale and the sound velocity in the medium 62 [8]. The amplitude of the phenomenon is measured through 63 the HS fraction, which is the proportion of molecules in the 64 HS state, denoted here by X HS . Yet more spectacular, the 65 elastic step can be followed, for appropriate photoexcitation 66 rates and system sizes, by the so-called thermal step, which 67 is a further increase of X HS , at a larger timescale due to the 68 propagation of heat deposited by laser excitation throughout 69 the lattice [21].

In order to understand this complex out of equilibrium 71 behavior, where different degrees of freedom equilibrate over 72 different spatial and time scales, we should first clarify the role 73 of independent propagations of elastic waves and heat on the 74 recovered within 1 µs; in some cases a thermal step can be noticed, even in the absence of a thermoelastic step.

The increase of the high spin fraction subsequent to photoexcitation takes place in two steps, as we explain below. During ultrafast laser irradiation, a large amount of energy is locally deposited in every molecule photoexcited from LS to HS state, which causes both the local switching of electronic state and the molecular structural reorganization, together with the coherent and incoherent activation of molecular vibrational modes [23][24][25]. The LS-HS switching is accompanied by the increase of the molecular volume (metal-ligand distance) within less than 1 ps, which first determines the building-up and then the propagation of elastic interactions between spin crossover molecules within the crystalline lattice. This finally leads to the global lattice expansion, which produces negative lattice pressures on individual molecules, favoring the HS state of higher volume and resulting in a fast switching of more molecules towards HS state through a cooperative elastic process [26]. After this thermoelastic expansion (thermoelastic step), HS fraction decreases for a while. On the other hand, the energy deposited into photoexcited molecules, which far exceeds the amount needed for the LS-HS conversion (20 meV), is transferred to the surrounding lattice through vibrational cooling [START_REF] Von Allmen | Laser-Beam Interactions with Materials[END_REF]. The distance between "hot spots (photoexcited sites)" is small for homogeneous excitation: namely, for 4% excitation one molecule over 25 absorb energy, which corresponds to one hot spot in a box of 3 × 3 × 3 molecules. Therefore the energy redistribution leads to a rapid homogenization (on the order of a dozen of ps) of the temperature of the lattice through molecular vibrations and phonon population [START_REF] Ruello | [END_REF]29]. This lattice warming can also result in its thermal expansion. However, it should be noted that the high spin state population becomes thermally equilibrated with the lattice in a time of the order of 10 µs, which corresponds to the "thermal step" in Fig. 1, much slower than the thermalization of the lattice. Thus the difference of the time scale of the thermalizations of lattice and spin state plays an important role for the elastic step and the thermal step.

B. Models and discussion

The thermo-mechano-elastic model

In order to discuss both elastic and thermal steps in the framework of a unified model, we use here a modified mechanoelastic model, which was previously applied to simulate the elastic self-amplification alone. In this model, the molecules, represented as rigid spheres, are situated in open boundary lattices. The elastic interactions are simulated by springs linking a molecule to its closest neighbors. When a molecule changes its state, its volume varies, which results in an immediate elongation or a compression of its closest springs, determining first the motion of neighboring molecules and then the propagation of the initial perturbation towards all the lattice. The elastic simulations of spin crossover materials imply two processes: the switch of spins and the change in molecular positions [30]. The spin change may be performed either by the way of Metropolis criterion or by an Arrhenius approach, which is used here, as faster and more appropriate to dynamic phenomena.

P i HS→LS = 1 τ exp - E A -κ p i k B T i , P i LS→HS = 1 τ exp - D -k B T i ln g k B T i × exp - E A + κ p i k B T i , (1)
where τ is a constant scaling factor ensuring that the probabil- as in the case of fast phenomena subsequent to femtoseconds photoexcitation experiments [30].

197 m d 2 x i dt 2 = F i,x -µ dx i dt m d 2 y i dt 2 = F i,y -µ dy i dt , (2)
The steady state of the system before photoexcitation is generated by computing at least 5000 MCES at the given temperature (i.e., at 145 K where most of the simulations have been performed to reach the steady state corresponding to a X HS of around 0.06). The material parameters used in simulations, in line with standard experimental calorimetric data for Fe(II) spin crossover systems [32,33], were D = 1100 K, g = 1096, thus giving a thermal transition temperature [D/(k B ln(g))] of around 157 K. The radius of HS molecules is considered to be 0.22 nm, and that of LS molecules 0.2 nm. The distance between centres of molecules is 1.04 nm for two molecules in HS state and 1 nm for two molecules in LS state. These values correspond to x-ray experimental measurements for typical spin crossover compounds [32,33]. The elastic constant of intermolecular springs is 7 N/m, which generates a moderate cooperativity in the system [34]. As in previous works, κ is taken 1450 × 10 -14 J/N while τ = 1000 s -1 [11,13].

Photoirradiation randomly transforms more molecules to the HS state in addition to those in the steady state. The temperature of the photoexcited molecules is increased by T = 100 K. The simulation starts immediately after the first thermalization and deals with the second and third thermalization, as described in the introduction. The heat propagates through the sample following the Fourier's law of heat conduction:

dQ dt = -k∇T (3)
with k denoting the thermal conductivity. Using the continuous equation of the heat and the relation Q = ρc p T , the equation of the distribution of the temperature is written after a few transformations [START_REF] Sands | Pulsed laser heating and melting[END_REF] as

dT dt = k ρc p ∇ 2 T = D T ∇ 2 T, ( 4
)
where ρ is the mass density, c p the specific heat at constant pressure and D T denotes the thermal diffusion constant.

Equation ( 4) can be approximately transposed by a finite difference method [START_REF] Sands | Pulsed laser heating and melting[END_REF][START_REF] Shneidman | [END_REF], considering two kinds of heat transfer, i.e., the diffusion in the bulk of system and the heat transfer to the bath (external), to the following equation of lattice temperature diffusion:

dT i L dt = -α T i L -T i L -β T i L -T B , (5)
where T i L is the lattice temperature of the ith molecule, ⟨T i L ⟩ is the average of lattice temperatures for all neighbors of the ith molecule, α is the diffusion coefficient, T B is the external thermal bath temperature and β is the heat transfer coefficient to the bath. By the second term at the right-hand side, which is considered for only edge molecules, the lattice temperature approaches the bath temperature from the edges to the inner part.

The algorithm for simulations is then completed within the thermoelastic model with the computation of temperatures of every molecule considering Eq. ( 5), after every update of all the spin state MCES after photoexcitation. Typical and 2(e)]. However, this slow thermalization over the whole lattice is not compatible to the experimental fact that the energy is distributed very rapidly to the whole lattice uniformly as mentioned above. This problem will be resolved in the following. Finally, due to heat transfer to the external bath, the temperature in the system slowly decreases to the initial one and most of the molecules switch back to their ground LS state.

As we have discussed in the previous paragraph, all the macroscopic phenomena are successfully reproduced within the present model. However, in the introductory part, we have claimed that the thermalization of the lattice is very fast, and happens long time before the thermal step, as discussed even in early works concerning the photoinduced molecular switching [37]. Within this aspect, the results presented in 300 Fig. 2(c) indicates a weakness of the present approach: we 301 notice that at the maximum of the thermoelastic step, the 302 distribution of the temperature in the lattice is very large; ac-303 tually, the thermalization of the lattice is produced in the same 304 time as the thermal step. In order to avoid this drawback, we 305 propose in the following two alternatives of the current model 306 corresponding to different thermalization processes of differ-307 ent subsystems. In previous works on ultrafast magnetism, 308 the different thermalization timescales and energy transfert 309 between electron, spin and lattice have been described through 310 different temporal evolution of their respective temperatures. 311 It was shown that they are equilibrated in the irradiated 312 material after a very short time (about 5 ps) [38][39][40][41] the 313 photoinduced out-of-equilibrium dynamics in Bismuth was 314 discussed in terms of equilibration of electronic and lattice 315 subsystems [42,43]. In the present work, we use an analogy to 316 describe multiscale thermal equilibration.

317

The thermo-mechano-elastic model including a two

318 thermalization sequence 319

The multiscale out-of-equilibrium dynamics induced by a 320 femtosecond laser pulse in spin crossover materials and the 321 subsequent thermalizations of the different degrees of free-322 dom are different from what occurs in ultrafast magnetism or 323 coherent phonons, for which the energy is deposited on delo-324 calized electrons which thermalize rapidly as heated up by the 325 laser pulse. Indeed, in the case of spin crossover materials, 326 the energy is locally deposited on switched molecules that 327 absorbed photons. We therefore consider different subsystems 328 -the "hot" molecules locally photoswitched from low to high 329 spin state due to optical excitation, the lattice heated up by 330 energy transfert and expanding due to molecular swelling and 331 the spin state of the molecules on the lattice, which equilibrate 332 on a different timescale. (1) Ultrafast photoexcitation: the light transforms some of the LS molecules into hot photoexcited molecules with T * S .

(2) Photoexcited molecules release their energy towards neighboring lattice: first thermalization within t ML ≈ 10 ps [7]. After this stage the molecular spin state temperature of photoexcited molecules and the lattice temperature are equal:

T S (t ML ) = T L (t ML ) = T B + T , while for the other molecules (3) Lattice towards spin state: the energy transfert from lattice to molecule is responsible for a slower spin state switching, which is responsible for slower thermal population of the HS state within t SS ≈ 10 µs. In this stage, the spin temperature becomes close to the lattice temperature, and the conversion from LS to HS at nonphotoexcited molecules takes place. After this step, the global thermalization of the subsystems is reached with equal temperatures [44].

T L (t ML ) = T B + T > T S (t ML ),
(4) System towards bath: third thermalization within t B ≈ 1 µs-1 ms. In this step, the energy injected by photoirradiation is released to the bath. This step ensures that the "temperatures" of the ensemble reaches the temperature of the bath (this time varies as a function of system size, and the process takes place from the border to the inner places of the sample), so that the system recovers initial equilibrium state prior to photoexcitation. S is the temperature of converted molecules just after photoexcitation; T L , lattice temperature, corresponds to all vibrational modes, T C is the critical temperature. Before excitation, all subsystems have the same temperature, which is the temperature of the bath T B . (Bottom) Simulations using the mechanoelastic model for the 2D system and correspondence between the thermalizations and the elastic, thermal step and relaxation towards initial state.

To be in accord with these statements, we need to extend 394 Eq. ( 5) to describe the equilibration of the spin state through 395 the equilibration of its temperature. Thus we introduce a 396 model in which Eq. ( 5) is phenomenologically coupled to the 397 equation of temperature diffusion from the lattice towards the 398 breathing (spin state) modes:

399 dT i S dt = -γ T i S -T i L , (6)
where γ is the spin-lattice interaction constant.

400

Using this model with appropriate values for the ther-401 mal coefficients (α = 0.005 MCS -1 , β = 0.005 MCS -1 , γ = 402 510 -5 MCS -1 ), we have simulated the evolution of the system 403 presented in Fig. 3 (bottom). Let us analyze the evolution of 404 the system. Only a small spin state relaxation can be observed 405 just after photoexcitation (b mark in Fig. 2); due to the fact 406 that most of the photoexcited molecules are found in a high-407 pressure environment, some of the HS molecules switch back 408 to the smaller volume LS state. However, the increase of the 409 temperature of the photoexcited molecule due to the laser 410 pulse prevents the switching back of more molecules. The 411 subsequent expansion of the volume of the whole lattice will 412 decrease the pressure on all molecules in the system, therefore 413 a large part of molecules, especially those on edge or corners, 414 004300-5

444

After this high-pressure state, the lattice expands and some 445 low spins converted to high spin state in the low pressure 446 due to the expansion. Then, the reorganization of the lattice takes place, forming small HS molecules clusters at edges or corners, corresponding to the thermoelastic step as depicted in Fig. 4(c). After this point, the spin temperature increases and many low spins begin to be converted to the high spin state. At the beginning of the thermal step [Fig. 4(d)] the clusters formed during previous steps develop from corners and expand towards the center of the lattice. Figure 4(e) shows the configuration at the maximum of thermal step; at this stage we notice that the system starts to cool down from the edges to central part due to heat flow towards the environment. Finally, at the end of the relaxation, the temperature of all molecules came back to their initial temperatures, and only few molecules thermally populate the HS state at T B , below T C . In the system studied here, a local spin state is coupled to lattice and patterns of HS and LS molecules appear. Similarly, in a ferroelastic system [45], a coupling between strain and electrical local dipole moments leads to a pattern formation about polarity. Moiré-type patterns on the arrangements of HS and LS systems, have been recently observed on thin layers of spin crossover systems on substrates in appropriate elastic models [46].

The dependence of the thermal coefficients on the thermal step are presented in Fig. 5. As we explained before, because the temperature variation during the thermoelastic step is small, these coefficients do not influence the amplitude or the position of the thermoelastic step.

The role of the internal coefficient γ is explained in Fig. 5 (top). When γ is large, the relaxation of T S to T L is fast, and thus the displacement of the thermal step appears at smaller values of the time. If the γ value is extremely large, the thermal step merges with the thermoelastic step. Reversely, a 004300-6 No elastic step and no thermal step take place (if the amount of energy transferred to the lattice is low). No elastic step, but thermal step takes place (for a higher amount of energy inside the lattice). (Bottom) Size effects: the thermal step is observed only in the case of a larger system, for which the exchange with the thermal bath is slower. on photoexcited molecules is high enough, then the temperature of the whole lattice increases enough in order to allow the presence of the thermal step, as experimentally shown in Fig. 1 for the smaller laser power. A trivial case of neither thermoelastic step nor thermal steps may be also obtained in experiments. The effect of a small size of the crystal which evacuate too fast the heat was also discussed in Ref. [12].

In Fig. 6 (bottom), we present the size effects, simulating the behavior of two systems with different sizes (13 824 and 3456 molecules), while keeping all other parameters in the system the same. In a smaller sample the heat escape to exchange with the bath is faster, and thus, in small samples, the temperature will not surpass the critical temperature. This corresponds to experimental data shown in Fig. 1, where the thermal step is observed only in the case of crystals, and there is no thermal step in nanoparticles. In addition, we notice that in the case of the small sample the height of the elastic step has the same value as in the case of the large sample, but it shifts towards shorter times. This is due to the fact that in a smaller sample the elastic wave reaches sooner the borders of the sample. For the same reason, the relaxation after photoexcitation is less intense. Simulations using the mechanoelastic model for the 3D system for two initial photoexcitation rates.

3D simulations of the two thermalization scales model 522

The 2D models can be considered to give a reasonable 523 approximation for conceptual understanding of the phenom-524 ena, especially for layered spin crossover compounds with 525 small interactions between the layers. However, several spin 526 crossover systems present a three-dimensional structure, with 527 the strength of interactions between planes on the same order as within planes and therefore it is important to extend the model towards 3D systems. Due to difficulties related to complexity of the system and the increased number of molecules, less studies have been devoted to modeling 3D systems except some works with open boundaries [47][48][49] or on surfaces [50][51][52]. In this section we use a rectangular cuboid system composed of 11 layers of molecules in a face-centred-cubic configuration. Each layer is composed of 1900 molecules in a triangular configuration; every bulk molecule is linked to its twelve closest neighboring molecules (six on the same plane, three below and three above) by springs; molecules on surface layers and those situated on the edge have less neighbors. The probabilities and the thermal diffusion equations are similar as for the 2D case; we have to note that in the 3D systems the molecules are allowed to move outside their initial planes to produce more favourable energetic conditions, leading to so-called buckling effects [47], which have been treated in previous papers for 3D systems of various shapes [48] or for monolayers of spin crossover molecules on substrates, leading to moiré patterns [46]. In Fig. 7, we present the whole curves obtained after 15% and 20% percentage photoexcitation, using the following transfer coefficient parameters α = 5 × 10 -7 MCES -1 , β = 5 × 10 -7 MCES -1 , and γ = 10 -6 MCES -1 . As in the case of 2D systems we notice the well-defined presence of both the thermoelastic and thermal step. The amplitude of the thermoelastic step is larger in the case of higher initial photoexcitation, as previously observed in experiments [7] and simulations [10].

We refer now to the configuration of the system during the thermoelastic step. As specified above, in a 2D system a reorganization of the molecules in the system is observed at the maximum of the elastic step, consisting of the accumulation of HS molecules towards edges and corners. We analyze here the molecular configuration at the maximum of the elastic step for a 3D configuration. In Fig. 8(a), we present a snapshot of the whole system after 30000 MCS, while the snapshots of the first and middle layers are presented in Fig. 8(c

The case of inhomogenous photoexcitation 575

In this section, we refer to the simple thermo-mechano-576 elastic model, including only the lattice temperature, but we 577 consider an inhomogenous photoexcitation, which can be ex-578 pected in the case of sample with size of the order or larger 579 than the light penetration depth. The photoexcitation probabil-580 ity can be written as P = exp(-εx), where ε is the absorption 581 coefficient and x is the distance of a molecule from the front 582 of irradiation. The effect of the absorption coefficient on 583 the distribution of the photoexcited molecules is depicted in 584 the right panel of the Fig. 9. In the same figure, we present 585 the evolution of the system keeping constant the heat transfer 586 parameters and varying the absorption coefficient. We find 587 that the homogenous case (ε = 0) gives a similar dependence 588 correlated regions, as it was theoretical stated in a review on the crackling noise in crystals, ferroelastic and porous materials [START_REF] Salje | Annual Review of Condensed Matter Physics[END_REF]. Within this aspect, it would be interesting to determine a possible power law function describing the current multiscale phenomena. Similar approaches can be applied for other molecular systems with fast temperature variations or for spin crossover systems heated by plasmonic nanodevices [START_REF] Hu | [END_REF]. It should be also noted that in the current method, the spin-transition is done by Monte Carlo sampling. Here we assume that the spin transition itself is much faster than the time scale of lattice motion. The spin transitions occur by contact with the thermal bath. Thus the comparison between the time scale of lattice motion and the that of contact with the bath. In principle, the ratio may result in different aspects of dynamics. This problem was studied by Nishino et al. [30], but so far, in our recent works, we used the present scheme to catch general aspects of dynamical phenomena. In the Monte Carlo simulations, the effect of degeneracies of HS and LS is treated in a thermodynamic approach, but not in a dynamical way. This must be investigated more carefully in the future. In principle, we may study the spin transition in the relation of lattice dynamics [55], but several fundamental problems remain to be studied. Thus, in particular for studies of dynamical properties, we should keep in mind these problems. But we still believe that the present work captures many important aspects of the system.

52 LAURENTIU

 52 FIG. 2. Evolution of the system after ultrafast photoexcitation showing relaxation, thermoelastic step, thermal step and the final relaxation to the ground state. (a) Snapshots of the system at the maximum of the thermoelastic step [(b) and (c)] and at the maximum of the thermal step [(d) and (e)]. Spin states maps [(b) and (d)]: red circles: HS molecules, yellow circles: LS molecules; down: temperature map). Temperature maps [(c) and (e)] blue: low temperature, red: high temperature

  photoexcited molecules and T B the temperature of the external thermal bath (cryostat, polymer, etc.) which corresponds to the initial (prior to photoexcitation) and final (after the relaxation of the thermal step) temperatures. The energy transfers and "temperatures" (we use quotes as the concept of temperature, necessary in the model, is different at the level of a single molecule) can be summarized as follows [see Fig. 3 (top)].

  photoexcited molecules and E is the energy injected by a photon.

FIG. 3 .

 3 FIG. 3. (Top)Schema for the behavior of the system after ultrafast photoexcitation: T S -like, molecular (spin state) temperature, corresponds only to breathing modes; T * S is the temperature of converted molecules just after photoexcitation; T L , lattice temperature, corresponds to all vibrational modes, T C is the critical temperature. Before excitation, all subsystems have the same temperature, which is the temperature of the bath T B . (Bottom) Simulations using the mechanoelastic model for the 2D system and correspondence between the thermalizations and the elastic, thermal step and relaxation towards initial state.

FIG. 4 .

 4 FIG.4. Snapshots of a 2D system just after photoexcitation (a), during the small relaxation following the photoexcitation (b), at the maximum of the elastic step (c), at the beginning of the thermal step (d), at the maximum of the thermal step (e) and at the end of the final relaxation (f). Yellow: LS molecules, Blue, red, orange: HS molecules (from low to high temperatures, the scale is not the same for all figures). The letters correspond to those marked in Fig.3, bottom.

  FIG. 5. Dependence of time evolution of X HS on coefficients (top) γ , (middle) α, and (bottom) β.

52 LAURENTIU

 52 FIG.7. Simulations using the mechanoelastic model for the 3D system for two initial photoexcitation rates.

FIG. 8 .

 8 FIG. 8. A snapshot of 3D system at the maximum of the elastic step (a), snapshots of the first (c) and middle layer (d). The dependence of X HS on the position of layers are shown in (b). 004300-8

  ) and respectively Fig.8(d). At first sight, we notice a higher proportion of HS molecules in the first layer comparing to the middle layer. A quantitative analysis showing the X HS for every layer is presented in Fig.8(b) and confirms the initial observation: there are more HS molecules in the layers near the surface and their number decreases in the layers situated towards the center of the sample. Therefore we can conclude that in 3D systems the molecular reorganization is similar to the one 573 visible in 2D systems.

  574

FIG. 10

 10 FIG. 10. (Left) Evolution of the X HS after photoexcitation for ε = 0.088 (a). Snapshots of the system at the maximum of the thermoelastic step [(b) and (c)] and at the maximum of the thermal step [(d) and (e)], spin states map [(b) and (d)], and temperature map [(c) and (e)], using the same scale as in Fig. 2.

  Because we study local spring interactions

	BY14107	PRB	July 10, 2023	13:52
		MULTISCALE OUT-OF-EQUILIBRIUM DYNAMICS DRIVEN …	PHYSICAL REVIEW B 00, 004300 (2023)
	166	among molecules, the strain effect (elastic interaction) is au-
	167	tomatically taken into account,
	168	Therefore every spin change will influence both neigh-
	169	boring and, in a smaller extent, far-away molecules. In this
	170	way, a single elastic constant stands for both short-range and
	171	long-range interactions.
	172	Most of the simulations in this paper have been per-
	173	formed for a 2D rectangular shape system composed of 13 824
	174	molecules in a triangular configurational bonding. The Arrhe-
	175	nius molecular switching probabilities of the ith spin between
	176	the LS and HS states depend on the "temperature" T , on
	177	intrinsic material parameters (the HS-LS energy difference D,
	178	the degeneracy ratio g, the effective activation energy E A ) and
	179	on the effective elastic interactions between molecules, repre-
	180	sented here by the way of the local pressure force p i acting
	181	on molecule i, and are explicitly given by in the following
	182	equations:
				004300-2

  where x i and y i are the Cartesian coordinates of the molecule

	198
	199	i, m is the mass of the molecule, µ is a damping constant,
	200	preventing the system to enter into an uncontrolled oscillatory
	201 202	motion, and F xi , F yi are the components of the instantaneous force -→ F i acting on molecule i given by the sum of the forces
	203	from the neighboring springs.
	204	The number of steps r to solve the system of coupled
	205	differential equations after every MCES is a key ingredient
	206	for the transient evolution of the system [26,31]. The value of
	207	r, tunes the different time scales of spin dynamics and lattice
	208	relaxation and is a measure for how fast the lattice relaxation
	209	is in comparison with individual molecular switching. A large
		r favours equilibrium distribution, specific for static phenom-

210

ena, whereas a smaller r favours nonequilibrium distribution 211

ACKNOWLEDGMENTS

The authors thank professor Hervé Cailleau for fruitful discussions during the manuscript preparation. This work was supported by a grant of the Romanian Ministry of Research, Innovation and Digitization, CNCS/CCCDI-UEFISCDI, FQ Project No. PN-III-P4-ID-PCE-2020-1946, within PNCDI III. The authors gratefully acknowledge Agence Nationale de la Recherche for financial support undergrant, ANR-19-CE30-0004 ELECTROPHONE, ANR-19-CE29-0018 MUL-TICROSS. E.C. thanks the University Rennes 1 and the Fondation Rennes1 for funding. The present work was supported by the Elements Strategy Initiative Center for Magnetic Materials (ESICMM) (Grant No. 12016013) funded by the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan, and was partially supported by Grants-in-Aid for Scientific Research C (Grant No. 18K03444 and No. 20K03809) from MEXT.

to that presented in Fig. 3. If the absorption is present (ε ̸ = 0), 589 then a clearer disentangling between thermoelastic and elastic 590 steps begins to be visible. Let us refer now to the case when to experimental data. Actually, in some previous paper [7,10],

the height of the simulated thermoelastic was smaller than the 597 experimental one, so an inhomogenous photoexcitation could 598 approach the simulated data to the experimental ones.

599

In the right panel of Fig. 10, we present snapshots of the