# Modeling of high harmonic generation in the C60 fullerene using ab initio,
DFT-based and semiempirical methods

Aleksander P. Woźniak<EMAIL_ADDRESS>Robert Moszyński  [

###### Abstract

We report calculations of the high harmonic generation spectra of the C60
fullerene molecule, employing a diverse set of real-time time-dependent
quantum chemical methods. All methodologies involve expanding the propagated
electronic wavefunction 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 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.

###### keywords:

American Chemical Society, LaTeX

Unknown University] Faculty of Chemistry, University of Warsaw, Pasteura 1,
02-093 Warsaw, Poland IR,NMR,UV

##  1 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, 2. It has enabled
the routine production of ultrashort coherent electromagnetic pulses, making
it possible to explore electron dynamics on previously unattainable
timescales. Attosecond impulses have a wide range of applications, including,
among others, molecular imaging 3, 4, 5, monitoring chemical and photochemical
reactions in real time 6, 7, determining molecular structures 8, 9, studying
photoionization 10, 11, and investigating quantum coherence and electronic
correlation effects 12, 13, 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 remaining the most commonly
used sources of harmonic radiation. In simple mono- or few-atomic 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 35. Thus, 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, 38.

Progress in experimental discoveries in attophysics also necessitates the
development of new theoretical methods capable of describing ultrafast
electron dynamics in increasingly larger systems. In the past decade, there
has been a significant rise in the popularity of employing quantum chemistry
methods extended to the real-time domain for this purpose 48, 49, 50, 51, 52.
These approaches are characterized by moderate computational costs typical of
their time-independent counterparts, as well as a reasonable level of
accuracy, allowing, for example, the consideration of multi-electron effects.
One of the most popular methods of this kind is the real-time time-dependent
configuration interaction singles (RT-TDCIS), in which the time-dependent
electronic wavefunction is represented as a linear combination of time-
independent ground and singly excited electronic eigenstates of the examined
system 53, 54, 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 73,
72, 76. This 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 C60 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, 81. This property is attributed to their high
polarizability, as well as the occurence of the plasmon resonance at the
fullerene surface 82. From 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, 91. However, 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 92, 93,
94, 95, 96, 66 and small molecules 97, 95, 98, 70, 71, in the case of C60 it
is known that the single-determinant restricted HF (RHF) wavefunction is not a
stable ground state due to global correlations in the $\pi$ orbital space 99.
Therefore, limiting the calculations solely to RT-TDCIS may not be sufficient
to obtain reliable HHG spectra. Among several quantum chemical approaches that
incorporate dynamical correlation into the real-time dynamics, the most
popular is the real-time time-dependent density functional theory (RT-TDDFT),
in which the propagation is applied directly to Kohn-Sham spinorbitals forming
the electron density 100, 101, 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 51. The 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, 98, 66, 71 and RT-
TDCIS(D)54, 57, 58. An 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 C60 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 107. INDO/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 108. Thus, 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 C60 fullerene. Finally, in Sec. IV we conclude our work.

##  2 Methods

###  2.1 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 excitation-based methods (RT-TDSEM). The common feature
among all RT-TDSEM is providing an approximate solution to the TDSE,

$i\frac{\partial}{\partial t}\Psi(t)=\hat{H}(t)\Psi(t),$ (1)

by expanding the time-dependent electronic wavefunction in the basis of the
time-independent eigenstates $\Psi_{m}$ of the examined system, which include
the ground state $\Psi_{0}$ and the excited states $\Psi_{k}$,

$\Psi(t)=\sum_{m}C_{m}(t)\Psi_{m}=C_{0}(t)\Psi_{0}+\sum_{k>0}C_{k}(t)\Psi_{k}.$
(2)

In all considered RT-TDSEM $\Psi_{0}$ is assumed to be the ground state
closed-shell Slater determinant $\Phi_{0}$ built from real, orthonormal
occupied molecular orbitals (MOs) $\phi_{i}$ represented in the linear
combination of atomic orbitals (LCAO) approximation,

$\phi_{i}=\sum_{\mu}c_{\mu i}\chi_{\mu}$ (3)

In most electronic structure methods, Gaussian functions are chosen for the
atomic orbital basis set $\chi_{\mu}$ due to computational efficiency. The
molecular orbitals $\phi_{i}$ are the solutions to the Hartree-Fock or Kohn-
Sham self-consistent field (SCF) equations,

$\mathbf{Fc}=\mathbf{ScE}$ (4)

where $\mathbf{F}$ is either the Fock or Kohn-Sham matrix, $\mathbf{S}$ is the
overlap matrix, $\mathbf{c}$ is the matrix of coefficients $c_{\mu i}$, and
$\mathbf{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
$\phi_{a}$.

The excited states $\Psi_{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-TDDFT eigenequation reads 109

$\begin{pmatrix}\mathbf{A}&\mathbf{B}\\\ -\mathbf{B}&-\mathbf{A}\\\
\end{pmatrix}\begin{pmatrix}\mathbf{X}\\\ \mathbf{Y}\\\
\end{pmatrix}=\omega\begin{pmatrix}\mathbf{X}\\\ \mathbf{Y}\\\ \end{pmatrix}$
(5)

Here, the matrices $\mathbf{A}$ and $\mathbf{B}$ are commonly referred to as
the excitation and deexcitation matrix, respectively. Their elements are
defined as

$\displaystyle
A_{ia,jb}=\delta_{ij}\delta_{ab}(\epsilon_{a}-\epsilon_{i})+(ia|jb)+(ia|f_{xc}|jb)$
(6) $\displaystyle B_{ia,jb}=(ia|bj)+(ia|f_{xc}|bj)$ (7)

where $i$ and $j$ denote occupied MOs (hole states) $\phi_{i}$ and $\phi_{j}$,
$a$ and $b$ denote virtual MOs (particle states) $\phi_{a}$ and $\phi_{b}$,
$\epsilon_{i}$ and $\epsilon_{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 $\Psi_{k}$ with the excitation
energy $\omega_{k}=E_{k}-E_{0}$ is defined by the vectors $\mathbf{X}_{k}$ and
$\mathbf{Y}_{k}$ as

$\Psi_{k}=\sum_{i}^{occ}\sum_{a}^{vir}(X_{ia}^{k}+Y_{ia}^{k})\Phi_{i}^{a},$
(8)

where $\Phi_{i}^{a}$ constitute the set of singly excited Slater determinants.
The vectors $\mathbf{X}_{\nu}$ and $\mathbf{Y}_{\nu}$ also fulfill the
normalization condition

$\sum_{i}^{occ}\sum_{a}^{vir}(X_{ia}^{k}X_{ia}^{l}-Y_{ia}^{k}Y_{ia}^{l})=\delta_{kl}.$
(9)

If the orbitals used to construct the matrices $\mathbf{A}$ and $\mathbf{B}$
are obtained from the Hartree-Fock method, the exchange-correlation kernel
$f_{xc}$ is replaced by the (non-local) HF exchange kernel $f_{\mathrm{HF}}$,
$(ia|f_{xc}|bj)=(ia|f_{\mathrm{HF}}|bj)=-(ib|aj)$, and Eq. (5) is reduced to
the linear-response time-dependent Hartree-Fock (LR-TDHF) method, also known
as the random phase approximation (RPA) 109. It 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,
112. However, 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. Moreover, some authors use the term RPA interchangeably not only
with LR-TDHF, but also with LR-TDDFT 113, 114, 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 $\mathbf{B}$ to zero in the LR-TDDFT (LR-TDHF)
equation leads to the well-known TDA (CIS) approximation 109,

$\mathbf{A}\mathbf{X}=\omega\mathbf{X},$ (10)

The properties (8) and (9) also apply to TDA and CIS states, with the
exception that $\mathbf{Y}_{k}=0$ for every $\Psi_{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 116. However, TDA (CIS)
still necessitates calculations of all the one- and two-electron integrals
required to solve the ground state Kohn-Sham (Hartree-Fock) equations and to
construct the excitation matrix $\mathbf{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, 118. The
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

$F_{\mu\nu}=h_{\mu\nu}+\sum_{\lambda\sigma}P_{\lambda\sigma}\left((\mu\nu|\lambda\sigma)-\frac{1}{2}(\mu\nu|\lambda\sigma)\right),$
(11)

where $h_{\mu\nu}$ are the elements of the one-electron Hamiltonian matrix and
$P_{\lambda\sigma}$ are the elements of the one-electron density matrix,
$P_{\lambda\sigma}=2\sum_{i}^{occ}c_{\lambda i}c_{\sigma 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
$\mathbf{S}$ in Eq. (4) to an identity matrix, $S_{\mu\nu}=\delta_{\mu\nu}$.
As a consequence, all three- and four-center two-electron integrals between
basis functions $(\mu\nu|\lambda\sigma)$ vanish. The two-center two-electron
are also set equal to zero with the exception of Coulomb-like integrals
involving only two basis functions,
$(\mu^{A}\nu^{A}|\lambda^{B}\sigma^{B})=\delta_{\mu^{A}\nu^{A}}\delta_{\lambda^{B}\sigma^{B}}(\mu^{A}\mu^{A}|\lambda^{B}\lambda^{B})$,
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_{\mu\nu}$ are also approximated using
combinations of one-electron core integrals, resonance integrals, modified
overlap integrals, and two-center two-electron integrals 107. All
approximations used during the construction of the semiempirical Fock matrix
are also applied when constructing the INDO/S excitation matrix $\mathbf{A}$.

Having obtained the set of electronic states $\Psi_{k}$ using the linear-
response theory of choice, one can proceed to propagate the time-dependent
wavefunction $\Psi(t)$. The time-dependent Hamiltonian $\hat{H}(t)$ in Eq. (1)
consists of the ground state molecular Hamiltonian $\hat{H}_{0}$ and the
interaction operator coupling the molecule to the external laser field,
represented in the length gauge and within the dipole approximation,

$\hat{H}(t)=\hat{H}_{0}-\vec{\hat{\mu}}\vec{\mathcal{E}}(t),$ (12)

where $\vec{\hat{\mu}}$ is the molecular dipole operator and
$\vec{\mathcal{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 $\mathbf{C}(t)$ of time-dependent coefficients
$C_{k}(t)$,

$i\frac{\partial}{\partial
t}\mathbf{C}(t)=\left(\mathbf{H}_{0}+\sum_{\alpha=x,y,z}\mathcal{E}_{\alpha}(t)\boldsymbol{\mu}_{\alpha}\right)\mathbf{C}(t).$
(13)

Here, $\mathbf{H}_{0,mn}=\delta_{mn}E_{m}$ and $\boldsymbol{\mu}_{\alpha}$ is
the matrix of the $\alpha$-th spatial component of the dipole moment operator,
$\boldsymbol{\mu}_{\alpha,mn}=\langle\Psi_{m}|\mu_{\alpha}\Psi_{n}\rangle$.
The elements of $\boldsymbol{\mu}_{\alpha}$, in particular the dipole moment
expectation values of excited states $\boldsymbol{\mu}_{\alpha,kk}$ and the
dipole transition moments between two excited states
$\boldsymbol{\mu}_{\alpha,kl}$, can be determined with high accuracy using
quadratic response theory 123. However, since we aim at propagating the time-
dependent wavefunction using the full eigenspectrum of Eq. (5), which for
systems as large as C60 may consist of tens of thousands of states, this
approach would be practically unfeasible. Therefore, we calculate the elements
of $\boldsymbol{\mu}_{\alpha}$ in an approximate manner, by inserting Eq. (8)
directly into $\langle\Psi_{m}|\mu_{\alpha}\Psi_{n}\rangle$ 124:

$\displaystyle\langle\Psi_{0}|\mu_{\alpha}\Psi_{k}\rangle=\sum_{i}^{occ}\sum_{a}^{vir}\left(X_{ia}^{k}+Y_{ia}^{k}\right)\langle\Phi_{0}|\mu_{\alpha}\Phi_{i}^{a}\rangle,$
(14)
$\displaystyle\langle\Psi_{k}|\mu_{\alpha}\Psi_{l}\rangle=\sum_{ij}^{occ}\sum_{ab}^{vir}\left(X_{ia}^{k}X_{jb}^{l}+Y_{ia}^{k}Y_{jb}^{l}\right)\langle\Phi_{i}^{a}|\mu_{\alpha}\Phi_{j}^{b}.\rangle$
(15)

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

$\displaystyle\langle\Psi_{0}|\mu_{\alpha}\Psi_{k}\rangle=\sum_{i}^{occ}\sum_{a}^{vir}X_{ia}^{k}\langle\Phi_{0}|\mu_{\alpha}\Phi_{i}^{a}\rangle,$
(16)
$\displaystyle\langle\Psi_{k}|\mu_{\alpha}\Psi_{l}\rangle=\sum_{ij}^{occ}\sum_{ab}^{vir}X_{ia}^{k}X_{jb}^{l}\langle\Phi_{i}^{a}|\mu_{\alpha}\Phi_{j}^{b}.\rangle$
(17)

To solve Eq. (13) we introduce time discretization and propagate the
wavefunction using the second-order split-operator technique,

$\mathbf{C}(t+\Delta t)=e^{i\mathbf{H_{0}}\Delta
t/2}\left[\prod_{\alpha=x,y,z}\mathbf{U}_{\alpha}e^{i\mathcal{E}_{\alpha}(t)\mathbf{d}_{\alpha}\Delta
t}\mathbf{U}_{\alpha}^{\dagger}\right]e^{i\mathbf{H_{0}}\Delta
t/2}\mathbf{C}(t),$ (18)

where $\Delta t$ denotes the time step and the unitary matrix
$\mathbf{U}_{\alpha}$ diagonalizes the $\alpha$-th dipole component matrix,
$\mathbf{U}_{\alpha}\boldsymbol{\mu}_{\alpha}\mathbf{U}_{\alpha}^{\dagger}=\mathbf{d}_{\alpha}$.

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, 69, 72, 63, 73, 70, 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 electronic energies $E_{k}$ (the diagonal entries of
$\mathbf{H}_{0}$) of excited states beyond the ionization threshold are
modified by adding imaginary ionization rates,

$E_{k}\rightarrow\overline{E}_{k}=E_{k}-i\Gamma_{m}/2\qquad\text{for}\;E_{k}\geq
E_{0}+I_{p}.$ (19)

The finite lifetime model was originally developed for RT-TDCIS 125, and later
extended to RT-TDCI with higher excitations 126, 66. The ionization rates of
CIS states are defined as

$\Gamma_{k}=\sum_{i}^{occ.}\sum_{a}^{vir.}|X_{ia}^{k}|^{2}\theta(\epsilon_{a})\sqrt{2\epsilon_{a}}/d,$
(20)

where $\theta(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

$\Gamma_{k}=\sum_{i}^{occ.}\sum_{a}^{vir.}\left(|X_{ia}^{k}|^{2}-|Y_{ia}^{k}|^{2}\right)\theta(\epsilon_{a})\sqrt{2\epsilon_{a}}/d.$
(21)

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 $|X_{ia}^{k}|^{2}$ multiplied by the ionization
rates of these orbitals $\sqrt{2\epsilon_{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_{ia}^{k}|^{2}$ that an electron becomes deexcited from the
virtual MO $\phi_{a}$ back to the occupied MO $\phi_{i}$.

Once the time-dependent wavefunction $\Psi(t)$ is known, the HHG spectrum of
the examined system $I_{\mathrm{HHG}}$ can be calculated from the Fourier
transform of the optical response, which in this work is taken to be the
dipole acceleration $\ddot{\langle\mu\rangle}$,

$I_{\mathrm{HHG}}(\omega)=\left|\frac{1}{T}\int_{0}^{T}\left(\frac{d^{2}}{dt^{2}}\sum_{mn}C_{m}^{*}(t)C_{n}(t)\langle\Psi_{m}|\hat{\mu}\Psi_{n}\rangle\right)e^{i\omega
t}dt\right|^{2},$ (22)

where $T$ is the total propagation time.

###  2.2 Computational details

We perform the calculations of the HHG spectra of the C60 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 C60 was optimized at
the B3LYP-D3/cc-pVTZ level of theory using Gaussian16 (Rev. C.01) software
127. The RT-TDRPA/DFT and RT-TDTDA calculations are perfomed in two variants,
utilizing two different exchange-correlation 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 $\mathbf{A}$ and $\mathbf{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 132. The
semiempirical $\mathbf{A}$ matrices used in the RT-TDINDO/S calculations and
the dipole moment integrals between INDO/S configurations are generated using
a modified version of MOPAC22 133, 134, 135. The solution of the linear-
response equation, the construction and diagonalization of the dipole moment
matrices $\boldsymbol{\mu}$, 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.

In the simulations, the external laser field is represented by a linearly
polarized electric field pulse with a sine-squared envelope:

$\vec{\mathcal{E}}(t)=\begin{cases}\vec{\mathcal{E}}_{0}\sin(\omega_{0}t)\sin^{2}(\omega_{0}t/2n_{c})&\text{if}\;0\leq
t\leq 2\pi n_{c}/\omega_{0},\\\ 0&\text{otherwise.}\end{cases}$ (23)

The cycle-averaged laser intensity $I_{0}$, related to the field amplitude
$\mathcal{E}_{0}$ via $I_{0}=\mathcal{E}_{0}^{2}/2$, is set to $5\times
10^{13}$ W/cm2, and the number of optical cycles $n_{c}=10$. We use two values
of carrier frequency $\omega_{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. ($\approx$ 26.7 fs), while the 1300 nm pulse has the
duration of 1793 a.u. ($\approx$ 43.4 fs). Due to the high symmetry of C60 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 $\Delta t=0.01$ a.u., a value that ensured convergence of obtained
results. Every propagation starts from the single-determinant ground state,
which serves as a reference for the corresponding RT-TDSEM. 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
$-\epsilon_{\mathrm{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.

##  3 Results and discussion

###  3.1 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, 137. It 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 H2 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 71. At 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}/\omega_{0}^{2}$, 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 C60 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 double-zeta Dunning cc-pVDZ basis set,
and its singly augmented variant, aug-cc-pVDZ. Due to the icosahedral symmetry
of C60, 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 $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\rightarrow\infty$).

Figure 1: MO energy diagram of 35 highest occupied (thick bars) and 6 lowest
unnocupied (thin bars) molecular orbitals of the C60 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

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. 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_{0}^{2}/4\omega_{0}^{2})$, 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.

Figure 2: Results of the benchmark calculations of HHG in C60 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$.

The two middle plots in Fig. 2 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}\approx 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
$\pi$ 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\approx$ 150 bohr. The spectra
obtained with $d=150$, 200, 250, and 300 bohr exhibit consistent peak shapes
and only marginal differences in overall intensity. Therefore, we pick an
intermediate value of $d=200$ bohr for further calculations. It is worth
mentioning that the spectra shapes, once converged with respect to $d$, are
not identical to the spectra shapes obtained without the finite lifetime
model. This implies that the absorber effectively elliminates the wavefunction
refelctions arising from the basis set incompleteness, as intended, without
hindering the HHG efficiency.

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 $\mathbf{A}$ and
$\mathbf{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 C60 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 C60 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 C60 ($\approx$ 13.4 bohr). The orbital
ionization rate in the finite lifetime model, $\sqrt{2\epsilon_{a}}/d$, is
interpreted as an inverse of the time required for an electron with kinetic
energy $\epsilon_{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 C60, 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.

###  3.2 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 129. While 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, 71.

Figure 3: Comparison of the HHG spectra of C60 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

The 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 interactions. 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.

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

###  3.3 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 $\mathbf{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.

Figure 5: Comparison of the HHG spectra of C60 calculated in the aug-cc-pVDZ
basis set: RT-TDCIS vs RT-TDRPA/HF (top row), RT-TDTDA/B3LYP vs RT-TDRPA/B3LYP
(middle row), and RT-TDTDA/CAM-B3LYP vs RT-TDRPA/CAM-B3LYP (bottom row). The
insets show the lowest-energy regions of the spectra, plotted in linear scale

Interestingly, RT-TDRPA most significantly reduces the intensity of HHG in the
lowest-energy 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 138.
This 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 109. Therefore, one should anticipate an improved
description of the ground-to-excited transitions in the former two methods. On
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 $\mathbf{B}$ matrix, allows for
a better description of electron dynamics.

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

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. This is the energy range in which the giant dipole
resonance (GDR) can be observed in the photoionization spectra of the C60
fullerene 139, 140, 141, 142, 143. The occurence of such dipole resonances in
known to enhance the harmonic intensity 82. Therefore, 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
C60 and 21-24 eV for C60+ 140, 141, 143. RT-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 C60, 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 C60. 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.

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

###  3.4 Semiempirical HHG calculations

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 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. 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.

Figure 8: Top row: comparison of the HHG spectra computed using RT-TDCIS in
the aug-cc-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

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 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.

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 $\sim$ 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.

Figure 9: a) Distributions of excitation energies calculated with CIS in the
aug-cc-pVDZ and STO-3G basis sets and with INDO/S b) Time-resolved dipole
acceleration during the first five optical cycles of the applied laser pulse
of 1300 nm, obtained from the RT-TDCIS propagations in the aug-cc-pVDZ and
STO-3G basis sets and from the RT-TDINDO/S propagation. Calculation involving
CIS in the STO-3G basis set and INDO/S employed an increased active occupied
space with $n_{o}=120$

##  4 Conclusion

In this study we conducted and analyzed a series of calculations of the HHG
spectra of the C60 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 compared to those on atoms and
smaller molecules, as most of the performed propagations involved expanding
the time-dependent wavefunction in several tens of thousands of states.
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 C60 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 $\pi$
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.

This work was supported by the Polish National Science Centre (NCN) through
Grant No. 2017/25/B/ST4/02698. The calculations have been carried out using
resources provided by Wroclaw Centre for Networking and Supercomputing, Grant
No. 567, and by University of Arizona Research Computing, courtesy of Prof.
Ludwik Adamowicz.

## References

  * Agostini and DiMauro 2004 Agostini, P.; DiMauro, L. F. The physics of attosecond light pulses. _Rep. Prog. Phys._ 2004, _67_ , 813–855
  * Krausz and Ivanov 2009 Krausz, F.; Ivanov, M. Attosecond physics. _Rev. Mod. Phys._ 2009, _81_ , 163
  * Vozzi et al. 2011 Vozzi, C.; Negro, M.; Calegari, F.; Sansone, G.; Nisoli, M.; De Silvestri, S.; Stagira, S. Generalized molecular orbital tomography. _Nat. Phys._ 2011, _7_ , 822–826
  * Salières et al. 2012 Salières, P.; Maquet, A.; Haessler, S.; Caillat, J.; Taïeb, R. Imaging orbitals with attosecond and Ångström resolutions: toward attochemistry? _Rep. Prog. Phys._ 2012, _75_ , 062401
  * Peng et al. 2019 Peng, P.; Marceau, C.; Villeneuve, D. M. Attosecond imaging of molecules using high harmonic spectroscopy. _Nat. Rev. Phys._ 2019, _1_ , 144–155
  * Wörner et al. 2010 Wörner, H. J.; Bertrand, J. B.; Kartashov, D. V.; Corkum, P. B.; Villeneuve, D. M. Following a chemical reaction using high-harmonic interferometry. _Nature_ 2010, _466_ , 604–607
  * Baykusheva et al. 2019 Baykusheva, D.; Zindel, D.; Svoboda, V.; Bommeli, E.; Ochsner, M.; Tehlar, A.; Wörner, H. J. Real-time probing of chirality during a chemical reaction. _Proc. Natl. Acad. Sci. U.S.A_ 2019, _116_ , 23923–23929
  * Torres et al. 2010 Torres, R.; Siegel, T.; Brugnera, L.; Procino, I.; Underwood, J. G.; Altucci, C.; Velotta, R.; Springate, E.; Froud, C.; Turcu, I. C. E. et al. Revealing molecular structure and dynamics through high-order harmonic generation driven by mid-IR fields. _Phys. Rev. A_ 2010, _81_ , 051802(R)
  * Wong et al. 2011 Wong, M. C. H.; Brichta, J.-P.; Spanner, M.; Patchkovskii, S.; Bhardwaj, V. R. High-harmonic spectroscopy of molecular isomers. _Phys. Rev. A_ 2011, _84_ , 051403
  * Wörner et al. 2009 Wörner, H. J.; Niikura, H.; Bertrand, J. B.; Corkum, P. B.; Villeneuve, D. M. Observation of Electronic Structure Minima in High-Harmonic Generation. _Phys. Rev. Lett._ 2009, _102_ , 103901
  * Bruner et al. 2016 Bruner, B. D.; Mašín, Z.; Negro, M.; Morales, F.; Brambila, D.; Devetta, M.; Faccialà, D.; Harvey, A. G.; Ivanov, M.; Mairesse, Y. et al. Multidimensional high harmonic spectroscopy of polyatomic molecules: detecting sub-cycle laser-driven hole dynamics upon ionization in strong mid-IR laser fields. _Faraday Discuss._ 2016, _194_ , 369–405
  * Wörner et al. 2010 Wörner, H. J.; Bertrand, J. B.; Hockett, P.; Corkum, P. B.; Villeneuve, D. M. Controlling the Interference of Multiple Molecular Orbitals in High-Harmonic Generation. _Phys. Rev. Lett._ 2010, _104_ , 233904
  * Shiner et al. 2011 Shiner, A.; Schmidt, B.; Trallero-Herrero, C.; Wörner, H. J.; Patchkovskii, S.; Corkum, P. B.; Kieffer, J.; Légaré, F.; Villeneuve, D. Probing collective multi-electron dynamics in xenon with high-harmonic spectroscopy. _Nat. Phys._ 2011, _7_ , 464–467
  * Marciniak et al. 2019 Marciniak, A.; Despré, V.; Loriot, V.; Karras, G.; Hervé, M.; Quintard, L.; Catoire, F.; Joblin, C.; Constant, E.; Kuleff, A. I. et al. Electron correlation driven non-adiabatic relaxation in molecules excited by an ultrashort extreme ultraviolet pulse. _Nat. Commun._ 2019, _10_ , 337
  * Monfared et al. 2022 Monfared, M.; Irani, E.; Lemell, C.; Burgdörfer, J. Influence of coherent vibrational excitation on the high-order harmonic generation of diatomic molecules. _Phys. Rev. A_ 2022, _106_ , 053108
  * McPherson et al. 1987 McPherson, A.; Gibson, G.; Jara, H.; Johann, U.; Luk, T. S.; McIntyre, I. A.; Boyer, K.; Rhodes, C. K. Studies of multiphoton production of vacuum-ultraviolet radiation in the rare gases. _J. Opt. Soc. Am. B_ 1987, _4_ , 595–601
  * Ferray et al. 1988 Ferray, M.; L’Huillier, A.; Li, X. F.; Lompre, L. A.; Mainfray, G.; Manus, C. Multiple-harmonic conversion of 1064 nm radiation in rare gases. _J. Phys. B: At. Mol. Opt. Phys._ 1988, _21_ , L31
  * Li et al. 1989 Li, X. F.; L’Huillier, A.; Ferray, M.; Lompré, L. A.; Mainfray, G. Multiple-harmonic generation in rare gases at high laser intensity. _Phys. Rev. A_ 1989, _39_ , 5751–5761
  * Sarukura et al. 1991 Sarukura, N.; Hata, K.; Adachi, T.; Nodomi, R.; Watanabe, M.; Watanabe, S. Coherent soft-x-ray generation by the harmonics of an ultrahigh-power KrF laser. _Phys. Rev. A_ 1991, _43_ , 1669–1672
  * Crane et al. 1992 Crane, J. K.; Perry, M. D.; Herman, S.; Falcone, R. W. High-field harmonic generation in helium. _Opt. Lett._ 1992, _17_ , 1256–1258
  * Faldon et al. 1992 Faldon, M. E.; Hutchinson, M.; Marangos, J.; Muffett, J.; Smith, R.; Tisch, J.; Wahlstrom, C. Studies of time-resolved harmonic generation in intense laser fields in xenon. _J. Opt. Soc. Am. B_ 1992, _9_ , 2094–2099
  * Miyazaki and Sakai 1992 Miyazaki, K.; Sakai, H. High-order harmonic generation in rare gases with intense subpicosecond dye laser pulses. _J. Phys. B: At. Mol. Opt. Phys._ 1992, _25_ , L83
  * L’Huillier and Balcou 1993 L’Huillier, A.; Balcou, P. High-order harmonic generation in rare gases with a 1-ps 1053-nm laser. _Phys. Rev. Lett._ 1993, _70_ , 774–777
  * Kondo et al. 1993 Kondo, K.; Sarukura, N.; Sajiki, K.; Watanabe, S. High-order harmonic generation by ultrashort KrF and Ti:sapphire lasers. _Phys. Rev. A_ 1993, _47_ , R2480–R2483
  * Macklin et al. 1993 Macklin, J. J.; Kmetec, J. D.; Gordon, C. L. High-order harmonic generation using intense femtosecond pulses. _Phys. Rev. Lett._ 1993, _70_ , 766–769
  * Wahlström et al. 1993 Wahlström, C.-G.; Larsson, J.; Persson, A.; Starczewski, T.; Svanberg, S.; Salières, P.; Balcou, P.; L’Huillier, A. High-order harmonic generation in rare gases with an intense short-pulse laser. _Phys. Rev. A_ 1993, _48_ , 4709–4720
  * Liang et al. 1994 Liang, Y.; Augst, S.; Chin, S.; Beaudoin, Y.; Chaker, M. High harmonic generation in atomic and diatomic molecular gases using intense picosecond laser pulses-a comparison. _J. Phys. B: At. Mol. Opt. Phys._ 1994, _27_ , 5119
  * Chin et al. 1995 Chin, S.; Liang, Y.; Augst, S.; Golovinski, P.; Beaudoin, Y.; Chaker, M. High Harmonic Generation in Atoms and Diatomic Molecules Using Ultrashort Laser Pulses in the Multiphoton Regime. _J. Nonlinear Opt. Phys. Mater._ 1995, _4_ , 667–686
  * Lyngå et al. 1996 Lyngå, C.; L’Huillier, A.; Wahlström, C.-G. High-order harmonic generation in molecular gases. _J. Phys. B: At. Mol. Opt. Phys._ 1996, _29_ , 3293
  * Krause et al. 1992 Krause, J. L.; Schafer, K. J.; Kulander, K. C. Calculation of photoemission from atoms subject to intense laser fields. _Phys. Rev. A_ 1992, _45_ , 4998
  * Krause et al. 1992 Krause, J. L.; Schafer, K. J.; Kulander, K. C. High-order harmonic generation from atoms and ions in the high intensity regime. _Phys. Rev. Lett._ 1992, _68_ , 3535–3538
  * Schafer et al. 1993 Schafer, K. J.; Yang, B.; DiMauro, L. F.; Kulander, K. C. Above threshold ionization beyond the high harmonic cutoff. _Phys. Rev. Lett._ 1993, _70_ , 1599
  * Corkum 1993 Corkum, P. B. Plasma perspective on strong field multiphoton ionization. _Phys. Rev. Lett._ 1993, _71_ , 1994–1997
  * Lewenstein et al. 1994 Lewenstein, M.; Balcou, P.; Ivanov, M. Y.; L’Huillier, A.; Corkum, P. B. Theory of high-harmonic generation by low-frequency laser fields. _Phys. Rev. A_ 1994, _49_ , 2117–2132
  * Constant et al. 1999 Constant, E.; Garzella, D.; Breger, P.; Mével, E.; Dorrer, C.; Le Blanc, C.; Salin, F.; Agostini, P. Optimizing High Harmonic Generation in Absorbing Gases: Model and Experiment. _Phys. Rev. Lett._ 1999, _82_ , 1668–1671
  * Ghimire et al. 2011 Ghimire, S.; DiChiara, A. D.; Sistrunk, E.; Agostini, P.; DiMauro, L. F.; Reis, D. A. Observation of high-order harmonic generation in a bulk crystal. _Nat. Phys._ 2011, _7_ , 138–141
  * Schubert et al. 2014 Schubert, O.; Hohenleutner, M.; Langer, F.; Urbanek, B.; Lange, C.; Huttner, U.; Golde, D.; Meier, T.; Kira, M.; Koch, S. W. et al. Sub-cycle control of terahertz high-harmonic generation by dynamical Bloch oscillations. _Nat. Photonics_ 2014, _8_ , 119–123
  * Luu et al. 2015 Luu, T. T.; Garg, M.; Kruchinin, S. Y.; Moulet, A.; Hassan, M. T.; Goulielmakis, E. Extreme ultraviolet high-harmonic spectroscopy of solids. _Nature_ 2015, _521_ , 498–502
  * Ndabashimiye et al. 2016 Ndabashimiye, G.; Ghimire, S.; Wu, M.; Browne, D. A.; Schafer, K. J.; Gaarde, M. B.; Reis, D. A. Solid-state harmonics beyond the atomic limit. _Nature_ 2016, _534_ , 520–523
  * You et al. 2017 You, Y. S.; Reis, D. A.; Ghimire, S. Anisotropic high-harmonic generation in bulk crystals. _Nat. Phys._ 2017, _13_ , 345–349
  * Luu et al. 2018 Luu, T. T.; Yin, Z.; Jain, A.; Gaumnitz, T.; Pertot, Y.; Ma, J.; Wörner, H. J. Extreme–ultraviolet high–harmonic generation in liquids. _Nat. Commun._ 2018, _9_ , 3723
  * Golde et al. 2008 Golde, D.; Meier, T.; Koch, S. W. High harmonics generated in semiconductor nanostructures by the coupled dynamics of optical inter- and intraband excitations. _Phys. Rev. B_ 2008, _77_ , 075330
  * Singhal et al. 2010 Singhal, H.; Ganeev, R.; Naik, P.; Srivastava, A.; Singh, A.; Chari, R.; Khan, R.; Chakera, J.; Gupta, P. Study of high-order harmonic generation from nanoparticles. _J. Phys. B: At. Mol. Opt. Phys._ 2010, _43_ , 025603
  * Han et al. 2016 Han, S.; Kim, H.; Kim, Y. W.; Kim, Y.-J.; Kim, S.; Park, I.-Y.; Kim, S.-W. High-harmonic generation by field enhanced femtosecond pulses in metal-sapphire nanostructure. _Nat. Commun._ 2016, _7_ , 13105
  * Liu et al. 2017 Liu, H.; Li, Y.; You, Y. S.; Ghimire, S.; Heinz, T. F.; Reis, D. A. High-harmonic generation from an atomically thin semiconductor. _Nat. Phys._ 2017, _13_ , 262–265
  * Sivis et al. 2017 Sivis, M.; Taucer, M.; Vampa, G.; Johnston, K.; Staudte, A.; Naumov, A. Y.; Villeneuve, D.; Ropers, C.; Corkum, P. Tailored semiconductors for high-harmonic optoelectronics. _Science_ 2017, _357_ , 303–306
  * Shcherbakov et al. 2021 Shcherbakov, M. R.; Zhang, H.; Tripepi, M.; Sartorello, G.; Talisa, N.; AlShafey, A.; Fan, Z.; Twardowski, J.; Krivitsky, L. A.; Kuznetsov, A. I. et al. Generation of even and odd high harmonics in resonant metasurfaces using single and multiple ultra-intense laser pulses. _Nat. Commun._ 2021, _12_ , 4185
  * Lopata and Govind 2011 Lopata, K.; Govind, N. Modeling Fast Electron Dynamics with Real-Time Time-Dependent Density Functional Theory: Application to Small Molecules and Chromophores. _J. Chem. Theory Comput._ 2011, _7_ , 1344–1355
  * Ishikawa and Sato 2015 Ishikawa, K. L.; Sato, T. A review on ab initio approaches for multielectron dynamics. _IEEE J. Sel. Top. Quantum Electron._ 2015, _21_ , 1–16
  * Goings et al. 2017 Goings, J. J.; Lestrange, P. J.; Li, X. Real‐time time‐dependent electronic structure theory. _WIREs Comput. Mol. Sci._ 2017, _8_ , e1341
  * Bedurke et al. 2021 Bedurke, F.; Klamroth, T.; Saalfrank, P. Many-electron dynamics in laser-driven molecules: wavefunction theory vs. density functional theory. _Phys. Chem. Chem. Phys._ 2021, _23_ , 13544–13560
  * Coccia and Luppi 2022 Coccia, E.; Luppi, E. Time-dependent ab initio approaches for high-harmonic generation spectroscopy. _J. Phys. Condens. Matter_ 2022, _34_ , 073001
  * Klamroth 2003 Klamroth, T. Laser-driven electron transfer through metal-insulator-metal contacts: Time-dependent configuration interaction singles calculations for a jellium model. _Phys. Rev. B_ 2003, _68_ , 245421
  * Krause et al. 2005 Krause, P.; Klamroth, T.; Saalfrank, P. Time-dependent configuration-interaction calculations of laser-pulse-driven many-electron dynamics: Controlled dipole switching in lithium cyanide. _J. Chem. Phys._ 2005, _123_ , 074105
  * Huber and Klamroth 2005 Huber, C.; Klamroth, T. Simulation of two-photon-photoelectron spectra at a jellium-vacuum interface. _Appl. Phys. A_ 2005, _81_ , 93–101
  * Rohringer et al. 2006 Rohringer, N.; Gordon, A.; Santra, R. Configuration-interaction-based time-dependent orbital approach for ab initio treatment of electronic dynamics in a strong optical laser field. _Phys. Rev. A_ 2006, _74_ , 043420
  * Krause et al. 2007 Krause, P.; Klamroth, T.; Saalfrank, P. Molecular response properties from explicitly time-dependent configuration interaction methods. _J. Chem. Phys._ 2007, _127_ , 034107
  * Schlegel et al. 2007 Schlegel, H. B.; Smith, S. M.; Li, X. Electronic optical response of molecules in intense fields: Comparison of TD-HF, TD-CIS, and TD-CIS(D) approaches. _J. Chem. Phys._ 2007, _126_ , 244110
  * Luppi and Head-Gordon 2013 Luppi, E.; Head-Gordon, M. The role of Rydberg and continuum levels in computing high harmonic generation spectra of the hydrogen atom using time-dependent configuration interaction. _J. Chem. Phys._ 2013, _139_ , 164121
  * Coccia et al. 2016 Coccia, E.; Mussard, B.; Labeye, M.; Caillat, J.; Taïeb, R.; Toulouse, J.; Luppi, E. Gaussian continuum basis functions for calculating high-harmonic generation spectra. _Int. J. Quantum Chem._ 2016, _116_ , 1120–1131
  * Coccia and Luppi 2016 Coccia, E.; Luppi, E. Optimal-continuum and multicentered Gaussian basis sets for high-harmonic generation spectroscopy. _Theor. Chem. Acc._ 2016, _135_ , 43
  * Pabst et al. 2016 Pabst, S.; Sytcheva, A.; Geffert, O.; Santra, R. Stability of the time-dependent configuration-interaction-singles method in the attosecond and strong-field regimes: A study of basis sets and absorption methods. _Phys. Rev. A_ 2016, _94_ , 033421
  * Coccia and Luppi 2019 Coccia, E.; Luppi, E. Detecting the minimum in argon high-harmonic generation spectrum using Gaussian basis sets. _Theor. Chem. Acc._ 2019, _138_ , 96
  * Coccia 2020 Coccia, E. How electronic dephasing affects the high-harmonic generation in atoms. _Mol. Phys._ 2020, _118_ , e1769871
  * Woźniak et al. 2021 Woźniak, A. P.; Lesiuk, M.; Przybytek, M.; Efimov, D. K.; Prauzner-Bechcicki, J. S.; Mandrysz, M.; Ciappina, M.; Pisanty, E.; Zakrzewski, J.; Lewenstein, M. et al. A systematic construction of Gaussian basis sets for the description of laser field ionization and high-harmonic generation. _J. Chem. Phys._ 2021, _154_ , 094111
  * Woźniak et al. 2022 Woźniak, A. P.; Przybytek, M.; Lewenstein, M.; Moszyński, R. Effects of electronic correlation on the high harmonic generation in helium: A time-dependent configuration interaction singles vs time-dependent full configuration interaction study. _J. Chem. Phys._ 2022, _156_ , 174106
  * Luppi and Head-Gordon 2012 Luppi, E.; Head-Gordon, M. Computation of high-harmonic generation spectra of H2 And N2 in intense laser pulses using quantum chemistry methods and time-dependent density functional theory. _Mol. Phys._ 2012, _110_ , 909–923
  * White et al. 2016 White, A. F.; Heide, C. J.; Saalfrank, P.; Head-Gordon, M.; Luppi, E. Computation of high-harmonic generation spectra of the hydrogen molecule using time-dependent configuration-interaction. _Mol. Phys._ 2016, _114_ , 947–956
  * Labeye et al. 2018 Labeye, M.; Zapata, F.; Coccia, E.; Véniard, V.; Toulouse, J.; Caillat, J.; Taïeb, R.; Luppi, E. Optimal Basis Set for Electron Dynamics in Strong Laser Fields: The case of Molecular Ion H2+. _J. Chem. Theory Comput._ 2018, _14_ , 5846–5858
  * Pauletti et al. 2021 Pauletti, C. F.; Coccia, E.; Luppi, E. Role of exchange and correlation in high-harmonic generation spectra of H2, N2, and CO2: Real-time time-dependent electronic-structure approaches. _J. Chem. Phys._ 2021, _154_ , 014101
  * Woźniak et al. 2023 Woźniak, A. P.; Lewenstein, M.; Moszyński, R. In _Polish Quantum Chemistry from Kołos to Now_ ; Musiał, M., Grabowski, I., Eds.; Advances in Quantum Chemistry; Academic Press, 2023; Vol. 87; pp 167–190
  * Bedurke et al. 2019 Bedurke, F.; Klamroth, T.; Krause, P.; Saalfrank, P. Discriminating organic isomers by high harmonic generation: A time-dependent configuration interaction singles study. _J. Chem. Phys._ 2019, _150_ , 234114
  * Luppi and Coccia 2021 Luppi, E.; Coccia, E. Probing the molecular frame of uracil and thymine with high-harmonic generation spectroscopy. _Phys. Chem. Chem. Phys._ 2021, _23_ , 3729–3738
  * Morassut et al. 2022 Morassut, C.; Luppi, E.; Coccia, E. A TD-CIS study of high-harmonic generation of uracil cation fragments. _Chem. Phys._ 2022, _559_ , 111515
  * Luppi and Coccia 2023 Luppi, E.; Coccia, E. Role of Inner Molecular Orbitals in High-Harmonic Generation Spectra of Aligned Uracil. _J. Phys. Chem. A_ 2023, _127_ , 7335–7343
  * 76 Woźniak, A. P.; Lewenstein, M.; Moszyński, R. submitted to Phys. Rev. A
  * Ganeev et al. 2009 Ganeev, R. A.; Bom, L. B. E.; Abdul-Hadi, J.; Wong, M. C. H.; Brichta, J. P.; Bhardwaj, V. R.; Ozaki, T. Higher-Order Harmonic Generation from Fullerene by Means of the Plasma Harmonic Method. _Phys. Rev. Lett._ 2009, _102_ , 013903
  * Ganeev et al. 2009 Ganeev, R. A.; Elouga Bom, L. B.; Wong, M. C. H.; Brichta, J.-P.; Bhardwaj, V. R.; Redkin, P. V.; Ozaki, T. High-order harmonic generation from ${\text{C}}_{60}$-rich plasma. _Phys. Rev. A_ 2009, _80_ , 043808
  * Ganeev et al. 2009 Ganeev, R.; Singhal, H.; Naik, P.; Chakera, J.; Srivastava, A.; Dhami, T.; Joshi, M.; Gupta, P. Influence of C60 morphology on high-order harmonic generation enhancement in fullerene-containing plasma. _J. Appl. Phys._ 2009, _106_
  * Ganeev et al. 2013 Ganeev, R. A.; Hutchison, C.; Witting, T.; Frank, F.; Weber, S.; Okell, W. A.; Fiordilino, E.; Cricchio, D.; Persico, F.; Zaïr, A. et al. High-order harmonic generation in fullerenes using few-and multi-cycle pulses of different wavelengths. _J. Opt. Soc. Am. B_ 2013, _30_ , 7–12
  * Ganeev et al. 2016 Ganeev, R. A.; Suzuki, M.; Kuroda, H. High-order harmonic generation in Ag, Sn, fullerene, and graphene nanoparticle-contained plasmas using two-color mid-infrared pulses. _Eur. Phys. J. D_ 2016, _70_ , 1–8
  * Ganeev 2011 Ganeev, R. Fullerenes: the attractive medium for harmonic generation. _Laser Phys._ 2011, _21_ , 25–43
  * Ciappina et al. 2007 Ciappina, M. F.; Becker, A.; Jaro ń Becker, A. Multislit interference patterns in high-order harmonic generation in ${\mathrm{C}}_{60}$. _Phys. Rev. A_ 2007, _76_ , 063406
  * Ciappina et al. 2008 Ciappina, M. F.; Becker, A.; Jaro ń Becker, A. High-order harmonic generation in fullerenes with icosahedral symmetry. _Phys. Rev. A_ 2008, _78_ , 063405
  * Topcu et al. 2019 Topcu, T.; Bleda, E. A.; Altun, Z. Drastically enhanced high-order harmonic generation from endofullerenes. _Phys. Rev. A_ 2019, _100_ , 063421
  * Zhang 2005 Zhang, G. P. Optical High Harmonic Generation in ${\mathrm{C}}_{60}$. _Phys. Rev. Lett._ 2005, _95_ , 047401
  * Zhang and George 2006 Zhang, G. P.; George, T. F. Ellipticity dependence of optical harmonic generation in ${\mathrm{C}}_{60}$. _Phys. Rev. A_ 2006, _74_ , 023811
  * Avetissian et al. 2021 Avetissian, H. K.; Ghazaryan, A. G.; Mkrtchian, G. F. High harmonic generation in fullerene molecules. _Phys. Rev. B_ 2021, _104_ , 125436
  * Redkin and Ganeev 2010 Redkin, P. V.; Ganeev, R. A. Simulation of resonant high-order harmonic generation in a three-dimensional fullerenelike system by means of a multiconfigurational time-dependent Hartree-Fock approach. _Phys. Rev. A_ 2010, _81_ , 063825
  * Redkin et al. 2011 Redkin, P. V.; Danailov, M. B.; Ganeev, R. A. Endohedral fullerenes: A way to control resonant high-order harmonic generation. _Phys. Rev. A_ 2011, _84_ , 013407
  * Zhang and Bai 2020 Zhang, G. P.; Bai, Y. H. High-order harmonic generation in solid ${\mathrm{C}}_{60}$. _Phys. Rev. B_ 2020, _101_ , 081412
  * Castro et al. 2015 Castro, A.; Rubio, A.; Gross, E. K. U. Enhancing and controlling single-atom high-harmonic generation spectra: a time-dependent density-functional scheme. _Eur. Phys. J. B_ 2015, _88_ , 191
  * Artemyev et al. 2017 Artemyev, A. N.; Cederbaum, L. S.; Demekhin, P. V. Impact of two-electron dynamics and correlations on high-order-harmonic generation in He. _Phys. Rev. A_ 2017, _95_ , 033402
  * Sato et al. 2018 Sato, T.; Pathak, H.; Orimo, Y.; Ishikawa, K. L. Communication: Time-dependent optimized coupled-cluster method for multielectron dynamics. _J. Chem. Phys._ 2018, _148_ , 051101
  * Neufeld and Cohen 2020 Neufeld, O.; Cohen, O. Probing ultrafast electron correlations in high harmonic generation. _Phys. Rev. Res._ 2020, _2_ , 033037
  * Reiff et al. 2020 Reiff, R.; Joyce, T.; Jaroń-Becker, A.; Becker, A. Single-active electron calculations of high-order harmonic generation from valence shells in atoms for quantitative comparison with TDDFT calculations. _J. Phys. Commun._ 2020, _4_ , 065011
  * Nguyen and Bandrauk 2006 Nguyen, N. A.; Bandrauk, A. D. Electron correlation of one-dimensional ${\mathrm{H}}_{2}$ in intense laser fields: Time-dependent extended Hartree-Fock and time-dependent density-functional-theory approaches. _Phys. Rev. A_ 2006, _73_ , 032708
  * Saalfrank et al. 2020 Saalfrank, P.; Bedurke, F.; Heide, C.; Klamroth, T.; Klinkusch, S.; Krause, P.; Nest, M.; Tremblay, J. C. In _Chemical Physics and Quantum Chemistry_ ; Ruud, K., Brändas, E. J., Eds.; Advances in Quantum Chemistry; Academic Press, 2020; Vol. 81; pp 15–50
  * Stück et al. 2011 Stück, D.; Baker, T. A.; Zimmerman, P.; Kurlancheek, W.; Head-Gordon, M. On the nature of electron correlation in C60. _J. Chem. Phys._ 2011, _135_ , 11B608
  * Runge and Gross 1984 Runge, E.; Gross, E. K. U. Density-functional theory for time-dependent systems. _Phys. Rev. Lett._ 1984, _52_ , 997–1000
  * Tong and Chu 1998 Tong, X.-M.; Chu, S.-I. Time-dependent density-functional theory for strong-field multiphoton processes: Application to the study of the role of dynamical electron correlation in multiple high-order harmonic generation. _Phys. Rev. A_ 1998, _57_ , 452–461
  * Castro et al. 2004 Castro, A.; Marques, M. A. L.; Rubio, A. Propagators for the time-dependent Kohn–Sham equations. _J. Chem. Phys._ 2004, _121_ , 3425–3433
  * Isborn and Li 2008 Isborn, C. M.; Li, X. Modeling the doubly excited state with time-dependent Hartree–Fock and density functional theories. _J. Chem. Phys._ 2008, _129_
  * Habenicht et al. 2014 Habenicht, B. F.; Tani, N. P.; Provorse, M. R.; Isborn, C. M. Two-electron Rabi oscillations in real-time time-dependent density-functional theory. _J. Chem. Phys._ 2014, _141_
  * Castro et al. 2012 Castro, A.; Werschnik, J.; Gross, E. K. U. Controlling the Dynamics of Many-Electron Systems from First Principles: A Combination of Optimal Control and Time-Dependent Density-Functional Theory. _Phys. Rev. Lett._ 2012, _109_ , 153603
  * Kvaal 2012 Kvaal, S. Ab initio quantum dynamics using coupled-cluster. _J. Chem. Phys._ 2012, _136_ , 194109
  * Ridley and Zerner 1973 Ridley, J.; Zerner, M. An intermediate neglect of differential overlap technique for spectroscopy: pyrrole and the azines. _Theor. Chim. Acta_ 1973, _32_ , 111–134
  * Ghosh et al. 2017 Ghosh, S.; Andersen, A.; Gagliardi, L.; Cramer, C. J.; Govind, N. Modeling optical spectra of large organic systems using real-time propagation of semiempirical effective Hamiltonians. _J. Chem. Theory Comput._ 2017, _13_ , 4410–4420
  * Dreuw and Head-Gordon 2005 Dreuw, A.; Head-Gordon, M. Single-reference ab initio methods for the calculation of excited states of large molecules. _Chem. Rev._ 2005, _105_ , 4009–4037
  * Casida 1995 Casida, M. E. _Recent Advances In Density Functional Methods: (Part I)_ ; World Scientific, 1995; pp 155–192
  * McLachlan and Ball 1964 McLachlan, A. D.; Ball, M. A. Time-Dependent Hartree-Fock Theory for Molecules. _Rev. Mod. Phys._ 1964, _36_ , 844–855
  * Rowe 1968 Rowe, D. Equations-of-motion method and the extended shell model. _Rev. Mod. Phys._ 1968, _40_ , 153
  * Fabiano and Della Sala 2012 Fabiano, E.; Della Sala, F. Accuracy of basis-set extrapolation schemes for DFT-RPA correlation energies in molecular calculations. _Theor. Chem. Acc._ 2012, _131_ , 1278
  * Ziegler et al. 2014 Ziegler, T.; Krykunov, M.; Autschbach, J. Derivation of the RPA (random phase approximation) equation of ATDDFT (adiabatic time dependent density functional ground state response theory) from an excited state variational approach based on the ground state functional. _J. Chem. Theory Comput._ 2014, _10_ , 3980–3986
  * Heßelmann 2015 Heßelmann, A. Molecular excitation energies from time-dependent density functional theory employing random-phase approximation hessians with exact exchange. _J. Chem. Theory Comput._ 2015, _11_ , 1607–1620
  * Čížek and Paldus 1967 Čížek, J.; Paldus, J. Stability Conditions for the Solutions of the Hartree—Fock Equations for Atomic and Molecular Systems. Application to the Pi-Electron Model of Cyclic Polyenes. _J. Chem. Phys._ 1967, _47_ , 3976–3985
  * Thiel 2014 Thiel, W. Semiempirical quantum–chemical methods. _WIREs Comput. Mol. Sci._ 2014, _4_ , 145–157
  * Voityuk 2013 Voityuk, A. A. Intermediate neglect of differential overlap for spectroscopy. _WIREs Comput. Mol. Sci._ 2013, _3_ , 515–527
  * Ridley and Zerner 1976 Ridley, J. E.; Zerner, M. C. Triplet states via intermediate neglect of differential overlap: benzene, pyridine and the diazines. _Theor. Chim. Acta_ 1976, _42_ , 223–236
  * Bacon and Zerner 1979 Bacon, A. D.; Zerner, M. C. An intermediate neglect of differential overlap theory for transition metal complexes: Fe, Co and Cu chlorides. _Theor. Chim. Acta_ 1979, _53_ , 21–54
  * Zerner et al. 1980 Zerner, M. C.; Loew, G. H.; Kirchner, R. F.; Mueller-Westerhoff, U. T. An intermediate neglect of differential overlap technique for spectroscopy of transition-metal complexes. Ferrocene. _J. Am. Chem. Soc._ 1980, _102_ , 589–599
  * Mataga and Nishimoto 1957 Mataga, N.; Nishimoto, K. Electronic structure and spectra of nitrogen heterocycles. _Z. Phys. Chem._ 1957, _13_ , 140–157
  * Cronstrand et al. 2000 Cronstrand, P.; Christiansen, O.; Norman, P.; Ågren, H. Theoretical calculations of excited state absorption. _Phys. Chem. Chem. Phys._ 2000, _2_ , 5357–5363
  * Yeager et al. 1975 Yeager, D. L.; Nascimento, M. A. C.; McKoy, V. Some applications of excited-state-excited-state transition densities. _Phys. Rev. A_ 1975, _11_ , 1168–1174
  * Klinkusch et al. 2009 Klinkusch, S.; Saalfrank, P.; Klamroth, T. Laser-induced electron dynamics including photoionization: A heuristic model within time-dependent configuration interaction theory. _J. Chem. Phys._ 2009, _131_ , 114304
  * Coccia et al. 2017 Coccia, E.; Assaraf, R.; Luppi, E.; Toulouse, J. Ab initio lifetime correction to scattering states for time-dependent electronic-structure calculations with incomplete basis sets. _J. Chem. Phys._ 2017, _147_ , 014106
  * Frisch et al. 2016 Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Petersson, G. A.; Nakatsuji, H. et al. Gaussian˜16 Revision C.01. 2016; Gaussian Inc. Wallingford CT
  * Becke 1993 Becke, A. D. Density‐functional thermochemistry. III. The role of exact exchange. _J. Chem. Phys._ 1993, _98_ , 5648–5652
  * Yanai et al. 2004 Yanai, T.; Tew, D. P.; Handy, N. C. A new hybrid exchange–correlation functional using the Coulomb-attenuating method (CAM-B3LYP). _Chem. Phys. Lett._ 2004, _393_ , 51–57
  * Rostov et al. 2010 Rostov, I. V.; Amos, R. D.; Kobayashi, R.; Scalmani, G.; Frisch, M. J. Studies of the ground and excited-state surfaces of the retinal chromophore using CAM-B3LYP. _J. Phys. Chem. B_ 2010, _114_ , 5547–5555
  * Sarkar et al. 2021 Sarkar, R.; Boggio-Pasqua, M.; Loos, P.-F.; Jacquemin, D. Benchmarking TD-DFT and wave function methods for oscillator strengths and excited-state dipole moments. _J. Chem. Theory Comput._ 2021, _17_ , 1117–1132
  * Sun et al. 2020 Sun, Q.; Zhang, X.; Banerjee, S.; Bao, P.; Barbry, M.; Blunt, N. S.; Bogdanov, N. A.; Booth, G. H.; Chen, J.; Cui, Z.-H. et al. Recent developments in the PySCF program package. _J. Chem. Phys._ 2020, _153_ , 024109
  * Stewart 1990 Stewart, J. J. MOPAC: a semiempirical molecular orbital program. _J. Comput. Aided Mol. Des._ 1990, _4_ , 1–103
  * Gieseking 2021 Gieseking, R. L. A new release of MOPAC incorporating the INDO/S semiempirical model with CI excited states. _J. Comput. Chem._ 2021, _42_ , 365–378
  * 135 Moussa, J. E.; Stewart, J. J. P. MOPAC. https://http://openmopac.net
  * Chi 1970 Chi, B. E. The eigenvalue problem for collective motion in the random phase approximation. _Nucl. Phys. A_ 1970, _146_ , 449–456
  * Morassut et al. 2023 Morassut, C.; Coccia, E.; Luppi, E. Quantitative performance analysis and comparison of optimal-continuum Gaussian basis sets for high-harmonic generation spectra. _J. Chem. Phys._ 2023, _159_
  * Berkelbach 2018 Berkelbach, T. C. Communication: Random-phase approximation excitation energies from approximate equation-of-motion coupled-cluster doubles. _J. Chem. Phys._ 2018, _149_
  * Gensterblum et al. 1991 Gensterblum, G.; Pireaux, J. J.; Thiry, P. A.; Caudano, R.; Vigneron, J. P.; Lambin, P.; Lucas, A. A.; Krätschmer, W. High-resolution electron-energy-loss spectroscopy of thin films of ${\mathrm{C}}_{60}$ on Si(100). _Phys. Rev. Lett._ 1991, _67_ , 2171–2174
  * Hertel et al. 1992 Hertel, I. V.; Steger, H.; de Vries, J.; Weisser, B.; Menzel, C.; Kamke, B.; Kamke, W. Giant plasmon excitation in free ${\mathrm{C}}_{60}$ and ${\mathrm{C}}_{70}$ molecules studied by photoionization. _Phys. Rev. Lett._ 1992, _68_ , 784–787
  * Yoo et al. 1992 Yoo, R.; Ruscic, B.; Berkowitz, J. Vacuum ultraviolet photoionization mass spectrometric study of C6. _J. Chem. Phys._ 1992, _96_ , 911–918
  * Ju et al. 1993 Ju, N.; Bulgac, A.; Keller, J. W. Excitation of collective plasmon states in fullerenes. _Phys. Rev. B_ 1993, _48_ , 9071–9079
  * Scully et al. 2005 Scully, S. W. J.; Emmons, E. D.; Gharaibeh, M. F.; Phaneuf, R. A.; Kilcoyne, A. L. D.; Schlachter, A. S.; Schippers, S.; Müller, A.; Chakraborty, H. S.; Madjet, M. E. et al. Photoexcitation of a Volume Plasmon in ${\mathrm{C}}_{60}$ Ions. _Phys. Rev. Lett._ 2005, _94_ , 065503