Modeling of High-Harmonic Generation in the C60 Fullerene Using Ab Initio, DFT-Based, and Semiempirical Methods

We report calculations of the high-harmonic generation spectra of the C60 fullerene molecule carried out by employing a diverse set of real-time time-dependent quantum chemical methods. All methodologies involve expanding the propagated electronic wave function in bases consisting of the ground and singly excited time-independent eigenstates obtained through the solution of the corresponding linear-response equations. We identify the correlation and exchange effect in the spectra by comparing the results from methods relying on the Hartree–Fock reference determinant with those obtained using approaches based on the density functional theory with different exchange–correlation functionals. The effect of the full random-phase approximation treatment of the excited electronic states is also analyzed and compared with the configuration interaction singles and the Tamm–Dancoff approximation. We also showcase the fact that the real-time extension of the semiempirical method INDO/S can be effectively applied for an approximate description of laser-driven dynamics in large systems.


Introduction
When atoms and molecules are subjected to extremely intense laser fields, their electron densities undergo rapid oscillations, emitting electromagnetic radiation that contains up to multiple hundreds of harmonic frequencies of the incident light.The discovery of this phenomenon, known as high-harmonic generation (HHG), has ushered in the era of attophysics. 1,2It has enabled the routine production of ultrashort coherent electromagnetic pulses, making it possible to explore electron dynamics on previously unattainable timescales.
[14][15] HHG was initially observed in noble [16][17][18][19][20][21][22][23][24][25][26] and molecular [27][28][29] gases, with the former re-maining the most commonly used sources of harmonic radiation.In simple mono-or fewatomic systems a single HHG event can be explained with the well-known three-step model (3SM), [30][31][32][33][34] according to which an electron (1) is detached from the atom or molecule via tunneling ionization, (2) is accelerated away by the driving field and then reaccelerated back toward the parent ion when the field switches its sign and (3) recombines with the parent ion which leads to the emission of radiation burst.High-harmonic generation from gas targets, although relatively easy to achieve, is however hindered by a low conversion efficiency attributed to the low density of the medium. 35Thus, there is an ongoing search for novel, more efficient HHG sources.In recent years, thanks to the rapid development of experimental techniques, HHG was demonstrated to occur also in bulk solids, [36][37][38][39][40] liquids, 41 and nanostructures, [42][43][44][45][46][47] with the harmonic yield greatly exceeding that of atomic and molecular gases.The mechanism of harmonic generation in bulk media, although not fully revealed, is suspected to differ substantially from the 3SM. 36,38ogress in experimental discoveries in attophysics also necessitates the development of new theoretical methods capable of describing ultrafast electron dynamics in increasingly larger systems.[50][51][52] These approaches are characterized by moderate computational costs typical of their timeindependent counterparts, as well as a reasonable level of accuracy, allowing, for example, the consideration of multi-electron effects.4][55][56][57][58] Thanks to its highly favorable scaling with the number of electrons, it can be routinely applied not only to atoms [59][60][61][62][63][64][65][66] and simple molecules, [67][68][69][70][71] but also to complex organic 72 and biological [73][74][75] compounds, often yielding results qualitatively or even quantitatively consistent with the experimental data. 72,73,76This raises the question of whether it can perform equally well for even larger systems, such as nanostructured materials.
In this work we report calculations of the HHG spectra of arguably the most well-known carbon nanostructure, the C 60 fullerene molecule, using quantum chemical approaches coupled to Gaussian basis sets.Fullerenes are currently of great interest of attophysics, as experimental studies report their exceptionally high HHG yield, [77][78][79][80][81] significantly surpassing that of bulk carbon. 77,79,81This property is attributed to their high polarizability, as well as the occurence of the plasmon resonance at the fullerene surface. 82From the theoretical point of view, HHG in fullerenes has been studied either using an extention of the three-step model, 83,84 or by real-time simulations employing a variety of different approximate Hamiltonians.These include SAE-based models, 80,85 tight-binding models, [86][87][88] jelliumlike approximation 89 and density functional theory combined with pseudopotentials. 90,91However, for most of these models to perform properly, some form of system-specific parametrization is typically necessary, such as the construction of effective potentials.On the other hand, all-electron quantum chemical methods are much less system-dependent and offer a more universal simulation framework, which may prove very useful in possible future studies, such as investigating the role of various chemical modifications on HHG in nanostructures.
Therefore, the first objective of this work is to assess whether the applicability of RT-TDCIS can be extended to large structures containing tens of atoms.RT-TDCIS, being an equivalent of the Hartree-Fock (HF) method for excited states, does not account for correlation effects.Although it has been demonstrated that dynamical electron correlation has little effect on the laser-driven dynamics in atoms 66,[92][93][94][95][96] and small molecules, 70,71,95,97,98 in the case of C 60 it is known that the single-determinant restricted HF (RHF) wavefunction is not a stable ground state due to global correlations in the π orbital space. 99Therefore, limiting the calculations solely to RT-TDCIS may not be sufficient to obtain reliable HHG spectra.1][102] While including correlation effects implicitly through the use of the exchange-correlation potential of choice, RT-TDDFT suffers from its own limitations, such as an inability to describe single excitations in closed-shell systems, 103,104 non-linearity of the time-evolution equations, 105 and, last but not least, computational costs greatly exceeding those of RT-TDCIS. 51The latter drawback is also shared with more sophisticated multideterminant methods such as real-time time-dependent coupled cluster (RT-TDCC), 106 RT-TDCISD 57,58,66,71,98 and RT-TDCIS(D). 54,57,58An alternative approach to RT-TDDFT, proposed by Pauletti et al. 70 and also utilized in this work, is to add the exchange-correlation potential directly to the RT-TDCIS Hamiltonian, effectively turning it into the real-time time-dependent counterpart of the Tamm-Dancoff approximation (TDA).Therefore, we compute the HHG spectra of the C 60 molecule using RT-TDCIS and RT-TDTDA, as well as their respective generalizations that include deexcitation terms: the real-time time-dependent extensions of the linear-response time-dependent Hartee-Fock method and the linear-response time-dependent density functional theory.This allows us to assess the role of multi-electron effects in the laser-driven dynamics in fullerenes, which constitutes the second goal of this paper.Finally, we also calculate the HHG response using the semiempirical INDO/S Hamiltonian, which has been parameterized to reproduce excitation energies obtained with CIS in a limited active space. 107INDO/S has been recently extended to the real-time domain, but so far, it has only been employed in the modeling of absorption spectra of in the perturbative regime, giving somewhat promising results. 108Thus, the third and final aim of this paper is to investigate if INDO/S can also be applied for simulating strong field dynamics in large systems and serve as a less expensive alternative to ab initio and DFT-based methods.
The paper is constructed as follows.In Sec.II we present a brief theoretical background of the used methods, as well as provide computational details of the simulations.In Sec.III we present and discuss the results of the HHG calculations on the C 60 fullerene.Finally, in Sec.IV we conclude our work.

Theory
In this section we provide an overview of the theoretical foundations for the methods used in the present work, collectively referred to as the real-time time-dependent single excitationbased methods (RT-TDSEM).The common feature among all RT-TDSEM is providing an approximate solution to the TDSE, by expanding the time-dependent electronic wavefunction in the basis of the time-independent eigenstates Ψ m of the examined system, which include the ground state Ψ 0 and the excited states Ψ k , In all considered RT-TDSEM Ψ 0 is assumed to be the ground state closed-shell Slater determinant Φ 0 built from real, orthonormal occupied molecular orbitals (MOs) ϕ i represented in the linear combination of atomic orbitals (LCAO) approximation, In most electronic structure methods, Gaussian functions are chosen for the atomic orbital basis set χ µ due to computational efficiency.The molecular orbitals ϕ i are the solutions to the Hartree-Fock or Kohn-Sham self-consistent field (SCF) equations, where F is either the Fock or Kohn-Sham matrix, S is the overlap matrix, c is the matrix of coefficients c µi , and E is the diagonal matrix of MO energies.In addition to the occupied MOs, solving the SCF problem also provides the set of virtual MOs ϕ a .
The excited states Ψ k in Eq. ( 2) are obtained by solving a linear-response equation specific to the particular RT-TDSEM.The arguably most general linear-response theory considered in this work is the linear-response time-dependent density functional theory (LR-TDDFT), based on the Kohn-Sham reference determinant.The (non-Hermitian) LR- Here, the matrices A and B are commonly referred to as the excitation and deexcitation matrix, respectively.Their elements are defined as where i and j denote occupied MOs (hole states) ϕ i and ϕ j , a and b denote virtual MOs (particle states) ϕ a and ϕ b , ϵ i and ϵ a are the occupied and virtual MO energies, respectively, and f xc is the exchange-correlation kernel used in the SCF procedure.The two-electron Coulomb integrals (ia|bj) and exchange-correlation integrals (ia|f xc |bj) are written in the Mulliken notation.According to Casida, 110 every excited state Ψ k with the excitation energy where Φ a i constitute the set of singly excited Slater determinants.The vectors X ν and Y ν also fulfill the normalization condition If the orbitals used to construct the matrices A and B are obtained from the Hartree-Fock method, the exchange-correlation kernel f xc is replaced by the (non-local) HF exchange kernel f HF , (ia|f xc |bj) = (ia|f HF |bj) = −(ib|aj), and Eq. ( 5) is reduced to the linearresponse time-dependent Hartree-Fock (LR-TDHF) method, also known as the random phase approximation (RPA). 109It is worth mentioning that despite their similarities, LR-TDHF and LR-TDDFT have, in fact, been derived independently from each other, with the former being significantly older. 111,112However, in this work it is more convenient for us to treat LR-TDHF as a special case of LR-TDDFT, since we aim at investigating the effect of the exchange-correlation potential on the laser-driven electron dynamics.
The terms LR-TDHF and LR-TDDFT have traditionally been used for the methods solely for obtaining the excited eigenspectra of atoms and molecules.On the other hand, RT-TDHF and RT-TDDFT usually refer to the single-determinant approaches in which the real-time propagation is applied to molecular orbitals.4][115] In order to avoid possible confusion, we make use of this fact in this work and refer to the method of propagating the time-dependent wavefunction (2), expanded in the basis of states obtained through the solution of Eq. ( 5), as RT-TDRPA/xc, where xc may stand for HF, DFT, or any particular exchange-correlation functional.
Setting the deexcitation matrix B to zero in the LR-TDDFT (LR-TDHF) equation leads to the well-known TDA (CIS) approximation, 109 The properties (8) and ( 9) also apply to TDA and CIS states, with the exception that Y k = 0 for every Ψ k .TDA (CIS) approximation already provides a significant simplification of the linear-response problem compared to the full LR-TDDFT (LR-TDHF), as Eq.(10)   is Hermitian.This allows for a reduction in computational costs required to obtain excited states, as well as helps avoid various numerical issues, with the most infamous being the triplet instabilities. 116However, TDA (CIS) still necessitates calculations of all the oneand two-electron integrals required to solve the ground state Kohn-Sham (Hartree-Fock) equations and to construct the excitation matrix A, which can be prohibitive for very large systems.
To overcome this bottleneck, various semiempirical quantum chemical methods have been historically proposed, all of which rely on setting certain types of integrals to zero or replacing them with much simpler parameterized formulae.Although the parameters of most semiempirical methods have been adjusted to reproduce the results of ab initio ground state calculations, 117 one of them, named INDO/S, has been specifically constructed to yield excitation energies matching those of CIS in a truncated active orbital space. 107,118e detailed description of INDO/S can be found elsewhere, 107,[119][120][121] while here we briefly discuss only the key assumptions of this method.Unlike in the ab initio and DFT-based methods, and similarly to other semiempirical methods, the MOs in INDO/S are expanded not in a Gaussian basis set, but in the minimal basis set of Slater orbitals that describe only the valence orbitals of every element.INDO/S originates from the Hartree-Fock formalism, with the elements of the Fock matrix defined in the atomic orbital representation as where h µν are the elements of the one-electron Hamiltonian matrix and P λσ are the elements of the one-electron density matrix, P λσ = 2 occ i c λi c σi However, the solution of the SCF problem is heavily simplified due to the so-called zero differential overlap (ZDO) approximation, which sets the overlap matrix S in Eq. ( 4) to an identity matrix, S µν = δ µν .As a consequence, all three-and four-center two-electron integrals between basis functions (µν|λσ) vanish.The two-center two-electron are also set equal to zero with the exception of Coulomb-

like integrals involving only two basis functions, (µ
which are calculated from the Mataga-Nishimoto formula 122 (the superscripts A and B denote different atomic centers).The one-center two-electron integrals are estimated using the Slater-Condon Coulomb and exchange factors.The one-electron integrals h µν are also approximated using combinations of one-electron core integrals, resonance integrals, modified overlap integrals, and two-center two-electron integrals. 107All approximations used during the construction of the semiempirical Fock matrix are also applied when constructing the INDO/S excitation matrix A.
Having obtained the set of electronic states Ψ k using the linear-response theory of choice, one can proceed to propagate the time-dependent wavefunction Ψ(t).The time-dependent Hamiltonian Ĥ(t) in Eq. ( 1) consists of the ground state molecular Hamiltonian Ĥ0 and the interaction operator coupling the molecule to the external laser field, represented in the length gauge and within the dipole approximation, where ⃗ μ is the molecular dipole operator and ⃗ E(t) is the time-dependent external electric field.Inserting Eqs. ( 2) and ( 12) into Eq.( 1) leads to the equation for the time-evolution of the state vector Here, H 0,mn = δ mn E m and µ α is the matrix of the α-th spatial component of the dipole moment operator, µ α,mn = ⟨Ψ m |µ α Ψ n ⟩.The elements of µ α , in particular the dipole moment expectation values of excited states µ α,kk and the dipole transition moments between two excited states µ α,kl , can be determined with high accuracy using quadratic response theory. 123However, since we aim at propagating the time-dependent wavefunction using the full eigenspectrum of Eq. ( 5), which for systems as large as C 60 may consist of tens of thousands of states, this approach would be practically unfeasible.Therefore, we calculate the elements of µ α in an approximate manner, by inserting Eq. ( 8) directly into ⟨Ψ m |µ α Ψ n ⟩: 124 The dipole moment integrals between the ground and excited Slater determinants can be readily evaluated using Slater-Condon rules.For CIS, TDA, and INDO/S states these expressions reduce to To solve Eq. ( 13) we introduce time discretization and propagate the wavefunction using the second-order split-operator technique, where ∆t denotes the time step and the unitary matrix Similarly to our previous works 65,66,71,76 and to the works of other authors utilizing RT-TDSEM to simulate strong field electron dynamics, [59][60][61]63,69,70,[72][73][74][75] we employ the heuristic finite lifetime model of Klinkusch et al. 125 to compensate for the incompleteness of the atomic orbital basis sets. The electronicenergies E k (the diagonal entries of H 0 ) of excited states beyond the ionization threshold are modified by adding imaginary ionization rates, The finite lifetime model was originally developed for RT-TDCIS, 125 and later extended to RT-TDCI with higher excitations. 66,126The ionization rates of CIS states are defined as where θ(x) is the Heaviside step function and the empirical parameter d represents a maximum distance from the molecule that a (semiclassical) electron can travel before undergoing ionization.Naturally, Eq. ( 20) applies to TDA and INDO/S states as well.In this work, we also extend the finite lifetime model to RT-TDRPA.By analogy with Eq. ( 20) we define the heuristic ionization rates of LR-TDDFT and LR-TDHF as The motivation behind Eq. ( 21) is as follows.The ionization rate of every CIS state (20)   can be interpreted as a sum of excitation probabilities to individual virtual orbitals multiplied by the ionization rates of these orbitals √ 2ϵ a /d (the Heaviside function assures that only the virtual MOs with positive energies are ionizable).Since the RPA theory accounts for both excitations and deexcitations, the ionization rates of LR-TDHF and LR-TDDFT have to be accordingly reduced by the probabilities |Y k ia | 2 that an electron becomes deexcited from the virtual MO ϕ a back to the occupied MO ϕ i .
Once the time-dependent wavefunction Ψ(t) is known, the HHG spectrum of the examined system I HHG can be calculated from the Fourier transform of the optical response, which in this work is taken to be the dipole acceleration ⟨ µ⟩, where T is the total propagation time.

Computational details
We perform the calculations of the HHG spectra of the C 60 molecule subjected to short, intense laser pulses at the RT-TDRPA/DFT, RT-TDRPA/HF, RT-TDTDA, RT-TDCIS, and RT-TDINDO/S levels of theory.The geometry of C 60 was optimized at the B3LYP-D3/cc-pVTZ level of theory using Gaussian16 (Rev.C.01) software. 127The RT-TDRPA/DFT and RT-TDTDA calculations are perfomed in two variants, utilizing two different exchangecorrelation functionals.The first one is the standard B3LYP hybrid functional, 128 frequently employed in calculations involving fullerenes and their derivatives.The second one is its Coulomb-attenuated version, CAM-B3LYP, 129 which is reported to perform better than B3LYP in describing excited states 130,131 -a feature that may be important for the correct description of HHG.Since we consider the time-evolution of a closed-shell system in the absence of spin-dependent perturbations, we are only interested in the singlet excited manifold, so the linear-response equations are constructed using singlet configuration state functions rather than pure Slater determinants.We obtain the A and B matrices and the dipole moment integrals used in the RT-TDRPA, RT-TDTDA and RT-TDCIS calculations using a modified version of PySCF 2.4 package. 1324][135] The solution of the linear-response equation, the construction and diagonalization of the dipole moment matrices µ, and the real-time propagation is performed using a home-made program.All utilized approaches employ a full diagonalization of either Eq.(10) or Eq. ( 5).The LR-TDDFT and LR-TDHF equations are solved using the Cholesky decomposition technique. 136 the simulations, the external laser field is represented by a linearly polarized electric field pulse with a sine-squared envelope: The cycle-averaged laser intensity I 0 , related to the field amplitude E 0 via Also, the value of I p used in the finite lifetime model is estimated from the Koopmans' theorem, as the negative energy of the highest occupied molecular orbital −ϵ HOMO within the reference state corresponding to a particular RT-TDSEM.
The HHG spectra are computed from the dipole acceleration obtained by numerically differentiating the time-resolved dipole moment twice with respect to time.We also apply the Hann window to the dipole acceleration before taking the Fourier transform in order to account for the finite simulation time.

Optimization of the HHG simulation framework
When performing quantum-chemical calculations in the strong-field non-perturbative regime, three factors have the most significant influence on the accuracy of reproducing laser-driven electron dynamics: (a) the employed Gaussian basis set.Given that multiple transitions to highly excited and unbound electronic states are an inherent part of HHG, it is desired to simulate the electron dynamics using the most accurate available representation of the electronic continuum.In smaller systems, this typically entails using large numbers of highly diffuse functions, often tailored for this specific purpose. 60,61,65,137It is worth noting that in our case, the choice of the basis set pertains only to the ab initio and DFT-based methods, as INDO/S has been designed to work solely with the minimal Slater basis set; (b) the size of the active orbital space.When the size of the simulated system becomes considerable, a full configurational space with excitations from all occupied MOs to all virtual MOS can no longer be employed in the calculations, especially if a large basis set is used and the propagation involves the entire eigenspectrum of the linear-response equation.In our earlier calculations on the H 2 molecule we demonstrated that truncating the virtual orbital space may be suboptimal, because excluding the highest lying virtual MOs may negatively affect the description of all excited states, not just those with the highest energies. 71At the same time, other studies indicate that HHG is usually dominated by transitions from several of the highest lying occupied MOs, 75 so excluding excitations from core orbitals may be a preferred option for reducing the number of configurations; (c) the parameterization of the applied wavefunction absorber.In the finite lifetime model employed in our calculations the absorption rate is governed by the escape length value d.
Ideally, the absorber should selectively eliminate components of the wavefunction that cannot be accurately represented by the basis set, without interfering with the HHG process.In practical applications involving atoms and moderately-sized molecules, known to generate harmonics in accordance with the three-step model, a value of d close to the maximum electron excursion distance in the laser field, E 0 /ω 2 0 , is usually selected.However, since HHG in nanostructures is known to deviate from the 3SM, the proper choice of d requires additional investigation.
Therefore, before comparing the HHG spectra obtained using different RT-TDSEM methods, we must first ensure that the results are converged with respect to all three above parameters.To achieve this, we conduct a series of benchmark HHG calculations for the C 60 molecule at both considered carrier frequencies, at the RT-TDTDA/B3LYP level of theory.
We test the performance of three basis sets: the minimal STO-3G basis set, the doublezeta Dunning cc-pVDZ basis set, and its singly augmented variant, aug-cc-pVDZ.Due to the icosahedral symmetry of C 60 , almost all of its MOs belong to degenerate energy levels, with the HOMO level exhibiting fivefold degeneracy (Fig. 1).We determine the optimal configurational space by employing a full virtual orbital space and gradually expanding the active occupied orbital space, starting from five HOMO orbitals and adding one occupied MO shell at a time.Finally, we perform every calculation using seven different values of the Specifically, it predicts an abrupt drop in the HHG intensity near the cutoff energy predicted from the three step model, E cut = I p + 3.17(E 2 0 /4ω 2 0 ), at both wavelengths.In contrast, both cc-pVDZ and aug-cc-pVDZ spectra show a sizeable cutoff extension, accompanied by a more gradual reduction in the HHG intensity relative to the harmonic order.This result is more in line with previous theoretical calculations of HHG in fullerenes and with the experimental observations.Interestingly, the inclusion of the diffuse basis functions has a noticeably less pronounced effect on the description of HHG compared to the addition of the polarization functions.This is evident as the spectra obtained in the cc-pVDZ and aug-cc-pVDZ basis sets do not differ significantly in terms of their overall shape.The only distinctions are that the latter basis yields slightly higher peak intensities in the high-energy part of the spectrum and is capable of describing several additional peaks.Since these differences are, nonetheless, noticeable, in the further calculations we use the largest aug-cc-pVDZ basis set.
The two middle plots in Fig. 2   At this point we would like to note that unlike in calculations for smaller systems, we do not supplement the employed basis sets with additional diffuse functions designed to describe the electron dynamics in the continuum.This is mainly dictated by practical considerations.The rank of the A and B matrices constructed in the aug-cc-pVDZ basis set, using the full virtual space and active occupied space with n o = 35, is equal to 42000.Fully diagonalizing even larger matrices, though technically feasible, would demand significant time and resources.To the best of our knowledge, this is already the highest number of excited states reported for use in RT-TDSEM-based calculations.Moreover, in case of the C 60 fullerene, the aug-cc-pVDZ basis set is already on the verge of linear dependencies.However, judging by the small differences between the cc-pVDZ and aug-cc-pVDZ spectra, we can infer that HHG in C 60 mainly occurs on the surface of the fullerene.Therefore, addition of more diffuse functions is not necessary.This is further evidenced by the high optimal d value that greatly exceeds the maximum electron excursion distance under the considered laser conditions (equal to 11.6 bohr at 800 nm and 30.7 bohr at 1300 nm), even when extended by the largest internuclear distance in C 60 (≈ 13.4 bohr).The orbital ionization rate in the finite lifetime model, √ 2ϵ a /d, is interpreted as an inverse of the time required for an electron with kinetic energy ϵ a to travel a distance d.In this context, the parameter d can be considered equivalent to the distance between the molecule and the starting point of the complex absorbing potential.However, from Fig. 2, it is evident that setting d similar to the electron excursion distance results in an overestimation of absorption.This overestimation occurs due to the assignment of excessively short lifetimes to virtual MOs, akin to placing the complex absorbing potential too close to the molecule.This effect implies that in C 60 , the proportionality between the virtual MO energy and the distance that an electron occupying this MO can travel within a given amount of time is no longer valid.The probable cause is that even virtual MOs with relatively high energies are still localized in the vicinity of the molecular surface.This indirectly confirms that fullerene HHG results from oscillations of the electron density within the molecule and cannot be fully described by the three-step model.

Comparison of different exchange-correlation potentials
In this section we compare the HHG spectra computed using the B3LYP functional with those obtained using the CAM-B3LYP functional and calculated at the HF level.At this point, we adhere to the CIS/TDA approximation.Let us briefly remind that while B3LYP and CAM-B3LYP share the same correlation functional, they differ in terms of the exchange functional.B3LYP contains a fixed portion of 20% HF exchange, which is known to lead to incorrect behavior of this functional at the long range.Therefore, CAM-B3LYP has been proposed as a range-separated version of B3LYP, in which the percentage of HF exchange depends on interelectronic distance and varies from 19% at the short range to 65% at the long range. 129While still not making it asymptotically correct, this considerably improves the long-range behavior of the functional.
The HHG spectra obtained from RT-TDCIS, RT-TDTDA/B3LYP, and RT-TDTDA/CAM-B3LYP are shown on Fig. 3.All three methods provide spectra of comparable quality in terms of overall shape and intensity.However, we can distinguish two regions where certain systematic differences can be observed depending on f xc used.The first of them is the lowest-energy part of the spectrum, where individual peaks are most clearly visible.It can be seen that RT-TDCIS predicts significantly higher intensities of peaks up to 20th harmonic order compared to both variants of RT-TDTDA, especially at 800 nm.In our opinion, this part of the spectrum is primarily infuenced by the short-range interactions, specifically the short-range correlation.The peaks with the lowest energy levels represent transitions to and from the least energetic excited states, during which the excited electron remains in close proximity to the molecule.Therefore, the dynamical correlation effects are expected be most prominent in this region.Since the correlation part of B3LYP is the same as of CAM-B3LYP, the peak intensities in this part of the spectrum predicted by these two functionals are much more similar to each other.We have already discussed the effect of decreased HHG intensity due to correlation effects in our works on smaller systems. 66,71The second region encompasses the highest-energy part of the spectrum, extending beyond the 40th harmonic at 800 nm and beyond the 80th harmonic at 1300 nm.By analogy, we suspect that the description of this region is mainly governed by the long-range interac-tions.It can be noticed that the intensity of the spectrum "tail" is more or less proportional to the amount of the HF exchange in the exchange-correlation functional.RT-TDCIS, which provides an asymptotically correct description of exchange effects, predicts the most intense peaks in this region.In contrast, RT-TDTDA/B3LYP, which is the least accurate at the long-range limit, predicts peaks of the lowest intensity.Additionally, at 1300 nm, both RT-TDCIS and RT-TDTDA/CAM-B3LYP are able to describe a few more peaks than RT-TDTDA/B3LYP.A similar pattern can also be seen when comparing the distributions of excited states obtained using different exchange-correlation potentials (Fig. 4).Introducing correlation by replacing the HF potential with B3LYP lowers nearly all excitation energies, but the correction of the long-range exchange has a counteracting effect, causing a slight shift back towards higher values.

Influence of deexcitation effects on the HHG response
In Fig. 5, we compare the HHG spectra obtained in the TDA/CIS approximation with those obtained using the full RPA framework, separately for each of the three considered correlation-exchange potentials.Including the B matrix when solving the linear-response problem has a much smaller impact on the description of HHG than the choice of the exchange-correlation potential.In all cases the spectra computed using RPA and TDA/CIS eigenstates nearly overlap with each other.The only systematic effect is that RT-TDRPA predicts a slightly lower HHG intensity across the entire spectrum, which is more visible at 800 nm.Also, this effect is subtly more pronounced for RT-TDRPA/DFT than for RT-TDRPA/HF.The only exception can be observed for CAM-B3LYP, where the backgroud level of the RT-TDRPA spectrum is upshifted compared to the RT-TDTDA spectrum, resulting in a reduction of the number of described peaks in the highest-energy region.
Interestingly, RT-TDRPA most significantly reduces the intensity of HHG in the lowestenergy region of the spectra, where the intensities of first few peaks are 2-4 times lower compared to RT-TDTDA and RT-TDCIS.Although we explained the similar differences between RT-TDTDA and RT-TDCIS spectra by the correlation effects, such an explanation may not be appropriate in this case.While RPA is recognized to be a correlated method for the ground state, it is a matter of debate whether this property applies to excited states as well.Recently, Berkelbach has shown that LR-TDHF is equivalent to a variant of EOM-CCD, in which the CCD ground state is taken as a reference state, but only single excitations are considered when solving the equations of motion. 138This implies that the RPA excitation energies do not include any additional correlation effects beyond those already accounted for in the exchange-correlation kernel, especially given that in all flavors of RT-TDRPA we use a single-reference ground state.We reach a similar conclusion when analysing the excited state distributions obtained from RPA and TDA/CIS calculations, which are almost identical to each other (Fig. 6).We believe that the differences between the peak intensities may instead come from the differences in the dipole moments and dipole transition moments between TDA/CIS and RPA states.It is known that LR-TDHF and LR-TDDFT obey the Thomas-Reiche-Kuhn sum rule, unlike CIS and TDA. 109Therefore, one should anticipate an improved description of the ground-to-excited transitions in the former two methods.On The insets show the lowest-energy regions of the spectra, plotted in linear scale the other hand, we calculate the excited-to-excited transition moments between LR-TDHF and LR-TDDFT using the approximate formula (15), whereas the analogous Eq. ( 17) is a valid expression for CIS and TDA eigenstates.Due to these two potential sources of errors, it is not possible to definitively determine which approach, with or without including the B matrix, allows for a better description of electron dynamics.When analysing the computed HHG signals, we noticed that some of the spectra calculated at 1300 nm display a local maximum of intensity between 15th and 30th harmonic order.At this wavelength, this region corresponds to photon energies of 14-29 eV.0][141][142][143] The occurence of such dipole resonances in known to enhance the harmonic intensity. 82Therefore, we take a closer look at this region of the spectra in Fig. 7. Surprisingly, the GDR is most clearly visible in the RT-TDCIS and RT-TDRPA/HF spectra, where a single peak at 25th harmonic order (corresponding to approximately 23.8 eV) is significantly enhanced.RT-TDTDA/CAM-B3LYP provides a not too dissimilar picture, as an enhancement of a group of peaks between 24th and 27th harmonic order (22.9-25.7 eV) can be observed.Both of these results are in reasonable agreement with experimental works that predict the GDR maximum at 20-22 eV for C 60 and 21-24 eV for C 60 + . 140,141,143RT-TDRPA/B3LYP predicts the strongest enhancement of 19th peak, corresponding to a somewhat lower energy of 18.1 eV.Finally, practically no enhancement can be seen in the RT-TDTDA/B3LYP and RT-TDRPA/CAM-B3LYP spectra.
This simple test allows us to draw a conclusion that the full RPA treatment of the excited states may indeed result in a less accurate description of the electron dynamics in C 60 , at least when a Kohn-Sham determinant is employed as a reference, as evidenced by the case of CAM-B3LYP calculations.Furthermore, the fact that B3LYP is unable to predict the GDR at the RT-TDTDA level and provides its incorrect position at the RT-TDRPA level may indicate that long-range exchange effects are particularly important for describing the global resonances in C 60 .Unfortunately, we are unable to detect the GDR-related enhancement at 800 nm, as at this wavelength it should be located at lower harmonic orders, where the overall HHG intensity is several orders of magnitude higher.provides at higher energies closely resembles that of RT-TDTDA/B3LYP in the minimal basis set (shown in Fig. 2), with a sudden decrease in the HHG intensity at too low energies.

Semiempirical HHG calculations
Nevertheless, before we deem RT-TDINDO/S unsuitable for HHG modeling, we must recall that the INDO/S Hamiltonian was parameterized to replicate CIS excitation energies in an already truncated active orbital space.Therefore, further constraining the active occupied space of INDO/S itself may yield suboptimal results.Having that in mind, we perform a second series of calculations, this time allowing excitations from all 120 INDO/S occupied MOs (light-blue curves on Fig. 8).This hugely improves the HHG depiction in the high-energy range.The positions of the last described peaks are upshifted from 37th to 57th harmonic order at 800 nm and from 57th to 87th harmonic order at 1300 nm.Additionally, the overall intensity of HHG now aligns more closely with that predicted by RT-TDCIS over a much broader energy range.Interestigly, RT-TDINO/S performs notably better at 1300 nm, where a good agreement with the RT-TDCIS curve is observed up to approximately 50th harmonic order and the intensities of the lowest-energy peaks are very closely reproduced.At 800 nm, RT-TDINDO/S predicts two regions where the intensity of high-harmonic generation (HHG) significantly decreases, the first one around 23th harmonic order and the second around 43rd harmonic order, leading to larger discrepancies with the RT-TDCIS predictions.
The above results give prompt a question about whether the observed improvement of RT-TDINDO/S performance is genuinely attributable to the effectiveness of the semiempirical Hamiltonian, or is merely a consequence of expanding the configurational space.To address it, we conduct analogous calculations at the RT-TDCIS using the STO-3G basis set, which  The superior performance of RT-TDINDO/S is also evident when comparing the excitation energies obtained from the respective linear-response equations (Fig. 9a).The distribution of excited states obtained using INDO/S, particularly those below 20 eV, relatively closely matches the one calculated with CIS in the Dunning basis set.In contrast, CIS in the minimal basis set significantly overestimates the excitation energies.It is worth highlighting that such a notable discrepancy can be observed despite the fact that exactly the same number of excited states is available in both Slater and Gaussian minimal basis set, totaling 14400.These disparities directly impact the time-resolved optical response, as seen on Fig. 9b.RT-TDINDO/S reasonably reproduces the oscillations in the dipole acceleration, with a moderate overestimation of amplitudes but accurate preservation of their temporal positions.This overestimation can be attributed to the underestimation of the excitation energies below ∼ 15 eV, resulting in an overpopulation of the respective excited states.On the other hand, the curve corresponding to CIS in the STO-3G basis set exhibits a smoother, less oscillatory character due to the limited availability of energetically accessible electronic transitions states.This translates to a reduction in the number of harmonic peaks in the HHG spectrum.

Conclusion
In this study we conducted and analyzed a series of calculations of the HHG spectra of the C 60 fullerene in the non-perturbational regime, employing various quantum chemical methods.These include the real-time time-dependent counterparts of CIS, TDA and RPA based on either the Hartree-Fock or Kohn-Sham determinant, as well as the semiempirical INDO/S method.
The main conclusion that can be drawn from this work is that ab initio and DFT-based methods coupled to Gaussian bases can be successfully applied to model strong-field processes in nanoscale systems.Naturally, these calculations are much more resource-intensive Nevertheless, owing to the linearity of the real-time propagation equations, the computational costs remain reasonable, and represent a modest price for an effectively all-electron picture of the laser-driven dynamics.Our calculations on C 60 correctly prectict that the high harmonic generation in fullerenes primarily arises from the oscillations of the electronic denisty at the molecular surface, with a predominant contribution from the 60 π electrons to this process.We also demonstrated that RT-TDSEM are able to detect more subtle features of the attosecond processes, such as the giant dipole resonance.
Our results indicate that the quantum chemical description of the electron dynamics primarily depends on the chosen basis set and the active orbital space employed.The influence of the exchange-correlation potential, on the other hand, is relatively minor.Nevertheless, we were able to identify some features of the spectra that can be attributed to both correlation and exchange effects.Based on our findings, we can recommend the use of range-separated DFT functionals in the calculations involving systems with a significant role of dynamical correlation.In our case, CAM-B3LYP successfully combines the good approximation to the short-range correlation of B3LYP with the proper treatment of the long-range exchange characteristic to HF-based methods.On the other hand, a full RPA description of the excited electronic states provides practically no improvement over CIS, and may, in fact, lead to inferior results compared to TDA based on the Kohn-Sham reference.
Finally, we have demonstrated that the semiempricial INDO/S Hamiltonian, traditionally recognized for offering reasonably accurate approximations to the excitation energies, can also be successfully used for modeling strong-field processes in large systems.The performance of RT-TDINDO/S surpasses that of RT-TDCIS and RT-TDA in a basis set of comparable size, yielding results in qualitative agreement with all-electron approaches employing larger basis sets, particularly at lower harmonic orders and for relatively long laser wavelengths.This positions it as a potentially valuable simulation framework for studying even larger systems, which are gaining increasing interest in attosecond science but exceed the capabilities of the ab ab initio and DFT-based methods.The applicability of RT-TDINDO/S will be further explored in future works.

is set to 5 ×
10 13 W/cm 2 , and the number of optical cycles n c = 10.We use two values of carrier frequency ω 0 : 1.55 eV and 0.95 eV, corresponding to the wavelengths of 800 nm and 1300 nm, respectively.The 800 nm pulse has the duration of 1103 a.u.(≈ 26.7 fs), while the 1300 nm pulse has the duration of 1793 a.u.(≈ 43.4 fs).Due to the high symmetry of C 60 and computational limitations, we consider only one polarization vector parallel to one of the S 6 improper axis.The wavefunctions are propagated using a timestep ∆t = 0.01 a.u., a value that ensured convergence of obtained results.Every propagation starts from the singledeterminant ground state, which serves as a reference for the corresponding RT-TDSEM.

d
parameter: 10, 50, 100, 150, 200, 250 and 300 bohr, and with the finite lifetime model turned off (which is equivalent to setting d → ∞).The results of these preliminary calculations are presented in Fig.2.The two top plots show the comparison of three considered basis sets.It can be seen that the STO-3G basis set provides a significantly different depiction of HHG in comparison to the Dunning basis sets.

Figure 1 :
Figure 1: MO energy diagram of 35 highest occupied (thick bars) and 6 lowest unnocupied (thin bars) molecular orbitals of the C 60 molecule, calculated at B3LYP/aug-cc-pVDZ level of theory.The irreps of the I h symmetry group corresponding to individual orbital groups are also provided depict the influence of the number of the active occupied MOs n o on the obtained HHG spectra.When increasing the active occupied space, we were able to achieve convergence at approximately n o ≈ 30, and, as can be seen, at both wavelengths the spectra obtained with n o = 31 and n o = 35 are nearly identical to each other.Such an outcome is consistent with the somewhat intuitive reasoning that the π band, consisting of 60 p-type orbitals of sp 2 -hybridized carbon atoms, should be most strongly involved in the HHG process because the electrons occupying it can move freely throughout the entire molecule.However, to ensure that the computed spectra are truly converged with respect to n o , in the subsequent calculations we use the highest considered number of 35 occupied MOs.In the plots, we also show the results for n o = 5, which in our case is the closest to the SAE model.Although not dramatically different from the converged ones, the spectra obtained with n o = 5 are of noticeably inferior quality, with the decreased HHG intensity and artificial local minima in the high-energy region.While some of the earliest HHG calculations on fullerenes based on the strong field approximation considered only electronic transitions from the HOMO level,83,84 our results indicate that this may not be the optimal strategy.Finally, in the two bottom plots we present the dependence between the value of the d parameter in the finite lifetime model and the computed HHG signal.It is evident that applying the two lowest d values leads to a significant reduction in HHG intensity and a decrease in the number of peaks, indicating that the absorption model interferes with the HHG process.The calculated HHG response stabilizes at much higher values of d ≈ 150

Figure 2 :
Figure 2: Results of the benchmark calculations of HHG in C 60 at the RT-TDTDA/B3LYP level of theory.Top row: comparison of the HHG spectra obtained using three different basis sets, with n o = 35 and d = 200 a.u.Black vertical arrows denote the positions of the HHG cutoff as predicted by the 3SM.Middle row: comparison of the HHG spectra obtained in the aug-cc-pVDZ basis set using different numbers of active occupied MOs, with d = 200.For better clarity, only the results for n o = 5, 31, and 35 are shown.Bottom row: comparison comparison of the HHG spectra obtained in the aug-cc-pVDZ basis set using different values of the d parameter, with n o = 35.

Figure 3 :
Figure 3: Comparison of the HHG spectra of C 60 calculated in the aug-cc-pVDZ basis set at the RT-TDCIS, RT-TDTDA/B3LYP, and RT-TDTDA/CAM-B3LYP level of theory.The insets show the lowest-energy regions of the spectra, plotted in linear scale

Figure 4 :
Figure 4: Distributions of the excitation energies obtained using CIS, TDA/B3LYP and TDA/CAM-B3LYP in the aug-cc-pVDZ basis set

Figure 6 :
Figure 6: Distributions of the excitation energies obtained using TDA/CAM-B3LYP and RPA/CAM-B3LYP in the aug-cc-pVDZ basis set.The inset shows the zoom of 600 lowest excitation energies

Finally, in this
section we analyze the results obtained from the RT-TDINDO/S propagations.As mentioned earlier, in RT-TDINDO/S simulations we employ the Slater minimal basis set compatible with semiempirical methods, while the size of the active occupied space and the value of the d parameter remains the same as in the ab initio and DFT-based calculations.Since INDO/S can be considered a semiempirical counterpart of CIS, a natural choice is to compare the computed spectra with RT-TDCIS ones.Such a comparison is shown in two upper plots in Fig. 8, where the purple curves represent RT-TDCIS results and the green curves represent RT-TDINDO/S results.RT-TDINDO/S successfully reproduces most of the peaks in the low-energy portion of the spectra.This is a commendable result, considering the approximate nature of the INDO/S Hamiltonian.However, the picture it

Figure 7 :
Figure 7: Excerpts of the spectra computed at 1300 nm from Fig. 5, covering the region in which the giant dipole resonance is expected to occur is our closest analogue to the INDO/S Slater minimal basis set.Here, we also use two different active occupied spaces, with n o = 35 and n o = 120.The results are presented on the two bottom plots of Fig. 8.While an increase in n o modestly extends the cutoff at both intensities, the magnitude of this extension is nowhere near what is observed in the RT-TDINDO/S spectra.The position of the last described peak is shifted upwards by about 5 harmonic orders at 800 nm and by about 12 harmonic orders at 1300 nm.

Figure 8 :
Figure 8: Top row: comparison of the HHG spectra computed using RT-TDCIS in the augcc-pVDZ basis set, and using RT-TDINDO/S within two different active occupied spaces.Bottom row: comparison of the HHG spectra computed using RT-TDCIS in the aug-cc-pVDZ basis set, and in the STO-3G minimal basis set within two different active occupied spaces