# Fileset

[nanomaterials14(2024)01775.pdf](https://mdr.nims.go.jp/filesets/91f490ef-b2ca-4015-a114-03fb66fbc969/download)

## Creator

[Aaditya Manjanath](https://orcid.org/0000-0003-2223-5954), [Ryoji Sahara](https://orcid.org/0000-0003-0788-2985), Yoshiyuki Kawazoe, Kaoru Ohno

## Rights

[Creative Commons BY Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0/)

## Other metadata

[Non-Adiabatic Excited-State Time-Dependent GW (TDGW) Molecular Dynamics Simulation of Nickel-Atom Aided Photolysis of Methane to Produce a Hydrogen Molecule](https://mdr.nims.go.jp/datasets/b853a986-2618-4875-b225-f165fb68e98b)

## Fulltext

Non-Adiabatic Excited-State Time-Dependent GW (TDGW) Molecular Dynamics Simulation of Nickel-Atom Aided Photolysis of Methane to Produce a Hydrogen MoleculeCitation: Manjanath, A.; Sahara, R.;Kawazoe, Y.; Ohno, K. Non-AdiabaticExcited-State Time-Dependent GW(TDGW) Molecular DynamicsSimulation of Nickel-Atom AidedPhotolysis of Methane to Produce aHydrogen Molecule. Nanomaterials2024, 14, 1775. https://doi.org/10.3390/nano14221775Academic Editors: Luis Alberto RuizPestana, Wenjie Xia and Zhaoxu MengReceived: 30 September 2024Revised: 30 October 2024Accepted: 2 November 2024Published: 5 November 2024Copyright: © 2024 by the authors.Licensee MDPI, Basel, Switzerland.This article is an open access articledistributed under the terms andconditions of the Creative CommonsAttribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).ArticleNon-Adiabatic Excited-State Time-Dependent GW (TDGW)Molecular Dynamics Simulation of Nickel-Atom AidedPhotolysis of Methane to Produce a Hydrogen MoleculeAaditya Manjanath 1 , Ryoji Sahara 1 , Yoshiyuki Kawazoe 2,3,4 and Kaoru Ohno 1,5,*1 Research Center for Structural Materials, National Institute for Materials Science (NIMS), 1-2-1 Sengen,Ibaraki, Tsukuba 305-0047, Japan; manjanath.aaditya@nims.go.jp (A.M.); sahara.ryoji@nims.go.jp (R.S.)2 New Industry Creation Hatchery Center, Tohoku University, 6-6-4 Aramaki Aza Aoba, Aoba-ku, Miyagi,Sendai 980-8579, Japan; kawazoe@e-workshop.co.jp3 School of Physics, Institute of Science, Suranaree University of Technology, 111 University Avenue,Nakhon Ratchasima 30000, Thailand4 Physics and Nanotechnoloy, SRM Institute of Science and Technology, Tamil Nadu,Kattankurathur 603203, India5 Department of Physics, Graduate School of Engineering, Yokohama National University (YNU),79-5 Tokiwadai, Hodogaya-ku, Yokohama 240-8501, Japan* Correspondence: ohno@ynu.ac.jpAbstract: Methane photolysis is a very important initiation reaction from the perspective of hydrogenproduction for alternative energy applications. In our recent work, we demonstrated using our re-cently developed novel method, non-adiabatic excited-state time-dependent GW (TDGW) moleculardynamics (MD), how the decomposition reaction of methane into a methyl radical and a hydrogenatom was captured accurately via the time-tracing of all quasiparticle levels. However, this processrequires a large amount of photoabsorption energy (PAE ∼10.2 eV). Moreover, only one hydrogenatom is produced via a single photon absorption. Transition metal atoms can be used as agents forphotochemical reactions, to reduce this optical gap and facilitate an easier pathway for hydrogenproduction. Here, we explore the photolysis of methane in the presence of a Ni atom by employingTDGW-MD. We show two possibilities for hydrogen-atom ejection with respect to the location of theNi atom, towards the Ni side or away from it. We demonstrate that only the H ejection away from theNi side facilitates the formation of a hydrogen molecule with the quasiparticle level correspondingto it having an energy close to the negative ionization potential of an isolated H2 molecule. This isachieved at a PAE of 8.4 eV which is lower compared to that of pristine methane. The results obtainedin this work are an encouraging step towards transition metal-mediated hydrogen production viaphotolysis of hydrocarbons.Keywords: quasiparticle; Koopmans’ theorem; TDDFT; GW approximation; Ni atom; CH4; chemicalreaction; photochemistry; surface hopping1. IntroductionMethane is a very important fuel gas as (a) it is a main constituent of liquefied naturalgas (LNG) [1], which is useful for long-distance transport, (b) it is the main component ofbiofuel or biogas [2], making it a source of clean energy, and (c) it is useful as a direct coolantin jet engines. In addition, methane serves as a precursor gas for hydrogen production. H2is considered the “fuel of the future” because it produces three times the amount of energy(39.4 kWh kg−1) compared to other fuels such as liquid hydrocarbons (13.1 kWh kg−1) [3].It has been reported extensively that the endothermic decomposition of methane leadsto the production of hydrogen. Cracking [4] (heating of methane in the absence of air),photolysis [5–9] (photoexcitation of methane), and steam reforming [10,11] (reaction ofmethane with steam), as well as thermocatalytic decomposition [12–16] and solar-aidedNanomaterials 2024, 14, 1775. https://doi.org/10.3390/nano14221775 https://www.mdpi.com/journal/nanomaterialshttps://doi.org/10.3390/nano14221775https://doi.org/10.3390/nano14221775https://creativecommons.org/https://creativecommons.org/licenses/by/4.0/https://creativecommons.org/licenses/by/4.0/https://www.mdpi.com/journal/nanomaterialshttps://www.mdpi.comhttps://orcid.org/0000-0003-2223-5954https://orcid.org/0000-0003-0788-2985https://orcid.org/0000-0002-2369-7045https://orcid.org/0000-0002-1980-5971https://doi.org/10.3390/nano14221775https://www.mdpi.com/journal/nanomaterialshttps://www.mdpi.com/article/10.3390/nano14221775?type=check_update&version=1Nanomaterials 2024, 14, 1775 2 of 18decomposition [17–20], both of which produce hydrogen along with “carbon black” ornano-sized carbon clusters, without emitting greenhouse gases, are the most commonlyemployed processes for this.From the perspective of hydrogen production, it is interesting and important to investi-gate the photochemical reactions of hydrocarbons. The computational approach is more fa-vorable than the typical experimental approach when investigating ultrafast reactions suchas photolysis. However, density functional theory (DFT) [21] is, in principle, a ground-statetheory and cannot be applied to any photoexcited state. Time-dependent density functionaltheory (TDDFT) [22] relying on adiabatic local density approximation (ALDA) [23] hasbeen the standard approach to investigate the excited-state (ES) dynamics; however, itfaces a longstanding difficulty of ALDA not being applicable to the time-dependent (TD)molecular dynamics (MD) for an initially excited state (ES). We have recently overcome thisdifficulty [24] by employing extended quasiparticle theory (EQPT) [25,26], in which the fullcorrespondence is achieved between the ES surfaces and corresponding total energies, withsatisfying extended Koopmans’ theorem [27,28]. In EQPT, each quasiparticle (QP) levelis related to a total energy difference. The QP energies of an occupied level (εQPi ) and anunoccupied level (εQPa ) are defined, respectively, asεQPi = E(N)ref − E(N−1)i→vac , (1a)εQPa = E(N+1)vac→a − E(N)ref , (1b)where E(N)ref , E(N−1)i→vac , and E(N+1)vac→a are the total energies of the reference N-electron system,the (N − 1)-electron system formed by removing an electron from the ith occupied levelto the vacuum level (vac), and the (N + 1)-electron system formed by adding an electronat the ath unoccupied level from vac. Here, we emphasize that the reference N-electronsystem is not necessarily the ground state (GS), but can be any of the excited eigenstates.The QP energies εQPi and εQPa in Equation (1) can be referred to as ‘negative ionizationpotentials’ and ‘negative electron affinities’, respectively, and have a direct correspondencewith the observations from photoemission/inverse photoemission experiments. The GWapproximation (GWA) is in full conformity with EQPT. We have applied EQPT within theGWA to the mixed quantum-classical Ehrenfest dynamics with surface hopping (SH) anddeveloped the NA-ES-TDGW (TDGW) MD method [24]. The merit of this method is thatwe can trace all the QP energies as well as the QP wavefunctions, which we will refer toas “charge densities”, during the simulation. Using this method, we have successfullyinvestigated methane photolysis CH4h̄ω−−→ CH ·3 + H at the lowest photoexcited state [24],which is considered as the initiation reaction of a complex multi-step process of a variety ofmethane combustion and hydrogen production reactions.Nevertheless, the photolysis of methane requires a large photoexcitation energy. Inthis regard, transition metal atoms can pave the way for reducing such large optical gapsand facilitate an easier pathway for hydrogen production. Here, we focus on the effect of atransition atom in this reaction. The aim of the present study is to unravel new reactionpathways in nickel atom-mediated methane photolysis. To the best of our knowledge, therehas been no direct molecular dynamics study to search the reaction pathway of CH4 inthe presence of a transition metal atom at any of the photoexcited states, although therehave been several computational studies involving systems with transition metal atoms.One such study is the TDDFT molecular dynamics investigation of the hydrogen spilloverprocess via Ni dimers by Sahara et al. [29]. They showed that a hydrogen molecule canbe dissociated into two hydrogen atoms near the Ni dimer by a small excitation energy.Another is a DFT study on the potential energy surfaces of a CH4 molecule with a Ni atomby Burghgraef et al. [30] or with an Fe atom by Sun et al. [31]. From their results, reactionssuch as CH4 + M −−→ CH3MH (M = Ni or Fe) can be expected. Concerning the late-stagedynamics of the thermocatalytic decomposition of methane, there is a DFT study on thepotential energy surface of a Ni55 cluster with an open (semi-capped) carbon nanotube byNanomaterials 2024, 14, 1775 3 of 18Wang et al. [32]. Although the earlier stage of methane decomposition was not consideredin their work, the mechanism of the catalytic behavior of the Ni cluster was clarified.The presence of a Ni atom breaks the tetrahedral symmetry of CH4, rendering a non-unique choice of the photoexcited state. The highest occupied molecular orbital (HOMO) ofthis system is no more that of the methane fragment but is now one of the highest occupiedNi 3d levels. Moreover, the lowest unoccupied molecular orbital (LUMO) of this system isno more than the negative electron affinity level of pure CH4 above the vacuum level butis the lowest unoccupied Ni 3d level, which is below the vacuum level. This means thatthere are several possible reaction pathways depending on the choice of the photoexcitedexcited state, for which the excitation energy is lower than 10.2 eV [5–7] that is needed inthe photolysis of the pristine methane case. In the course of our study, we show that thereare at least two possible pathways of methane decomposition. In one case, hydrogen atomson the side opposite Ni are ejected from the methane molecule, while in the other case,hydrogen atoms facing the Ni atom are ejected from the methane molecule. We present theresults of the two different pathways in detail in Section 3.2. Method2.1. Time-Dependent QP EquationThe wavepackets ϕλ(r; R(t), t), where r, R(t) = {RI(t)}, and t represent the positionof the QP, a set of nuclear coordinates (I - atomic index), and time, respectively, satisfy thetime-dependent QP (TDQP) equation [24]i∂∂tϕλ(r; R(t), t) = HQPmixed(r; R(t))ϕλ(r; R(t), t), (2)which is similar to the time-dependent Kohn–Sham equation [22]. Here, HQPmixed(r; R(t)) de-notes the QP (GW) Hamiltonian for the mixed state constructed by these wavepackets [24].We introduce the QP wavefunctions φQPn (r; R(t)) and the QP energies εQPn (R(t)) [n = i(occupied) or a (unoccupied), QP level indices], which satisfy the eigenvalue equationHQPmixed(r; R(t))φQPn (r; R(t)) = εQPn (R(t))φQPn (r; R(t)) (3)during the simulation. Then, using the spectral method [24,29,33], we can expand thewavepackets ϕλ(r; R(t), t) asϕλ(r; R(t), t) = ∑nφQPn (r; R(t))cnλ(t) (4)with the expansion coefficientscnλ(t) = ⟨ φQPn (R(t)) | ϕλ(R(t), t) ⟩. (5)The TDQP Equation (2) can be numerically solved by propagating the wavepackets insmall timesteps ∆t. The wavepackets ϕλ(r; R(t + ∆t), t + ∆t) at a slightly later time t + ∆twith a small time interval ∆t can be expressed asϕλ(r; R(t + ∆t), t + ∆t) = exp[−i∫ t+∆ttHQPmixed(r; R(t))dt]ϕλ(r; R(t), t)≈ ∑nexp[−iεQPn (R(t + ∆t))]φQPn (r; R(t + ∆t))cnλ(t + ∆t) (6)Nanomaterials 2024, 14, 1775 4 of 18withcnλ(t + ∆t) = ⟨ φQPn (R(t + ∆t)) | ϕλ(R(t), t) ⟩≈ ⟨ φQPn (R(t)) | ϕλ(R(t), t) ⟩+ ∆t∂∂t′⟨ φQPn (R(t′)) | ϕλ(R(t), t) ⟩t′=t= cnλ(t)− ∆t ∑mcmλ(t)∑IṘI(t)〈φQPn (R(t))∣∣∇RI∣∣ φQPm (R(t))〉. (7)In the last equality of Equation (7), we used Equation (4) again. The second termof the rightmost expression in Equation (7) is proportional to the nuclear velocityṘI(t) = dRI(t)/dt and represents a non-adiabatic interaction [34,35]. The matrix ele-ments dnm =〈φQPn (R(t))∣∣∇RI∣∣ φQPm (R(t))〉represent the non-adiabatic coupling vectors.2.2. One-Shot GW Within TDQPWe apply this formalism to the one-shot GW approach [36], which is the simplest andmost reliable approach in the GWA. In the one-shot GW approach, the QP wavefunctionsare replaced by Kohn–Sham orbitals in the local density approximation (LDA) [36]φQPn (r; R(t)) ≈ φLDAn (r; R(t)), (8)and the QP energy eigenvalues φQPn (r; R(t)) are calculated from the LDA eigenvaluesφLDAn (r; R(t)) and the exchange-correlation potential µLDAxc (r) asεQPn (R(t)) ≈ εLDAn (R(t)) +〈φLDAn (R(t))∣∣ [Σxc(εQPn (R(t)))− µLDAxc ]∣∣ φLDAn (R(t))〉, (9)where Σxc(εQPn (R(t))) is the exchange-correlation part of the self-energy, which does notinclude the Hartree term. The self-energy Σxc(εQPn (R(t))) in Equation (9) explicitly dependson the resulting QP energy εQPn (R(t)). The QP energy obtained in a previous time stepcan be used as the argument for self-energy in the present time step. Therefore, therenormalization procedure using the Z factor [36] is not required in the present time-dependent approach. The usage of the QP energies εQPn (R(t)), which are obtained byEquation (9) in Equation (6), is the distinguishing feature of the TDGW approach. Exceptfor this difference, everything else remains the same as in the TDDFT dynamics formulation.The Newtonian equation of motion for NA-ES-TDGW-MD is approximated asMId2RI(t)dt2 = −∇RI E(N)ref (R(t), t)≈ −∇RI E(N)LDA(R(t), t), (10)where E(N)ref (R(t), t) and E(N)LDA(R(t), t) denote the GW total energy and the approximateLDA total energy, respectively, of an N-electron reference state. Although the exchange-correlation contribution to the total energy is approximated by its LDA form, the wavepack-ets used for TD charge densities are updated by Equation (6). Therefore, E(N)LDA(R(t), t) isnot the typical LDA total energy, but it includes the QP effects via the QP HamiltonianHQPmixed(r; R(t)).2.3. Ab Initio Cloning in NA-ES-TDGW-MDIn our approach, the nuclear trajectory evolves as the gradient of the average potentialgenerated by the electrons, which is a mean-field approximation in the Ehrenfest framework.In this mean-field approximation, the correlation (also known as quantum coherence)between the electron “motion” and nuclear trajectory is neglected. This suffers from theproblem of strong non-adiabatic mixing particularly when two Born–Oppenheimer (BO)surfaces cross or approach each other. A proper description of such correlations requires aNanomaterials 2024, 14, 1775 5 of 18distinct classical trajectory for each electronic state, which is provided by an SH strategy,such as that proposed by Tully [37,38] or Makhov et al. [35]. In this work, we adopt the SHstrategy proposed by Makhov et al. [35] named ab initio (multiple) cloning in our dynamicsformalism, although our approach does not include multiple trajectory basis functions.A quantity that is used as a measure of quantum (de)coherence, i.e., hopping from amixed surface (with index λ) to a pure BO surface (with index n) is called the “breakingforce” Fbrλ→n. The indices λ and n typically indicate the newly occupied/empty QP level(OCC1/EMP1; see Section 3.1) in the QP representation. The breaking force is defined asFbrλ→n = (1 − |cnλ(t)|2)∆Fnλ(t), (11)where cnλ(t) are the coefficients in the expansion of the wavepacket as defined inEquation (5) or (7) and ∆Fnλ(t) is the deviation of the force felt by the nuclei on a pure BOsurface (−∇RI E(N)n (R(t), t)) from the mean-field force in Equation (10) (∇RI E(N)λ (R(t), t)):∆Fnλ(t) = ∇RI E(N)λ (R(t), t)−∇RI E(N)n (R(t), t). (12)The condition for a hop (or “clone” in Ref. [35]) is that the breaking acceleration abrλ→nis greater than a pre-decided threshold δclone,abrλ→n = |M−1Fbrλ→n| > δclone. (13)In this work, we define δclone = 3 × 10−6 a.u. for exploring the surface hop time as in ourprevious study [24].2.4. Computational DetailsSince the eigenstates span the complete Hilbert space, we use the all-electron mixed-basis approach [24,29,33,39,40] implemented in our home-grown ab initio package namedTohoku mixed basis orbitals (TOMBO) [41], in which the one-electron orbitals are expressedby both plane waves (PWs) and atomic orbitals (AOs). We use a simple cubic unit cell of12 Å, where the Coulomb interaction is spherically cut to avoid interactions with adjacentunit cells. The 14147 PWs corresponding to the 17.3 Ry cutoff energy and minimal numberof numerical AOs have finite values only within each non-overlapping atomic sphere.AOs are smoothly truncated by subtracting a smooth quadratic function, which has thesame amplitude and derivative at the atomic sphere surface [41]. This quadratic functionsmoothly connecting to the tail of the true orbital can be successfully described by alinear combination of PWs. The cutoff energies corresponding to 69 Ry and 7.7 Ry are set,respectively, for the exchange and correlation parts of the self-energy. 190 levels are usedin the spectral decomposition [24] and in the calculation of the correlation part of the self-energy. The generalized plasmon pole model [36] is used to avoid frequency integration.We performed tests to ensure that all these parameters are sufficient for obtaining goodconvergence of results. However, we increased the cutoff energies to 44 Ry (for PWs), 111 Ry(for exchange), and 11 Ry (for correlation) as well as the number of levels to 3500 when aseparate one-shot GW calculation was performed.The precursor to performing a TDGW molecular dynamics simulation is to obtainconverged electronic states at a given electron configuration within the LDA. Once conver-gence is achieved, the dynamics simulation is performed by updating the wavepacketsstepwise in ∆t = 0.01 fs intervals, where the Hamiltonian is expected to change veryslightly. The wavepackets are the same as the QP (Kohn–Sham) eigenfunctions for the ESreference at t = 0. We perform 20 sub-loops of electronic state updates after every updateof the atomic positions, to ensure that the total energy is conserved. The GW calculation isperformed only during the final 20th sub-loop, and the QP energies and atomic positionsare updated subsequently.Nanomaterials 2024, 14, 1775 6 of 183. ResultsThe initial geometry of the CH4 +Ni system is shown in the leftmost panel in Figure 1.The Ni atom is placed at a C−Ni bond distance of 2.00 Å, which is the stable distance forcarbon chemisorption on a Ni(111) surface [42]. Two H atoms (H1 and H2) are away fromthe Ni atom while the other two (H3 and H4) face the Ni atom. We impose symmetrybreaking of the initial geometry of CH4 by slightly elongating one of the C−H bonds facingthe Ni atom (H4) by 0.05 Å. The CH4 + Ni system has a triplet ground state with 20 α-spin(↑-spin) and 18 β-spin (↓-spin) electrons. The charge density distribution of each QP levelfor the GS reference of this geometry is depicted in Figure 1. Here, and hereafter, thecharge density implies the absolute value of the QP wavefunction-squared irrespective ofits occupancy.From a single-shot GW calculation for the GS reference, the QP energies are obtained aspresented in Table 1. Since the 14th level, which corresponds to the HOMO of the methanefragment (indicated by ‘*’ in Table 1), is slightly shallower for β-spin than that for α-spin,we chose the β-spin for electron excitations to the LUMO. The QP energies of the 14thHOMO−4β (HOMO of CH4), 13th HOMO−5β (HOMO−1 of CH4), and LUMOβ levelsare, respectively, εHOMO−4β= −14.9 eV, εHOMO−5β= −15.1 eV, and εLUMOβ= −1.0 eV.Therefore, the ionization potential (IP) of the methane fragment at the GS of the CH4 + Nisystem, i.e., the energy required to remove one electron from the 14th HOMO−4β level, is−εHOMO−4β= 14.9 eV which is higher than that of pristine methane (GW: 13.7± 0.5 eV [24]and experiments: 12.6–14.8 eV [43–45]), while the electron affinity (EA) of the CH4 + Nisystem at the GS, i.e., the energy obtained by adding one electron to the 19th LUMOβ level,is −εLUMOβ= 1.0 eV. The higher IP and positive EA values in CH4 + Ni are attributed tothe presence of the Ni d orbitals.The slightly distorted CH4 still has a nearly two-fold degenerate HOMOβ (13th and14th β-spin levels of the CH4 + Ni system), whose QP wavefunctions exhibit two differentorientations with respect to the Ni atom, see Figure 1d,e. Therefore, the trajectory differsdepending on which level is excited. Here, we show two significantly different trajectorieswhen an electron is excited from the 13th/14th β-spin level to the 19th empty (LUMO)level. The 13th→19th level excitation involves the movement of the two H atoms facing Ni(H3 and H4 in Figure 1) while the 14th→19th level excitation involves the movement of thetwo H atoms away from Ni (H1 and H2 in Figure 1). We first present the latter case beforepresenting the former.HOMO−8 HOMO−7α-spinHOMO−6 HOMO−5 HOMO−4 HOMO−3 to HOMO LUMOβ-spinxyzxyzxyzxyzxyzxyzxyz(a) (b)GS structure(c) (d) (e) (f) (h)HOMO−3 HOMO(g)xyz2.00 ÅH4H2 H1H3xyzFigure 1. Initial geometry (leftmost panel) and (a–h) charge densities of the relevant α- and β-spinQP levels for the GS reference. The isosurface for the charge density plots was set at 1 × 10−8electrons/Å3. The directions of viewing are indicated for reference.Nanomaterials 2024, 14, 1775 7 of 18Table 1. The calculated QP energies (in units of eV) within the GWA for the GS reference using theinitial geometry. The 14th level (14*) corresponds to the CH4 HOMO.Level α-Spin (eV) β-Spin (eV)21 LUMO −1.0 LUMO+2 −0.320 HOMO −6.9 LUMO+1 −0.819 HOMO−1 −9.9 LUMO −1.018 HOMO−2 −10.3 HOMO −8.917 HOMO−3 −10.4 HOMO−1 −9.516 HOMO−4 −10.6 HOMO−2 −9.315 HOMO−5 −10.6 HOMO−3 −9.414* HOMO−6 −15.0 HOMO−4 −14.913 HOMO−7 −15.7 HOMO−5 −15.112 HOMO−8 −16.0 HOMO−6 −15.811 HOMO−9 −23.8 HOMO−7 −23.63.1. H Ejection Opposite to the Ni SideIn this case, we excite an electron from the slightly shallower HOMOβ level of CH4,which is the 14th HOMO−4β level of the CH4 + Ni system, to the empty Ni 3d level, whichis the 19th LUMOβ level. Now, the original HOMO−4β is called “EMP1β” and the originalLUMOβ, “OCC1β”. The labels such as OCC#α/β and EMP#α/β with # = 1, 2, ... indicate thelower occupied and higher empty levels, respectively. This nomenclature is based on thelevel ordering at t = 0 and will be constantly used even if the order of the levels changesduring the simulation.The time evolution of the QP energy eigenvalues εQPn (R(t)) of the α- and β-spin levelsis shown, respectively, in Figure 2a,b, while, the charge densities of the different levels ofinterest are shown at t = 3.8, 10.8, 25.8, and 32.6 fs in a tabular format below the QP energyplots in Figure 3.The levels of interest are (a) OCC9 (GS 12th HOMO−8 level), (b) OCC8 (GS 13thHOMO−7 level), (c) OCC7 (GS 14th HOMO−6 level) for α-spin, and (d) OCC7 (GS 12thHOMO−6 level), (e) OCC6 (GS 13th HOMO−5 level), (f) EMP1 (GS 14th HOMO−4 level),(g) OCC5 (GS 15th HOMO−3 level), (h) OCC2 (GS 18th HOMO level), (i) OCC1 (GS 19thLUMO level) for β-spin. The alphabets a-i followed by the simulation time in femtoseconds(fs) are used to represent the charge density panels.At t = 0, the QP energies of the β-spin OCC1 and EMP1 levels are, respectively,εOCC1β= −6.7 eV [thick orange dotted line in Figure 2b] and εEMP1β= −9.6 eV [thick violetdashed-dotted-dashed line in Figure 2b]. They are the Ni 3d orbital and the CH4 bonding or-bital as seen in Figure 1h,e, respectively. According to EQPT, εHOMO−4β= EGS − EHOMO−4βand εOCC1β= Ephoto − EHOMO−4β, where EGS, Ephoto, and EHOMO−4βare the total ener-gies of the GS, the photoexcited state, and the (N − 1)-electron state with one electronmissing at the the 14th HOMO−4β level, respectively. Note that εHOMO−4βis equal to− IP of the methane fragment in the CH4 + Ni system. Similarly, EA = EGS − Eanionand εEMP1β= Eanion − Ephoto, where Eanion is the total energy of the anionic state withone electron added to the LUMO level of the GS. From these relations, the photoabsorp-tion energy (PAE) for this excitation Ephoto − EGS can be obtained in two different ways,εOCC1β− εHOMO−4β= −6.7− (−14.9) eV = 8.2 eV and −EA − εEMP1β= −1.0− (−9.6) eV= 8.6 eV. The similarity in these two values clearly indicates the accuracy of our calculation.The resulting PAE (the averaged value is 8.4 eV) is about 1.8 eV lower than that for pristineCH4 (10.2 eV [5–7]).Nanomaterials 2024, 14, 1775 8 of 18−20−15−10−50OCC9OCC8OCC7OCC4OCC3OCC2OCC6OCC5OCC1EMP1EMP2α-spin quasiparticle energy (eV)(a)0−20−15−10−5010 20 30 40Time (fs)β-spin quasiparticle energy (eV)(b)OCC7OCC6OCC1OCC2OCC4OCC5OCC3EMP1EMP2EMP3Figure 2. Time evolution of (a) α-spin and (b) β-spin QP energy eigenvalues εQPn (R(t)).Nanomaterials 2024, 14, 1775 9 of 18OCC9 OCC8α-spinxyzxyz(a3.8) (b3.8)3.8 fsα/βt10.8 fs(a10.8) (b10.8)25.8 fs(a25.8) (b25.8)32.6 fs(a32.6) (b32.6)OCC7 OCC6 OCC1OCC5 to OCC2EMP1β-spinxyzxyzxyzxyzxyz(d3.8) (e3.8) (f3.8) (g3.8) (h3.8)(d10.8) (e10.8) (f10.8) (g10.8) (h10.8)OCC5 OCC2(i10.8)(d25.8) (e25.8) (f25.8) (g25.8) (h25.8) (i25.8)(d32.6) (e32.6) (f32.6) (g32.6) (h32.6) (i32.6)OCC5 OCC2OCC5 OCC2OCC5 OCC2(i3.8)xyzOCC7xyz(c3.8)(c10.8)(c25.8)(c32.6)xxyzz(f3.8)(f10.8)(f25.8)(f32.6)xxyz(a3.8)(a10.8)(a25.8)(a32.6)xxyzz(b3.8)(b10.8)(b25.8)(b32.6)xxyzz(c3.8)(c10.8)(c25.8)(c32.6)xxyz(d3.8)(d10.8)(d25.8)(d32.6)xxyzz(e3.8)(e10.8)(e25.8)(e32.6)xxyz(i10.8)(i25.8)(i32.6)(i3.8)(h3.8)(h10.8)OCC2(h25.8)(h32.6)OCC2OCC2OCC2xxyzxxyz(g3.8)(g10.8)OCC5(g25.8)(g32.6)OCC5OCC5OCC5Figure 3. Panels (a–i) represent the charge densities of the levels of interest for the H ejection oppositeto the Ni side for each time instant at t = 3.8, 10.8, 25.8, and 32.6 fs. The isosurface for the chargedensity plots was set at 1 × 10−8 electrons/Å3. The directions of viewing are indicated for referencein the uppermost panels (a).At t = 3.8 fs, we see that the charge densities of all QP levels [Figure 3a–i3.8] donot change significantly from those of the GS [Figure 1a–h], with the exception that theOCC8α and OCC7α levels for the ES reference [Figure 3b,c3.8] are the swapped HOMO−7α[Figure 1b] and HOMO−6α (not shown in Figure 1) levels for the GS reference. On theother hand for 4 ≤ t < 5 fs, the EMP1β level [thick violet dashed-dotted-dashed line inFigure 2b] crosses with the other OCC5β-OCC1β levels [OCC5 – brown dotted-dashed line,OCC4, light-green dotted-dashed line, OCC3, gray dotted-dashed-dotted line, OCC2, thinblue solid line, and OCC1, thick orange dotted line in Figure 2b]. The charge density ofEMP1β at t = 3.8 fs [Figure 3f3.8] is localized on the methane fragment but its population isshifted to the Ni side at t = 10.8 fs [Figure 3f10.8]. On the other hand, the charge density ofOCC1 at t = 3.8 and 10.8 fs [Figure 3i3.8,10.8] does not change and continues to exhibit a Ni3d character. Therefore, the electronic structure (mainly composed of the occupied orbitals)does not deviate significantly from the BO surface during level crossings, and the breakingacceleration does not exceed the SH threshold value δclone. Consequently, the simulation iscontinued without SH. As the simulation progresses, two H atoms opposite the Ni side (H1and H2 in Figure 1) begin to dissociate from the CH4 fragment of the CH4 + Ni system.Interesting parallels emerge in both the α- and β-spin QP energies in Figure 2a,b,respectively. While the QP energies of the OCC8α [red dashed line in Figure 2a] andthe EMP1β [thick violet dashed-dotted-dashed line in Figure 2b] levels increase (∼−17eV to ∼−10 eV for OCC8α and −9.6 eV to ∼−2 eV for EMP1β) with time, the QP ener-gies of the OCC7α [thick green dotted line in Figure 2a] and OCC6β [thick red dashedline in Figure 2b] levels oscillate analogously around ∼−15 eV. The resemblances in theNanomaterials 2024, 14, 1775 10 of 18temporal evolution of OCC9α and OCC7β, OCC8α and EMP1β, and OCC7α and OCC6βindicate that these orbitals are spin-paired with each other. This is further reflected in thesimilarities in the charge densities in Figure 3a3.8–32.6,d3.8–32.6 for OCC9α and OCC7β,Figure 3b3.8–32.6,f3.8–32.6for OCC8α and EMP1β, and Figure 3c3.8–32.6,e3.8–32.6 forOCC7α and OCC6β throughout the entire duration (t = 3.8, 10.8, 25.8, and 32.6 fs) ofthe simulation.The charge population in both OCC9α and OCC7β is initially spread toward the Niside (charge densities for t = 3.8, 10.8 fs) before being shifted towards the two ejected Hatoms (H1 and H2) at t = 25.8 fs [Figure 3a,d25.8]. Subsequently, the H1 and H2 atomsapproach each other and begin to form a H−H bond. Finally, at t = 32.6 fs, these α- andβ-spin levels become the completely isolated σ orbital of H−H [Figure 3a,d32.6]. At thistime, the H1−H2 bond length is 0.72 Å, which is very close to the experimental H2 bondlength, 0.74 Å [46], and the bond distances of C−H1 and C−H2 are 2.51 Å and 2.78 Å,respectively, which are both considerably large. We additionally perform a contour analysisof the charge densities shown in Figure 3a,d32.6 to obtain quantitative estimates of chargepopulations around the H2 fragment [Supplementary Information Figure S1a]. We comparethis with the charge contour of an isolated H2 molecule obtained from a single-point LDAcalculation [Supplementary Information Figure S1b]. We find an extremely good agreementbetween the TDGW-MD and the LDA results demonstrating the computational reliabilityof our predicted dynamics using TDGW-MD. Moreover, the QP energies of the OCC9αand OCC7β levels at t = 32.6 fs in Figure 2a,b are, respectively, −15.8 eV and −15.3 eV,which are also very close to the − IP of the hydrogen molecule 15.4 eV [47]. Therefore, weconclude that an isolated hydrogen molecule H2 is created as a product of this photolysisreaction in this very short time period. Based on this trajectory, we speculate that thephotolysis reaction would be CH4 + Ni h̄ω=8.4 eV−−−−−−→ CH2−Ni + H2.3.2. H Ejection Towards the Ni SideWe excite an electron from the slightly deeper HOMO−1β level of CH4, which is the13th HOMO−5β level of the CH4 + Ni system, to the 19th empty Ni 3d LUMOβ level. Theoriginal HOMO−5β is now empty and is called “EMP1β”, while the original LUMOβ isnow occupied and is called “OCC1β”. Figure 4a,b shows the early time behavior of the QPenergy eigenvalues εQPn (R(t)) of the α- and β-spin levels, respectively. The TDGW and theBO charge densities for the ES reference at t = 1.7 fs and t = 3.7 fs are shown in Figure 5.At t = 0, the QP energy of EMP1β [violet dashed-dotted-dashed line in Figure 4b]is εEMP1β= −9.8 eV and the QP energy of OCC1β [red dashed line in Figure 4b] isεOCC1β= −6.9 eV. They are the CH4 bonding orbital and the Ni 3d orbital as seen inthe charge density plots in Figure 5a5,7 at t = 1.7 fs. From EQPT, the required energyfor this photoexcitation can be calculated as either the difference between εOCC1βof thisphotoexcited state and εHOMO−5βof the GS, or the difference between εLUMOβof the GSstate, i.e., −EA, and εEMP1βof this photoexcited state. The former is −6.9 − (−15.1) eV= 8.2 eV and the latter is −1.0 − (−9.8) eV = 8.8 eV, which are close to each other. Thisagain demonstrates the accuracy of our simulation. The resulting PAE (the averaged valueis 8.5 eV) is about 1.7 eV lower than that for pristine CH4 (10.2 eV [5–7]).At t = 1.7 fs, the TDGW and BO charge densities are quite similar to each otheras seen in Figure 5a1–7,b1–7. They also resemble the GS charge densities in Figure 1a–h.This indicates that the non-adiabatic effect is very weak at this time. Indeed, the break-ing acceleration abrλ→n = 6.13 × 10−10 a.u. is much smaller than the threshold valueδclone = 3 × 10−6 a.u. In the 2.5 ≤ t ≤ 3.0 fs time interval, we see level crossings betweenthe EMP1β level [violet dashed-dotted-dashed line in Figure 4b] and the other OCC5β-OCC1β levels [OCC5, brown dotted-dashed line, OCC4, orange dotted-dashed line, OCC3,gray dotted-dashed-dotted line, OCC2, blue dashed line, and OCC1, red dashed line inFigure 4b]. In contrast, there is no such level crossing in the lower α-spin levels besides thatbetween the EMP1 [pink dashed line in Figure 4a] and EMP2 [light blue dotted-dashed linein Figure 4a] levels at t = 2.6 fs.Nanomaterials 2024, 14, 1775 11 of 18−20−15−10−50OCC9OCC8OCC7OCC5OCC3OCC2OCC6OCC4OCC1EMP1EMP2α-spin quasiparticle energy (eV)(a)OCC7OCC1OCC6OCC2OCC4OCC5OCC3EMP1EMP2EMP34 5 6 7−20−15−10−500 1 2 3Time (fs)β-spin quasiparticle energy (eV)(b)Figure 4. Early time behavior of (a) α-spin, (b) β-spin QP energy eigenvalues εQPn (R(t)) for the Hejection towards the Ni side.Nanomaterials 2024, 14, 1775 12 of 18xyzxyzxyzxyzxyzxyzxyz1.79 × 10−5a. u.6.13 × 10−10a. u.α-spin β-spinOCC9 OCC8 OCC7 OCC7 EMP1 OCC6 OCC1TDGWBOBO1.7 fs3.7 fstα/βQP level(a1) (a2) (a3) (a4) (a5) (a6) (a7)(b1) (b2) (b3) (b4) (b5) (b6) (b7)(c1) (c2) (c3) (c4) (c5) (c6) (c7)(d1) (d2) (d3) (d4) (d5) (d6) (d7)TDGWFigure 5. Comparison of TDGW charge densities (a,c) with BO charge densities (b,d) for the ESreference at t = 1.7 fs (a1–7) and (b1–7) and t = 3.7 f s (c1–7) and (d1–7). The isosurface for the chargedensity plots was set at 1 × 10−8 electrons/Å3. The directions of viewing are indicated for referencein the uppermost panels (a).At t = 3.7 fs, the breaking acceleration abrλ→n = 1.79 × 10−5 a.u. exceeds the thresholdvalue δclone = 3 × 10−6 a.u., similar to the case of pristine methane [24]. Indeed, the TDGWcharge densities at this time instant [Figure 5c1–7] are somewhat different from the BOcharge densities with the same geometry [Figure 5d1–7]. For example, those of the OCC9αand OCC8α levels are different as if the left and right charge lobes were reversed. Moreover,those of the OCC1β level are slightly different in terms of the Ni 3d character. Therefore,the simulation undergoes SH at t = 3.7 fs. In our Ehrenfest framework, we perform SH tothe BO surface for the GS reference.The time evolution of the QP energies (for the simulation with SH at t = 3.7 fs) areshown in Figure 6a,b, respectively, for α- and β-spins while the trajectory of this simulationis presented in Figure 7.At the SH time, EMP1β level [violet dashed-dotted-dashed line in Figure 6b] jumps upabove the vacuum level, while OCC1β level [red dashed line in Figure 6b] falls down to∼−14 eV. Around this SH time, the two H atoms in the Ni side (H3 and H4 in Figure 1)start to move away from the C atom and approach the Ni atom. Thereafter, OCC9α-OCC7α, OCC7β, OCC6β and OCC1β levels oscillate in their corresponding QP energiesaround −14 eV [Figure 6a,b]. The similarity in their temporal variation is reflected in thecharge densities shown in Figure 7a–e,g, implying that these orbitals are spin-paired witheach other.Nanomaterials 2024, 14, 1775 13 of 18−20−15−10−50OCC9OCC8OCC7OCC5OCC3OCC2OCC6OCC4OCC1EMP1EMP2α-spin quasiparticle energy (eV)(a)OCC7OCC1OCC6OCC2OCC4OCC5OCC3EMP1EMP2EMP340 50 60 70−20−15−10−500 10 20 30Time (fs)β-spin quasiparticle energy (eV)(b)Figure 6. Time evolution of (a) α-spin, (b) β-spin QP energy eigenvalues εQPn (R(t)).With time, the H3 and H4 atoms begin to move away from the Ni atom. At t = 22.2 fs,H3 and H4 return to the methane side, but around the [40.6, 57.6] fs time interval, H3 is boundto the Ni atom. However, this bonding is not strong with the H3−Ni bond being broken, anda subsequent return of H3 toward methane. Finally, at t = 72.6 fs, the original geometry isNanomaterials 2024, 14, 1775 14 of 18virtually retrieved. Indeed, the charge densities at t = 72.6 shown in Figure 7a–g are veryclose to those of the GS charge densities in Figure 1. At the same time, interestingly, the Niatom starts to move away from the CH4 fragment which is clearly evident at t = 72.6 fsfrom the C−Ni distance around 2.29 Å. From Figure 6a,b around t = 72.6 fs, we see that theQP energies εOCC9α-εOCC7α, εOCC7β, εOCC6β, and εOCC1βoscillate in their QP energies aroundthe mean value of −14 eV, which is very close to − IP = −13.7 ± 0.5 eV of the pristinemethane [24]. (Experimental IP is 12.6–14.8 eV [43–45].) In addition, all empty levels (EMP1,EMP2, etc.) are above the vacuum level or just close to the vacuum level regardless of the spin,as seen in Figure 6a,b, reflecting the negative EA of the pristine methane molecule. All theseobservations indicate that the combined CH4−Ni system is undergoing dissociation into theCH4 and Ni fragments, i.e, CH4 +Ni −−→ CH4−Ni h̄ω=8.5 eV−−−−−−→ CH4 +Ni.OCC9 OCC8α-spinxyzxyz(a)72.6 fsα/βt OCC7 OCC6 OCC1EMP1β-spinxyzxyzxyzxyz(f)OCC7xyz0.0 fs3.7 fs(SH)22.2 fs 40.6 fs 57.6 fs 72.6 fsxyzxxyzxxyzxxyzxxyzxxyzxxxyz(a)xxyz(f)(b) (d) (e) (g)(c)Figure 7. Time evolution of atomic geometry for the H ejection in the Ni side. Panels (a–g) representcharge densities of relevant QP levels at t = 72.6 fs. The isosurface for the charge density plots was set at1× 10−8 electrons/Å3. The directions of viewing are indicated for reference.4. DiscussionWe have considered two possible pathways of the photolysis of methane in the pres-ence of a Ni atom depending on either the β-spin 14th (HOMO−4β) level (HOMO of themethane fragment) or the 13th (HOMO−5β) level (HOMO−1 of the methane fragment)being excited to the 19th (LUMOβ) level. In both the cases, we estimated the required PAEin two different ways under EQPT. One approach is εOCC1βof the ES minus εHOMO−4β(orεHOMO−5β) of the GS and the other is −EA minus εEMP1βof the ES. (IP of the methanefragment is −εHOMO−4β.) The resulting PAEs via either approach are close to each other,suggesting the accuracy of the present calculation. The PAEs are 8.4 eV and 8.5 eV, respec-tively, for the 14th and the 13th level excitations which are lower than that of 10.2 eV forpristine methane, indicating the reduction in the PAE due to the existence of a Ni atom.This reduction in PAE is due to the excitation of an electron from a C 2p orbital [Figure 1d,e]to an intermediate 3d orbital of Ni [Figure 1h] instead of to the higher energy 3s orbitalof C (in the pristine CH4 case). As a matter of fact, photolysis of methane on a Pt(111)surface has experimentally been shown to occur at an excitation energy of around 6.4 eV(∼193 nm) [48], which is significantly lower than 10.2 eV (∼122 nm) in the absence of such asubstrate further exemplifying the crucial role played by transition metal atoms in makingphotochemical reactions more accessible.It is noteworthy that, in both pathways, two hydrogen atoms are simultaneouslyejected from the methane molecule at the beginning, which is a distinct characteristic of thepresent results. On the contrary, in the pristine methane case, only one H atom is ejectedby a single photon absorption. Why are two H atoms ejected simultaneously from CH4despite undergoing an excitation via a single photon absorption in the presence of a NiNanomaterials 2024, 14, 1775 15 of 18atom? The reason for this is that the C−H bonds in the CH4 molecule become considerablyweakened while the C−Ni bonding becomes stronger via p-d orbital hybridization.When an electron is excited from the 14th level to the 19th level, two H atoms oppositeNi are ejected from the methane molecule, and they eventually combine to become an isolatedhydrogen molecule. The QP energy of the HOMO level of this ejected hydrogen moleculeis ∼−15.5 eV irrespective of the spin, which is close to the − IP of the isolated hydrogenmolecule 15.4 eV [47]. The QP wavefunction corresponding to the ejected hydrogen moleculeis completely localized at the hydrogen molecule, and there is no charge fraction remaining inthe CH2−Ni side (see also Supplementary Information (SI) for the contour plots of the chargedensity). Thus, we can conclude that a hydrogen molecule is produced as a product. This is asensational result because two hydrogen atoms are ejected from one methane molecule by asingle photon absorption and they automatically combine to become a hydrogen molecule.This clearly demonstrates the important role played by Ni in providing a more efficientroute for photolysis. This pathway leads to the formation of the H2 molecule as there doesnot exist an ‘absorber’ such as Ni to impede the motion of the two H atoms. This reactionis closely related to the thermocatalytic or solar-aided decomposition of methane, where atransition metal atom or cluster mediates the growth of nanocarbon materials together withthe production of hydrogen molecules. For example, there is a possibility that, if two methanemolecules are attached to a Ni cluster, the remaining C atoms can combine to produce theC2 dimer subsequent to several photolysis reactions that may result in the dissociation offour hydrogen molecules. With the addition of more methane molecules, this process may becontinued to produce carbon nanomaterials.In contrast, when two hydrogen atoms facing the Ni atom are ejected, the ejected H atomscannot escape because of the strong attractive forces exerted by both the Ni and methylene(CH2) fragments. The two H atoms initially approach the Ni atom due to the inertia of motionbut eventually return to the CH4 fragment. There are three competing factors that govern theentire dynamic: the inertia of motion of the ejected H, the Coulombic interaction between theH and Ni atoms, and the covalent interaction between the H and C atoms. During the initialstage of the photolysis, as the ejected H atoms approach the Ni atom, their inertia of motioncarries them away from the C atom; see Figure 7. After SH to the GS-BO surface, the strengthof the covalent interaction exceeds this inertia and forces the ejected H atoms to return. Withthis return, the inertia of motion in the reverse direction increases leading to the H atomsovershooting their original positions at t = 22.2 fs. In this process, the overall CH4 geometrystrongly deviates from its stable configuration, which initiates the rebound of one of the Hatoms toward Ni. Mediated by a weak Coulomb interaction as a result of charge transferfrom this H atom to the Ni atom, this H temporarily bonds with Ni around t = 40.6–57.6 fs.However, the influence of the C atom continues to persist, leading to the return of this ejectedH to the original position at t = 72.6 fs [Figure 7a–g]. Because of these three competing factors,none of the H atoms can get dissociated after ejection, but rather prefer to return to (nearly)recover the original methane geometry. At the same time, during the reversal of the H atomfrom Ni towards CH4, the increase in the inertia of H in the opposite direction initiates theconcerted motion of CH4 away from Ni. This can be seen not only from the atomic trajectorybut also from the fact that methane-derived QP levels have energy values similar to −IP ofpristine methane. Therefore, in this case, we conclude that the combined CH4−Ni system isdissociated into CH4 and Ni.We next comment on the reliability of our simulation results despite the short MDsimulation times. Photolysis reactions are generally ultrafast [49] as they usually completewithin several tens of femtoseconds. Therefore, even though the simulation times seem tobe short in our work, the match between QP energies and experimental values for both thetrajectories is clearly evident. Therefore, we argue that the present results are reliable.In addition to the reliability and accuracy of our results, these were obtained inreasonable wall clock times. The simulation of H ejection away from/toward Ni took6/10 days using four nodes with 48 MPI processes on the Numerical Materials Simulatorsupercomputer at the National Institute for Materials Science.Nanomaterials 2024, 14, 1775 16 of 18Through our simulations, we show how two H atoms can be dissociated via a singlephoton absorption eventually leading to the formation of an isolated H2 molecule. Thisis achieved at a much lower PAE of 8.4 eV, compared to pristine methane (10.2 eV). Over-all, our work presents an opportunity for experimental verification of our observations.The results obtained from our simulations are accurate as they are obtained without anyadjustable/fitting parameters and are based on sound mathematical principles.Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/nano14221775/s1, Figure S1: Contour plots of the H2 charge density; MDTrajctory S1: away_from_Ni_MD.xyz; MD Trajectory S2: towards_Ni_MD.xyz; see Data AvailabilityStatement to view these .xyz files.Author Contributions: Conceptualization, K.O., R.S. and Y.K.; methodology, A.M. and K.O.; software,A.M., K.O., R.S. and Y.K.; validation, A.M., K.O., R.S. and Y.K.; investigation, A.M.; resources, R.S.;writing—original draft preparation, A.M. and K.O.; writing—review and editing, R.S. and Y.K.;visualization, A.M.; supervision, R.S.; project administration, R.S.; funding acquisition, R.S., Y.K. andK.O. All authors have read and agreed to the published version of the manuscript.Funding: This research was funded by the Japan Society for the Promotion of Science (JSPS) KAK-ENHI grant numbers 21H01607, 21H01877, 23K21094 and 24K01149, and by Asian Office of AerospaceResearch and Development (AOARD) grant number FA2386-22-1-4024.Data Availability Statement: The data that support the findings of this study are available within thearticle and its supplementary files, away_from_Ni_MD.xyz and towards_Ni_MD.xyz. The xyz filescan be visualized using the VMD software package [50]. The TOMBO executable used for performingNA-ES-TDGW-MD simulations and one-shot GW calculations is available upon request.Acknowledgments: The calculations in this study were performed using the Numerical MaterialsSimulator at the National Institute for Materials Science (NIMS); System B at the Institute for SolidState Physics (ISSP), the University of Tokyo, the MASAMUNE supercomputing system at the Centerfor Computational Materials Science, Institute for Materials Research (IMR), Tohoku University(Proposal Nos. 202308-SCKGE-0216, 202312-SCKXX-0206 and 202311-SCKXX-0501) and the JHPCNsystem project (Project IDs: jh220037 and jh230038). Y.K. was supported by (i) the National Science,Research, and Innovation Fund (NSRF) (NRIIS Project No. 90465); (ii) Thailand Science Research andInnovation (TSRI) and (iii) Suranaree University of Technology (SUT).Conflicts of Interest: The authors declare no conflicts of interest.AbbreviationsThe following abbreviations are used in this manuscript:DFT density functional theoryTD time-dependentTDDFT time-dependent density functional theoryLDA local density approximationALDA adiabatic local density approximationMD molecular dynamicsES excited stateGS ground stateQP quasiparticleEQPT extended quasiparticle theoryvac vacuum levelGWA GW approximationSH surface hoppingTDGW Non-adiabatic excited-state time-dependent GWPW plane waveAO atomic orbitalHOMO highest occupied molecular orbitalLUMO lowest unoccupied molecular orbitalIP ionization potentialhttps://www.mdpi.com/article/10.3390/nano14221775/s1https://www.mdpi.com/article/10.3390/nano14221775/s1Nanomaterials 2024, 14, 1775 17 of 18EA electron affinityPAE photoabsorption energyReferences1. Mokhatab, S.; Mak, J.Y.; Valappil, J.V.; Wood, D.A. Handbook of Liquefied Natural Gas, 3rd ed.; Gulf Professional Publishing: Boston,MA, USA, 2014; pp. 1–624.2. Gunaseelan, V.N. Anaerobic digestion of biomass for methane production: A review. Biomass Bioenergy 1997, 13, 83–114. [CrossRef]3. Züttel, A. Hydrogen storage methods. Naturwissenschaften 2004, 91, 157–172. [CrossRef] [PubMed]4. Yousefi, M.; Donne, S. Experimental study for thermal methane cracking reaction to generate very pure hydrogen in small ormedium scales by using regenerative reactor. Front. Energy Res. 2022, 10, 1. [CrossRef]5. Laufer, A.H.; McNesby, J.R. Photolysis of methane at 1236-Å: Quantum yield of hydrogen formation. J. Chem. Phys. 1968, 49,2272–2278. [CrossRef]6. Rebbert, R.; Ausloos, P. Photolysis of methane: Quantum yield of C(1D) and CH. J. Photochem. 1972, 1, 171–176. [CrossRef]7. Brownsword, R.; Hillenkamp, M.; Laurent, T.; Vatsa, R.; Volpp, H.-R.; Wolfrum, J. Quantum yield for H atom formation in themethane dissociation after photoexcitation at the lyman-α (121.6 nm) wavelength. Chem. Phys. Lett. 1997, 266, 259–266. [CrossRef]8. Lee, L.C.; Chiang, C.C. Fluorescence yield from photodissociation of CH4 at 1060-1420 Å. J. Chem. Phys. 1983, 78, 688–691.[CrossRef]9. Mebel, A.M.; Lin, S.-H.; Chang, C.-H. Theoretical study of vibronic spectra and photodissociation pathways of methane. J. Chem.Phys. 1997, 106, 2612–2620. [CrossRef]10. York, A.P.; Xiao, T.; Green, M.L. Brief overview of the partial oxidation of methane to synthesis gas. Top. Catal. 2003, 22, 345–358.[CrossRef]11. Matsumura, Y.; Nakamori, T. Steam reforming of methane over nickel catalysts at low reaction temperature. Appl. Catal. A Gen.2004, 258, 107–114. [CrossRef]12. Ashik, U.P.M.; Daud, W.M.A.W.; Abbas, H.F. Production of greenhouse gas free hydrogen by thermocatalytic decomposition ofmethane–A review. Renew. Sust. Energ. Rev. 2015, 44, 221–256. [CrossRef]13. Keipi, T.; Tolvanen, H.; Konttinen, J. Economic analysis of hydrogen production by methane thermal decomposition: Comparisonto competing technologies. Energy Convers. Manag. 2018, 159, 264–273. [CrossRef]14. Qian, J.X.; Chen, T.W.; Enakonda, L.R.; Liu, D.B.; Basset, J.-M.; Zhou, L. Methane decomposition to pure hydrogen and carbonnano materials: State-of-the-art and future perspectives. Int. J. Hydrogen Energy 2020, 45, 15721–15743. [CrossRef]15. Fan, Z.; Weng, W.; Zhou, J.; Gu, D. Xiao, W. Catalytic decomposition of methane to produce hydrogen: A review. J. Energy Chem.2021, 58, 415–430. [CrossRef]16. Pham, C.Q.; Siang, T.J.; Kumar, P.S.; Ahmad, Z.; Xiao, L.; Bahari, M.B.; Cao, A.N.T.; Rajamohan, N.; Qazaq, A.S.; Kumar, A.; et al.Production of hydrogen and value-added carbon materials by catalytic methane decomposition: A review. Environ. Chem. Lett.2022, 20, 2339–2359. [CrossRef]17. Abanades, S.; Kimura, H.; Otsuka, H. Hydrogen production from CO2-free thermal decomposition of methane: Design andon-sun testing of a tube-type solar thermochemical reactor. Fuel Process. Tech. 2014, 122, 153–162. [CrossRef]18. Pinilla, J.L.; Torres, D.; Lázaro, M.J.; Suelves, I.; Moliner, R.; Cañadas, I.; Rodríguez, J.; Vidal, A.; Martínez, D. Metallic andcarbonaceous-based catalysts performance in the solar catalytic decomposition of methane for hydrogen and carbon production.Int. J. Hydrogen Energy 2012, 37, 9645–9655. [CrossRef]19. Boretti, A. A perspective on the production of hydrogen from solar-driven thermal decomposition of methane. Int. J. HydrogenEnergy 2021, 46, 34509–34514. [CrossRef]20. Jia, S.; Cao, X.; Yuan, X.; Yu, K.-T. Multi-objective topology optimization for the solar thermal decomposition of methane reactorenhancement. Chem. Eng. Sci. 2021, 231, 116265. [CrossRef]21. Hohenberg, P.; Kohn, W. Inhomogeneous electron gas. Phys. Rev. 1964, 136, B864–B871. [CrossRef]22. Runge, E.; Gross, E.K.U. Density-functional theory for time-dependent systems. Phys. Rev. Lett. 1984, 52, 997–1000. [CrossRef]23. Petersilka, M.; Gossmann, U.J.; Gross, E.K.U. Excitation energies from time-dependent density-functional theory. Phys. Rev. Lett.1996, 76, 1212–1215. [CrossRef] [PubMed]24. Manjanath, A.; Sahara, R.; Ohno, K.; Kawazoe, Y. Non-adiabatic excited-state time-dependent GW molecular dynamics (TDGW)satisfying extended Koopmans’ theorem: An accurate description of methane photolysis. J. Chem. Phys. 2024, 160, 184102.[CrossRef] [PubMed]25. Ohno, K.; Ono, S.; Isobe, T. A simple derivation of the exact quasiparticle theory and its extension to arbitrary initial excitedeigenstates. J. Chem. Phys. 2017, 146, 084108. [CrossRef] [PubMed]26. Ohno, K.; Esfarjani, K.; Kawazoe, Y. Computational Materials Science: From Ab Initio to Monte Carlo Methods, 2nd ed.; Springer:Heidelberg, Germany, 2018; pp. 1–427.27. Morrell, M.M.; Parr, R.G.; Levy, M. Calculation of ionization potentials from density matrices and natural functions, and thelong-range behavior of natural orbitals and electron density. J. Chem. Phys. 1975, 62, 549–554. [CrossRef]28. Smith, D.W.; Day, O.W. Extension of Koopmans’ theorem. I. Derivation. J. Chem. Phys. 1975, 62, 113–114. [CrossRef]29. Sahara, R.; Mizuseki, H.; Sluiter, M.H.F.; Ohno, K.; Kawazoe, Y. Effect of a nickel dimer on the dissociation dynamics of ahydrogen molecule. RSC Adv. 2013, 3, 12307–12312. [CrossRef]http://doi.org/10.1016/S0961-9534(97)00020-2http://dx.doi.org/10.1007/s00114-004-0516-xhttp://www.ncbi.nlm.nih.gov/pubmed/15085273http://dx.doi.org/10.3389/fenrg.2022.971383http://dx.doi.org/10.1063/1.1670396http://dx.doi.org/10.1016/0047-2670(72)85005-6http://dx.doi.org/10.1016/S0009-2614(96)01526-6http://dx.doi.org/10.1063/1.444812http://dx.doi.org/10.1063/1.473410http://dx.doi.org/10.1023/A:1023552709642http://dx.doi.org/10.1016/j.apcata.2003.08.009http://dx.doi.org/10.1016/j.rser.2014.12.025http://dx.doi.org/10.1016/j.enconman.2017.12.063http://dx.doi.org/10.1016/j.ijhydene.2020.04.100http://dx.doi.org/10.1016/j.jechem.2020.10.049http://dx.doi.org/10.1007/s10311-022-01449-2http://dx.doi.org/10.1016/j.fuproc.2014.02.002http://dx.doi.org/10.1016/j.ijhydene.2012.03.075http://dx.doi.org/10.1016/j.ijhydene.2021.07.234http://dx.doi.org/10.1016/j.ces.2020.116265http://dx.doi.org/10.1103/PhysRev.136.B864http://dx.doi.org/10.1103/PhysRevLett.52.997http://dx.doi.org/10.1103/PhysRevLett.76.1212http://www.ncbi.nlm.nih.gov/pubmed/10061664http://dx.doi.org/10.1063/5.0202590http://www.ncbi.nlm.nih.gov/pubmed/38716844http://dx.doi.org/10.1063/1.4976553http://www.ncbi.nlm.nih.gov/pubmed/28249434http://dx.doi.org/10.1063/1.430509http://dx.doi.org/10.1063/1.430253http://dx.doi.org/10.1039/c3ra40928gNanomaterials 2024, 14, 1775 18 of 1830. Burghgraef, H.; Jansen, A.P.J.; van Santen, R.A. Electronic structure calculations and dynamics of methane activation on nickeland cobalt. J. Chern. Phys. 1994, 101, 11012–11020. [CrossRef]31. Sun, Q.; Li, Z.; Du, A.; Chen, J.; Zhu, Z.; Smith, S.C. Theoretical study of two states reactivity of methane activation on iron atomand iron dimer. Fuel 2012, 96, 291–297. [CrossRef]32. Wang, J.-T.; Chen, C.; Ohno, K.; Wang, E.; Chen, X.-L.; Wang, D.-S.; Mizuseki, H.; Kawazoe, Y. Atomistic nucleation and growthmechanism for single-wall carbon nanotubes. on catalytic nanoparticle. Nanotechnology 2010, 21, 115602. [CrossRef]33. Pham, T.N.; Ono, S.; Ohno, K. Ab initio molecular dynamics simulation study of successive hydrogenation reactions of carbonmonoxide producing methanol. J. Chem. Phys. 2016, 144, 144309. [CrossRef] [PubMed]34. Saalmann, U.; Schmidt, R. Non-adiabatic quantum molecular dynamics: Basic formalism and case study. Z. Phys. D Atom. Mol.Cl. 1996, 38, 153–163. [CrossRef]35. Makhov, D.V.; Glover, W.J.; Martinez, T.J.; Shalashilin, D.V. Ab initio multiple cloning algorithm for quantum nonadiabaticmolecular dynamics. J. Chem. Phys. 2014, 141, 054110. [CrossRef] [PubMed]36. Hybertsen, M.S.; Louie, S.G. Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies. Phys.Rev. B 1986, 34, 5390–5413. [CrossRef] [PubMed]37. Tully, J.C.; Preston, R.K. Trajectory surface hopping approach to non-adiabatic molecular collisions: The reaction of H+ with D2.J. Chem. Phys. 1971, 55, 562–572. [CrossRef]38. Tully, J.C. Molecular dynamics with electronic transitions. J. Chem. Phys. 1990, 93, 1061–1071. [CrossRef]39. Ohno, K.; Manjanath, A.; Kawazoe, Y.; Hatakeyama, R.; Misaizu, F.; Kwon, E.; Fukumura, H.; Ogasawara, H.; Yamada, Y.;Zhang, C.; et al. Extensive First-Principles Molecular Dynamics Study on the Li Encapsulation into C60 and Its ExperimentalConfirmation. Nanoscale 2018, 10, 1825–1836. [CrossRef]40. Ohtsuki, T.; Manjanath, A.; Ohno, K.; Inagaki, M.; Sekimoto, S.; Kawazoe, Y. Creation of Mo/Tc@C60 and Au@C60 and molecular-dynamics simulations. RSC Adv. 2021, 11, 19666–19672. [CrossRef]41. Ono, S.; Noguchi, Y.; Sahara, R.; Kawazoe, Y.; Ohno, K. TOMBO: All-electron mixed-basis approach to condensed matter physics.Comput. Phys. Commun. 2015, 189, 20–30. [CrossRef]42. Klinke, D.J.; Wilke, S.; Broadbelt, L.J. A Theoretical Study of Carbon Chemisorption on Ni(111) and Co(0001) Surfaces. J. Catal.1998, 158, 540–554. [CrossRef]43. Watanabe, K. Ionization potentials of some molecules. J. Chem. Phys. 1957, 26, 542–547. [CrossRef]44. Collin, J.E.; Delwiche, J. Ionization of methane and its electronic energy levels. Can. J. Chem. 1967, 45, 1875–1882. [CrossRef]45. Cook, P.A.; Ashfold, M.N.R.; Jee, Y.-J.; Jung, K.-H.; Harich, S.; Yang, X. Vacuum ultraviolet photochemistry of methane, silane andgermane. Phys. Chem. Chem. Phys. 2001, 3, 1848–1860. [CrossRef]46. Huber, K.P.; Herzberg, G. Molecular spectra and molecular structure. In Constants of Diatomic Molecules; Van Nostrand ReinholdCo.: New York, NY, USA, 1979; Volume IV.47. Herzberg, G.; Jungen, C. Rydberg series and ionization potential of the H2 molecule. J. Mol. Spectrosc. 1972, 41, 425–486. [CrossRef]48. Gruzdkov, Y.A.; Watanabe, K.; Sawabe, K.; Matsumoto, Y. Photochemical C−H bond activation of methane on a Pt(111) surface.Chem. Phys. Lett. 1994, 227, 243–247. [CrossRef]49. Troß, J.; Carter-Fenk, K.; Cole-Filipiak, N.C.; Schrader, P.; Word, M.; McCaslin, L.M.; Head-Gordon, M.; Ramasesha, K. Excited-state dynamics during primary C−I homolysis in acetyl iodide revealed by ultrafast core-level spectroscopy. J. Phys. Chem. A2023, 127, 4103-4114. [CrossRef]50. Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. J. Mol. Graph. 1996, 14, 33–38. [CrossRef]Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individualauthor(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury topeople or property resulting from any ideas, methods, instructions or products referred to in the content.http://dx.doi.org/10.1063/1.467852http://dx.doi.org/10.1016/j.fuel.2011.12.060http://dx.doi.org/10.1088/0957-4484/21/11/115602http://dx.doi.org/10.1063/1.4945628http://www.ncbi.nlm.nih.gov/pubmed/27083723http://dx.doi.org/10.1007/s004600050077http://dx.doi.org/10.1063/1.4891530http://www.ncbi.nlm.nih.gov/pubmed/25106573http://dx.doi.org/10.1103/PhysRevB.34.5390http://www.ncbi.nlm.nih.gov/pubmed/9940372http://dx.doi.org/10.1063/1.1675788http://dx.doi.org/10.1063/1.459170http://dx.doi.org/10.1039/C7NR07237Fhttp://dx.doi.org/10.1039/D0RA10196Fhttp://dx.doi.org/10.1016/j.cpc.2014.11.012http://dx.doi.org/10.1006/jcat.1998.2175http://dx.doi.org/10.1063/1.1743340http://dx.doi.org/10.1139/v67-299http://dx.doi.org/10.1039/b100248lhttp://dx.doi.org/10.1016/0022-2852(72)90064-1http://dx.doi.org/10.1016/0009-2614(94)00849-3http://dx.doi.org/10.1021/acs.jpca.3c01414http://dx.doi.org/10.1016/0263-7855(96)00018-5 Introduction Method Time-Dependent QP Equation One-Shot GW Within TDQP Ab Initio Cloning in NA-ES-TDGW-MD Computational Details Results H Ejection Opposite to the Ni Side H Ejection Towards the Ni Side Discussion References