{'fraction_non_alphanumeric': 0.03770543950800551, 'fraction_numerical': 0.00962782313646485, 'mean_word_length': 4.73813579946459, 'pattern_counts': {'":': 1, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Understanding real-time time-dependent density-functional theory simulations of ultrafast laser-induced dynamics in organic molecules

Real-time time-dependent density functional theory, in conjunction with the Ehrenfest molecular dynamics scheme, is becoming a popular methodology to investigate ultrafast phenomena on the nanoscale. Thanks to recent developments, it is also possible to explicitly include in the simulations a time-dependent laser pulse, thereby accessing the transient excitation regime. However, the complexity entailed in these calculations calls for in-depth analysis of the accessible yet approximate (either"dressed"or"bare") quantities in order to evaluate their ability to provide us with a realistic picture of the simulated processes. In this work, we analyze the ultrafast dynamics of three small molecules (ethylene, benzene, and thiophene) excited by a resonant laser pulse in the framework of the adiabatic local-density approximation. The electronic response to the laser perturbation in terms of induced dipole moment and excited-state population is compared to the results given by an exactly solvable two-level model. In this way, we can interpret the charge-carrier dynamics in terms of simple estimators, such as the number of excited electrons. From the computed transient absorption spectra we unravel the appearance of nonlinear effects such as excited-state absorption and vibronic coupling. In this way, we observe that the laser excitation affects the vibrational spectrum by enhancing the anharmonicities therein while the coherent vibrational motion contributes to stabilize the electronic excitation already within a few tens of femtoseconds.

Real-time time-dependent density functional theory, in conjunction with the Ehrenfest molecular dynamics scheme, is becoming a popular methodology to investigate ultrafast phenomena on the nanoscale. Thanks to recent developments, it is also possible to explicitly include in the simulations a time-dependent laser pulse, thereby accessing the transient excitation regime. However, the complexity entailed in these calculations calls for in-depth analysis of the accessible yet approximate (either "dressed" or "bare") quantities in order to evaluate their ability to provide us with a realistic picture of the simulated processes. In this work, we analyze the ultrafast dynamics of three small molecules (ethylene, benzene, and thiophene) excited by a resonant laser pulse in the framework of the adiabatic local-density approximation. The electronic response to the laser perturbation in terms of induced dipole moment and excited-state population is compared to the results given by an exactly solvable two-level model. In this way, we can interpret the charge-carrier dynamics in terms of simple estimators, such as the number of excited electrons. From the computed transient absorption spectra we unravel the appearance of nonlinear effects such as excited-state absorption and vibronic coupling. In this way, we observe that the laser excitation affects the vibrational spectrum by enhancing the anharmonicities therein while the coherent vibrational motion contributes to stabilize the electronic excitation already within a few tens of femtoseconds.

I. INTRODUCTION
Real-time time-dependent density-functional theory (RT-TDDFT) is a mature methodology to study atomistically the ultrafast dynamics of systems up to thousands of atoms. 1 Being intrinsically non-perturbative, this approach can access on the same footing both the linear and the nonlinear regime of excitations. 2,3 Recent developments of RT-TDDFT have added the possibility to explicitly include a time-dependent electric field, [4][5][6][7][8][9][10][11][12] paving the way for RT-TDDFT to simulate pump-probe experiments and thus to investigate quantitatively the dynamics of charge carriers on their natural femtosecond timescale. Merging of RT-TDDFT with the Ehrenfest molecular dynamics scheme [13][14][15][16] additionally allows us to describe the coupled electron-nuclear dynamics in the sub-picosecond time window, adopting a mixed quantum-classical treatment.
With the capabilities illustrated above, it is possible for RT-TDDFT to provide us with insights about excitation processes in the femtosecond regime which are now revealed by the most recent laser technologies. [17][18][19][20][21][22][23] The possibility to access the (transient) response of the system before, during, and after laser illumination, as well as the flexibility to decouple the effects of electronic and nuclear response, make the RT-TDDFT simulations complementary to laboratory observations. The interpretation of RT-TDDFT results, however, can be far from trivial. First, the single-particle orbitals used to reproduce the dynamics of the particle density do not have any stringent physical meaning unless the effects of the electronelectron interactions are somewhat small. Second, applications must unavoidably make use of approximations. Especially in dealing with large systems, restrictions to functional forms that can only capture local correlations, either in space, in time, or both, are often invoked to make calculations doable. Consequently, additional analysis is needed to eval-uate the validity of the corresponding results. Such a task is not simplified by the fact that accurate numerical benchmarks may not be available and that reference experiments often deal with novel phenomena that are non-trivial to interpret.
In this work, we show that some additional and yet minimalist modeling can help understanding and thus interpreting RT-TDDFT simulations. In particular, we investigate the laser-induced dynamics of three small molecules (ethylene, benzene, and thiophene) which can be viewed as the building blocks of larger organic complexes of experimental and technological interest. For this purpose, we employ RT-TDDFT in the adiabatic local-density approximation (ALDA), [24][25][26][27] which is widely used especially in calculations of large systems due its reduced computational costs compared to more advanced approximations. After the analysis of the linear absorption properties of the three molecules, we investigate their response to a resonant femtosecond (fs) laser pulse, monitoring relevant quantities such as the induced dipole moment and photo-excited charge-carrier population. The comparison with the results obtained in the same conditions from an exactly solvable two-level model offers a useful diagnostic tool to interpret the laser-induced electron dynamics from RT-TDDFT, and to highlight the physical meaning of the corresponding outcomes. The analysis of the transient absorption spectra allows to illustrate the capability of RT-ALDA to reproduce entailed nonlinear processes such as excited-state absorption. Finally, the effects of the coupled electron-nuclear dynamics are discussed in details to disclose the role of coherent nuclear motion and vibronic coupling in stabilizing the excitation in the fs time window. This paper is organized as follows. In Section II we review the methodology and report the computational details. Section III contains the body of the results, starting from the linear absorption properties of the systems (Section III A) and proceeding with the outcomes of the laser-induced dynamics of the electrons (Section III B), the transient absorption spectra (Section III C), excited-state absorption (Section III D), and vibronic coupling (Section III E). The conclusions are reported in Section IV.

A. Theoretical background
The results presented in this work are obtained in the framework of time-dependent DFT, 28 adopting its real-time implementation to access the electron dynamics. Within RT-TDDFT, the time-dependent (TD) electron density for a spinunpolarized system is calculated through the propagation of the Kohn-Sham orbitals φ i (r,t) according to the TD Kohn-Sham equation where (3) is the effective potential. The first two terms describe the interaction of the electrons with the nuclei and with external fields, respectively, while the last two account for the electronelectron interactions. The exchange-correlation (XC) potential, v xc , must be approximated, as its exact form is not known. This term is nonlocal in space and time, 29 meaning that it depends not only on the instantaneous density, but also on the density at earlier times. 10 While efforts have been directed towards taking into account memory effects, 29 in most calculations, the adiabatic approximation is still assumed, i.e., the XC potential is approximated using a ground-state (GS) XC potential with the TD density. Here, we adopt the ALDA, where v xc is assumed to be local in space and time. The locality in space, which is inherited from the local-density approximation for v xc in the GS calculation, is known to impair the ability of TDDFT to correctly predict the energy of charge-transfer excitations [30][31][32][33][34][35][36] and excitonic features. [37][38][39][40] The locality in time stems from the adiabatic approximation, i.e., the approximation of the TD XC functional as a GS XC functional of the instantaneous density. This assumption leads to a neglect of memory effects, which are thought to become crucial especially in real-time propagation schemes, where the system can be driven far away from its GS. 29,[41][42][43] For these reasons, the accuracy of the ALDA results is often questioned (see for example Refs. 44 and 45). Nonetheless, the ALDA remains a popular approach especially thanks to its numerical simplicity, which enables the treatment of systems otherwise out of reach by means of more accurate approaches.
The coupling of the electronic system to an external TD electric field E(t) is described semi-classically in the dipole approximation and within the length gauge, such that the second term in Eq. (3) is given by v ext (r,t) = r · E(t). Absorption spectra are computed according to the scheme proposed in Ref. 2, i.e., the system is probed by an instantaneous, broadband electric field, commonly referred to as a "kick": where n defines the polarization direction, κ the kick strength, and δ (t) is the Dirac delta function. Perturbing the system with a kick is equivalent to multiplying each occupied orbital with a phase factor: 2 To calculate the linear absorption spectrum, we start from the GS wave functions, φ i (r, 0 − ) = φ i (r). After the instantaneous perturbation given by the kick, the system is propagated for a certain amount of time. The total duration of the propagation determines the achieved spectral resolution. The resulting induced dipole moment, is a superposition of free oscillations with frequencies and amplitudes corresponding to the energies and oscillator strengths of the excitations of the system, respectively. Here, ρ(r,t), defined in Eq. (1), and ρ g (r) are the TD and GS electron densities, respectively. The absorption spectrum is obtained by taking the imaginary part of the Fourier transform of d(t) as obtained from Eq. (6). If κ is sufficiently small, only the linear response is obtained, while stronger perturbations excite in higher order. 3 In order to calculate the response of the material to an ultrafast laser pulse, we consider electric fields of the form Here, E 0 is the peak amplitude, f (t) an envelope function, and ω p the carrier frequency, which can be tuned to be in resonance with specific excitation energies of the system. It is possible to combine this realistically shaped laser pulse and the instantaneous kick mentioned above in order to simulate pump-probe experiments. 4 In this case, the electric field can be written as where the first term, referred to as the "pump", prepares the system in a non-stationary state, and the second one, the "probe", interrogates it at a later time t 0 . In this way, it is possible to gain information about excited and transient states of the system. The dipole moment induced by the two fields includes interference effects depending on their phase relation. These effects are not measured in typical pump-probe experiments since they vanish when taking the average over the response of the system to multiple laser shots. 46 To obtain a transition dipole moment that is insensitive to the aforementioned phase difference, we adopt the so-called "phase-cycling approach" 47 which consists in performing four simulations applying the field in Eq. (8), adding each time a constant phase of π/2 to the pump field. In this way, the applied laser field becomes: where j runs from 1 to 4. The four resulting total dipole moments are subsequently averaged and the pump-only contributions are subtracted. A sine transformation finally yields the absorption spectrum. In order to calculate the transient absorption spectrum (TAS), this procedure is repeated for different probe delay times t 0 , corresponding to instants before, during, and after the application of the pump. The non-perturbative approach provided by RT-TDDFT does not offer any information about the single-particle transitions that contribute to the absorption in linear regime. For this reason, this method is often complemented by the solution of the TD Kohn-Sham equations in the linear-response scheme first proposed by Mark Casida. 48 Within this framework, the problem is linearized and reformulated in terms of the matrix equation 49 Here, X and Y are vectors including the excitation and deexcitation coefficients X im and Y im , respectively, which correspond to transitions from the occupied orbitals φ i (r) to the unoccupied ones φ m (r), and vice versa. The eigenvalue ω represents the excitation energy, and, within the adiabatic approximation, the matrices A and B are given by A im, jn = δ i, j δ m,n (ε m − ε i ) + im| f Hxc |n j (11) and B im, jn = im| f Hxc | jn , respectively. In Eq. (11), ε k are the KS eigenvalues, and the Hartree-XC kernel appearing in Eqs. (11) and (12) reads: where v xc [n](r) is the GS XC potential and ρ g (r) the GS electron density [already introduced in Eq. (6)]. Eq. (10) yields the excitation energies ω e and the corresponding transition weights X (e) and Y (e) of the excited many-body states |e . Given the known limitations of ALDA in describing excitonic effects, 37-40 we employ many-body perturbation theory (MBPT) as a benchmark for the linear absorption spectra. To this end, we solve the Bethe-Salpeter equation (BSE) for a system of quasiparticles with energies calculated within the G 0 W 0 approximation. 50,51 Formally, the BSE is the equation of motion for the two-particle polarizability, L = L 0 +L 0 KL, 52 where the interaction kernel K includes the repulsive exchange term given by the bare Coulomb potential and the statically screened electron-hole Coulomb attraction. For further details about this formalism, see e.g., Refs. 40,51,53,54.
One of the advantages of RT-TDDFT is the possibility to interface it with molecular dynamics, e.g. adopting the Ehrenfest scheme. [13][14][15] The (classical) equation of motion for the Jth nucleus out of M nuclei in total is herein given by 14 where M K and R K are the mass and position of the Kth nucleus, respectively. The term v(r, R 1 , . . . , R M ) includes the electrostatic energy among the nuclei and between the nuclei and an electron at position r, as well as the energy related to the coupling of the nuclei to an external field. The Ehrenfest formalism is a mean-field approach, with the forces on the classical trajectory being calculated by averaging over the quantum-mechanical degrees of freedom. 16

B. Computational details
All RT-TDDFT calculations in this work are performed with the OCTOPUS code. [55][56][57] Core electrons are approximated by norm-conserving Troullier-Martins pseudopotentials 58 , calculated within the LDA in the Perdew-Zunger parameterization. 59 The pseudo-orbitals (hereafter referred to as orbitals for brevity) are represented on a realspace grid with a spacing of 0.3 bohr inside a simulation box built by interlocking spheres with a radius of 12 bohr centered at each atom. The orbitals are propagated using the approximated enforced time-reversal symmetry scheme 60 with a time step of 1.2 attoseconds. To obtain the absorption spectra, each system is propagated for 10 fs after the application of the kick, resulting in a full width at half maximum of the peaks of 435 meV. For the laser fields, we use a peak amplitude E 0 corresponding to an intensity of 5 × 10 11 W/cm 2 (unless stated otherwise), and a Gaussian envelope function f (t) with standard deviation of 2 fs. The carrier frequency ω p is chosen in resonance with specific excitations in each system, identified in their RT-ALDA spectra. TAS are calculated by probing the system at intervals of 1 fs. In the RT-ALDA calculations including coupled electron-nuclear dynamics, we start from a geometry relaxed using the FIRE algorithm 61 minimizing the forces until they are lower than 10 −3 Ha/bohr. No initial velocities are assigned to the nuclei unless stated otherwise.
Linear-response spectra are computed with the MOLGW code 53 , using Gaussian-type cc-pVQZ basis sets. 62 To calculate the MBPT spectra, the resolution-of-identity approximation 63 is employed. The range-separated hybrid functional CAM-B3LYP 64 is used in the underlying DFT calculation as an improved starting point for the single-shot G 0 W 0 approximation . 51,65 FIG. 1. Absorption spectra of a) ethylene, b) thiophene, and c) benzene calculated from MBPT and RT-ALDA. The excitations targeted by the resonant pulse are highlighted by green arrows. Insets: isosurfaces of the difference density (∆ρ (e) (r) ranges from 1 × 10 −3 Å −3 to 5 × 10 −3 Å −3 depending on the molecule) related to the these excitations on top of ball-and-stick representations of the molecules (black, white, and yellow spheres stand for carbon, hydrogen, and sulfur atoms, respectively). Blue and red colors in the isosurfaces denote electron depletion and accumulation, respectively.

A. Linear absorption spectra
We start our analysis by examining the linear absorption spectra of ethylene, thiophene, and benzene, depicted in Fig. 1, as obtained from RT-ALDA and MBPT. Focusing on the pronounced maximum at around 7.6 eV in the spectrum of ethylene [ Fig. 1a)], we find a good match between the RT-ALDA and the MBPT results. This peak is given by a π → π * excitation and has transition dipole moment (TDM) in the x direction [see coordinate systems in the inset of Fig. 1a)]. According to the linear-response calculations (see Supporting Information, Table S1 and Fig. S2), this excitation corresponds mainly to a transition from the highest occupied molecular orbital (HOMO) to the lowest unoccupied molecular orbital (LUMO). To visualize the character of this excitation, we plot the corresponding difference density, 66 calculated as where the transition weights X (e) im are obtained from the solution of the Casida's equation, Eq. (10). On the left-handside, ∆ρ (e) (r) represents the difference between the excitedstate and the ground-state electron densities. In the case of ethylene, we see that the electron density accumulates on the opposite ends of the molecule along the direction of the carbon double bond [inset of Fig. 1a)]. This characteristic is particularly important in the analysis of the coupled electronnuclear dynamics, which is discussed in Section III E below. The weak peak at lowest energy in the RT-ALDA spectrum [ Fig. 1a)] is given by the transition from the HOMO to the LUMO+1, and has a TDM oriented in the z direction. Due to the small overlap between the involved orbitals (see Fig.  S1), the excitation energy is underestimated by ALDA, which notoriously fails to correctly predict the energy of chargetransfer excitations. 31,32,[35][36][37] Indeed, in the MBPT spectrum, the excitation appears above the absorption onset. The incorrect order of the first two excitations given by the ALDA does, however, not affect the following analysis on the laser-induced electronic and nuclear dynamics, as we can selectively tune the laser frequency to be in resonance with a specific excitation. Spurious resonance effects between these two excitations, which could be anticipated due to their energetic proximity, can be ruled out due to the perpendicular orientation of the corresponding TDMs.
In the case of thiophene [ Fig. 1b)], the agreement between the MBPT and TDDFT spectra is very good, especially for the first peaks. The first maximum comprises two bright excitations energetically very close to each other. The lower energy one has a TDM oriented in the y direction, while the second one is polarized along x. The analysis of the linear-response calculation (see Table S2) shows that the first excitation corresponds mainly to the transition from the HOMO-1 to the LUMO, but has a significant contribution from the HOMO to the LUMO+3 as well. The x-polarized excitation is given mostly by the HOMO-LUMO transition.
Considering now benzene [see Fig. 1c)] we see that the RT-ALDA and MBPT spectra match qualitatively. Both methods predict an intense peak formed by a doubly degenerate excitation with contributions from transitions between the doubly degenerate HOMO and LUMO (see Table S3). This excitation is again a π → π * transition with the electron density accumulating on the outer part of the aromatic ring, as shown in the inset of Fig. 1c).

B. Laser-induced dynamics
With the knowledge of the linear absorption spectra, in the next step we investigate the electron dynamics induced in the molecules by an ultrafast laser pulse in resonance with the ex- citations highlighted in Fig. 1 and discussed above. In particular, we analyze the induced dipole moment and the excitedstate population. To interpret our findings, we use as a point of comparison an exactly solvable problem, consisting of a two-level system perturbed to the same laser pulse. While this model is by itself subject to the crude approximation of twolevel excitations, its solution is exact and its physical meaning transparent. The Hamiltonian of the two-level system is given bŷ where |g and |e are the GS and excited state (ES), respectively; the TDM d ge is taken from the Casida calculation and the excitation energy (E e − E g ) stems from the RT-ALDA spectrum. We neglect the terms related to the permanent dipole of the molecules in the GS (d g ) and in the ES (d e ), which would contribute the additional terms −d g · E(t)|g g| and −d e · E(t)|e e|, since they do not have a significant impact upon the time evolution. The coefficients c g (t) and c e (t) in the TD wavefunction, are calculated numerically using the Euler method. These coefficients remain constant once the field is turned off. To relate this model to the results provided by RT-ALDA, we analyze the two-level system in terms of the electron density, which is obtained from the many-body wavefunction Ψ(r 1 , . . . , r N ,t) = r 1 . . . r N |Ψ(t) as Inserting the wavefunction in the form of Eq. (17) into Eq. (18), we find where the the transition density is defined following Ref. 66: The difference density ∆ρ (e) (r) = ρ e (r) − ρ g (r) corresponds to Eq. (15), if the ES is represented as a superposition of Slater determinants with coefficients obtained from the linear-response calculation. From the analytic solution of the two-level problem, it is known that the excited-state population |c e (t)| 2 changes on a timescale determined by the intensity of the field, whereas the coherence c * g (t)c e (t)e −i(E e −E g )t oscillates with the frequency of the resonant field. Hence the induced density is the sum of two terms: the first one is related to ∆ρ (e) (r) and changing at a rate determined by the field strength, while the second one gives rise to the coherent electronic motion of a superposition state, with a spatial profile given by ρ The induced dipole moment obtained from the two-level model [Eq. (21)] can be compared to that from the RT-ALDA calculation [Eq. (6)].
In the case of ethylene, the results obtained with the two methods agree well [see Fig. 2a)]. In both cases, the dipole moment builds up during laser irradiation and remains constant afterwards. A careful inspection of Fig. 6a) shows a frequency shift which increases over time and becomes clearly visible towards the end of the simulation. This behavior likely due to a shift of the absorption peak caused by pumping the system (see also Sec. III C), which is an unphysical artifact of the adiabatic approximation. 67 In the case of thiophene, the dipole moment computed from RT-ALDA has an additional contribution during the active time window of the laser [ Fig.  2b)]. It represents a non-dissipative, temporary change of the electronic density and is related to the real part of the polarizability, which in the case of thiophene is relatively large in comparison to its imaginary part. This effect is not included in the two-level model. After the field is switched off, the results of the model and of RT-ALDA are in excellent agreement. In the case of benzene, the envelope of the dipole moment increases once the field is applied [ Fig. 2c)]. In the two-level model, the envelope stays constant afterwards, whereas in the RT-ALDA result a slight oscillatory modulation of the amplitude is observed. These quantum beats suggest that the evolution of the dipole moment for this system can involve a threestate superposition including the GS and two excited states, which must be close in energy since the beat frequency is low. Next, we compare the ES population obtained from the model system, with the number of excited electrons, N ex (t) computed from RT-ALDA in an independent-particle manner, i.e. by projecting all occupied time-dependent orbitals onto the static unoc-cupied ones, 9 where N is the total number of electrons and the normalization of the TD orbitals is used: In case only one single excitation is affected by the laser, P e (t) and N ex (t) should be comparable, as the transition from |g to |e is a single-electron excitation. The curves in Fig. 3 obtained from RT-ALDA exhibit rapid oscillations, which are due to the sloshing of the electron density, but generally agree well with the ones obtained from the model. For all three molecules, both RT-ALDA and the model feature an increasing number of excited electrons as the laser is applied, and remain nearly constant afterwards. An excellent agreement between the results obtained with the two methods is found for thiophene [ Fig. 3b)]. On the other hand, differences appear in the case of benzene [see Fig. 3c)]: N ex (t) (RT-ALDA, black) and P e (t) (two-level model, dashed blue) differ after the pulse, with N ex (t) decreasing slightly, while P e (t) remains constant. This behavior of N ex (t) mirrors the decrease observed in the envelope of the induced dipole moment [ Fig. 2c)]. Generally, we observe that deviations become significant as population inversion is reached, i.e., as 50% excited-state population is exceeded. Previous studies indicate that in this scenario memory effects play a crucial role and no adequate description within the adiabatic approximation is possible. 42,43,68 We observe that the number of excited electrons deduced directly from the bare KS quantities [Eq. (25)] compares well with the populations given by the interacting two-level system [Eq. (22)]. This is a clear indication that the effects of the electronelectron interaction are rather small in our cases and, thus, even bare KS quantities, such as N ex , can directly provide us with informative results for physical interpretations.

C. Transient absorption spectra
We next investigate transient absorption spectra, which illustrate the temporal change of the absorption characteristics of a system during and after laser irradiation. For ethylene and benzene, we use the same pulses as before, while for thiophene, we increase the intensity by a factor of 10 to reach a similar amount of excited-state population as in the other molecules.
When only the electronic part of the dynamics is considered [see Fig. 4a-c)], changes in the TAS mainly occur in the time window in which the field is acting. In ethylene the resonance energy is 7.6 eV, in thiophene 5.6 eV, and in benzene 6.9 eV, also marked in Fig. 1. A common characteristic in the spectra of the three molecules is the red-shift of the pumped excitation probed during the pulse. This characteristic persists also after the field is switched off and is likely an artifact of the adiabatic approximation 67 (see Sec. III B). Another general feature is the emergence of new peaks in the lower-energy region of the spectra during laser irradiation. The energies of these pre-peaks remain constant after the field is switched off, as expected. 69,70 On the other hand, their absorption strength is either unchanged over time, as in the case of ethylene [ Fig.  4a)], or has an oscillatory time dependence, as in the other two molecules [ Fig. 4b-c)]. These oscillations are related to the quantum beats resulting from the system being in a coherent superposition of electronic states including more than one excitation. 69 When taking into account the nuclear motion [ Fig. 4d-f)], the time evolution is no longer trivial after the end of the pulse. The absorption energies vary in time also beyond the duration of the laser illumination. In all cases, the pumped peak exhibits a more pronounced red-shift than the one observed when the nuclei are clamped. This behavior can be explained in terms of the Stokes shift, namely to the red-shift due to the change of the geometry upon electronic excitation. Unfortunately, we are unable to quantify this effect, as the corresponding spectral signatures are superimposed with the red-shift of the pumped peak discussed above for the dynamics with fixed ions.
The behavior of the emerging pre-peaks varies depending on the system. In the TAS of ethylene, the distinct peak at 2 eV in Fig. 4d) exhibits an oscillatory blue-shift, which mirrors the red-shift of the pumped peak. In thiophene, the first prepeak at about 1.5 eV [ Fig. 4e)] is broadened in energy around 30 fs before narrowing and red-shifting. The second peak at about 3 eV remains instead largely unchanged, and the third one at approximately 4.5 eV exhibits a blue-shift similar to that observed in the TAS of ethylene. Finally, in benzene [ Fig.  4f)] the first pre-peak around 1 eV disappears completely after 30 fs, whereas the next three maxima at 2 eV, 3 eV, and 4 eV become more defined and gain oscillator strength with time in comparison to their dynamics with fixed nuclei.

D. Excited-state absorption
To better understand the nature of the pre-peaks appearing in the TAS (Fig. 4), we analyze in more details the feature appearing at around 2.3 eV in the spectrum of ethylene [ Fig. 4a)]. To do so, we interrogate the system by applying a monochromatic probe pulse in resonance with the excitation energy of the pre-peak (2.3 eV), after the excitation with the first resonant pump at 7.6 eV. The total field in this case is a pulse train, depicted in the upper panel of Fig. 5. To understand the resulting electron dynamics, we monitor the populations of the individual orbitals contributing to the total number  25) and smoothed out, upon the application of the two laser pulses sketched in the top panel. The first one is tuned to the first π → π * transition, the second one has a frequency corresponding to the energy of the prepeak appearing in the TAS (Fig. 4). The solid curve in orange is obtained applying only the first pulse depicted in the uppermost panel, while the dotted-dashed one in black is computed when the two pulses are applied subsequently, as illustrated in the inset. of excited electrons: Upon the application of the pumping laser, we find mainly the LUMO gaining population, but also a higher unoccupied state participating. As these orbitals contribute to one single excitation, we add them up and consider their sum as the occupation of S 1 (Fig. 5). The increase of the population of S 1 is accompanied by a symmetric decrease of the population of S 0 , which is calculated from the population of the static HOMO only, since the TD HOMO is the only active orbital. When the probe pulse is applied, the populations of the orbitals contributing to S 1 decrease, while three higher-energy orbitals gain population. Assuming that these three singleparticle states are associated with one single ES, we add them and identify their sum as the population of S 2 (details in the Supporting Information). Notably, the population of S 0 is not affected by the probe pulse. We thus conclude that the prepeak in the TAS is related to excited-state absorption from S 1 to S 2 about 2.3 eV higher in energy, as previously speculated in Ref. 4.

E. Vibronic couplings
After clarifying the nature of the peaks in the TAS, we now examine their dynamical behavior when the vibronic coupling is taken into account. In this analysis, we focus again on ethylene as in this molecule, formed by a carbon double bond with four terminating H atoms, it is straightforward to make a connection between the temporal evolution of the absorption peaks, the vibrational properties, and the excited-state geometry. As we are also interested in a possible change of the vibrational modes upon excitation, in this study we assign to the nuclei initial velocities, taken randomly from a Maxwell-Boltzmann distribution with a temperature of 300 K. We then propagate the system freely (i.e., without laser perturbation) for about 50 fs, to enable a uniform distribution of the kinetic energy among the nuclear degrees of freedom. Finally, we extract coordinates and velocities and use them as a starting point for two simulations: The first one is another free propagation, corresponding to the nuclear motion on the GS potential energy surface (PES); in the second calculation, the molecule is excited with a laser pulse, targeting the HOMO-LUMO transition at 7.6 eV. In the latter case, the dynamics take place on both the GS PES and the ES PES. In both cases, we monitor the time evolution of the C=C bond length, plotted in Fig. 6a).
From the inspection of Fig. 6a), we immediately notice that the electronically excited molecule displays a general increase of the C=C bond length. This behavior can be understood in more detail from the induced density [Eq. (19)] and from the difference density ∆ρ (e) (r) of the pumped excitation computed in the linear-response framework [Eq. (15)]. According to Eq. (19), the induced density δ ρ(r,t) is the sum of two terms: The first one is related to ∆ρ (e) (r) and changes relatively slowly (∼10 fs for the considered laser intensity), and only while the field is active. The density redistribution related to ∆ρ (e) (r) causes static forces on the positively charged nuclei, driving them away from their equilibrium positions. For the specific case of ethylene, the electron depletion along the C=C bond creates a positive charge in the center of the molecule that repels the positively charged nuclei, while the electron accumulation on the sides causes attractive forces [see inset of Fig. 1a)]. As a result, the carbon atoms are driven apart and the C=C bond is elongated. The second contribution to δ ρ(r,t) is related to the transition density ρ (e) t (r) [Eq. (20)] of the electronic excitation and corresponds to the coherent electronic motion. The resulting forces on the nuclei, however, have a strong oscillatory time dependence (< 1 fs, determined by the electronic excitation frequency), and average to zero. As such, they cannot trigger the nuclear motion effectively. Hence, the difference density ∆ρ (e) (r) provides us with more information about the vibronic coupling than the transition density ρ fs and 50 fs, we find a maximum formed by a double peak. This feature suggests that, upon laser illumination the motion of the C=C bond becomes a superposition of two oscillations with different frequencies. To investigate this matter further, we perform a Fourier transform of both curves shown in Fig. 6a) after extending the propagation up to 600 fs, in order to obtain a spectrum with sufficient frequency resolution [see Fig. 6b)]. The free propagation yields a pronounced peak close to the frequency of the C=C stretching mode (1650 cm −1 ), in agreement with the harmonic behavior observed in time domain [Fig. 6a)]. Another local maximum is recognizable in the vicinity of the C=C scissoring mode around 1350 cm −1 . When the molecule is electronically excited, this peak undergoes a red-shift of about 150 cm −1 , which we attribute to the anharmonicity of the oscillation. The most striking feature in the orange curve in Fig. 6b) is however the dramatic increase in strength of the peak at about 1350 cm −1 as a consequence of vibronic coupling. The frequency of this peak re-mains unchanged regardless of the laser perturbation, suggesting low anharmonicity in this motion, likely due to the small amplitude of the oscillation prior to the electronic excitation. The presence of two maxima with similar oscillator strength in the spectrum of the laser-excited molecule [ Fig. 6b] clarifies the double peak in Fig. 6a as a superposition of two oscillations in time domain. The displayed behavior of the vibrational motion in ethene depends on the specific initial conditions, and the same electric field will cause a different response when the molecule is in a different nuclear configuration. For a more general analysis, a statistical approach is required, performing the same calculations for a number of different starting conditions. The present results correspond to a particular scenario, but show potential effects that a laser excitation can exert with respect to the vibrational properties. Through the vibronic coupling, the electronic excitation energy is partially transferred to the nuclei, increasing the energy in the normal modes. Meanwhile, the modes can exhibit a frequency shift, which is likely related to the anharmonicity of the potential energy associated with the normal modes.
To conclude this analysis of the electron-nuclear dynamics of ethylene, we inspect the feedback of the nuclear motion to the electronic excitation. To do so, in Fig. 7 we plot the time-dependent induced density computed for fixed and moving ions. We notice that, immediately after the excitation, the inclusion of the nuclear motion generates only a slight difference in the induced density compared to the purely electronic contribution [ Fig. 7a) and b)]: The induced electronic density, integrated over the plane perpendicular to the bond axis [see coordinate system in Fig. 1a)], is very similar for fixed and moving nuclei. Towards the end of the pulse, the electrostatic forces on the nuclei due to the laser-induced electronic redistribution become noticeable. When the nuclei are allowed to move, the carbon atoms are driven apart and the bond length increases [inset of Fig.  7b)]. After approximately 10 fs, the bond length reaches a maximum and δ ρ(x,t) deviates significantly from the case in which the nuclei are fixed [ Fig. 7c) and d)]: The electron depletion between the carbon atoms and the accumulation on the outside are enhanced considerably when the nuclei are allowed to move [ Fig. 7d)]. When the nuclei are clamped, δ ρ(r,t) looks similar at t 10 fs and t 20 fs [Fig. 7a) and c)], since the electronic system propagates freely after the excitation. We compare δ ρ(x,t) to the plane-integrated difference density of the excitation, ∆ρ (e) (x) = dydz ∆ρ (e) (r) = dydz [ρ e (r) − ρ g (r)] , (27) where ∆ρ (e) (r) is calculated in linear response from Eq. (15). This is the value that δ ρ(r,t) would assume if the laser evoked 100% ES population. As this scenario does no longer represent a superposition state, the density is stationary. When the nuclear motion is included, δ ρ(x,t) shows a clear bias toward ∆ρ (e) (x) at 20 fs [ Fig. 7d)]. This is evident when comparing The nuclear motion strongly enhances the charge redistribution of the excitation. It should be emphasized that at t=20.1 fs the C=C bond length exhibits a maximum. In the subsequent nuclear oscillations, the electron density continuously adapts to the vibrational phase, and in the minima of the bond length the enhancement is substantially weaker. However, since the mean bond length increases [ Fig. 6a)], there is a net enhancement of electron localization that remains after averaging over a vibrational cycle. For a more quantitative and transferable analysis of the electron localization caused by the moving nuclei, we define a function measuring the similarity between the induced density δ ρ(r,t) and the difference density ∆ρ (e) (r) = ρ e (r)−ρ g (r) of the considered excitation: where w(t) is a normalized window function with an extension such that the strong oscillations related to the coherent electronic motion are removed, but the trend of the curves is maintained. χ(t) is a dimensionless quantities corresponding to the excited-state population when the nuclei are kept fixed (further details in the Supporting Information). More generally, χ(t) measures the similarity between the TD density and the density of the vertically excited state. If the nuclei are kept fixed, only an electronic excitation contributes to χ(t); if they are released, the nuclear relaxation can cause an corresponding redistribution of charge. χ(t) = 0 indicates that the density displacements are uncorrelated with the excitation. Conversely, if χ(t) = 1, the electron density of the system is close to that of the molecule in a vertically excited state with 100% ES population. Comparing χ(t) obtained in the simulations with clamped and moving nuclei, we can determine to which extent the nuclear motion enhances the excitationinduced charge redistribution. We compare this result to an indicator for the deviation of the nuclear geometry from the equilibrium one: .
The individual atomic displacements |R J (t) − R J (0)| are weighted by the valence charge Z J of the corresponding pseudopotential, which is proportional to the electrostatic force exerted by the positive pseudo-nuclei and thus to the interaction strength with the electronic system. The quantity in Eq. (29) predicts quite well the positions of the maxima of χ(t) (see Fig. 8). In the case of ethylene [ Fig. 8a)], σ (t) essentially describes the change of the C=C bond length, as its HOMO-LUMO excitation is mainly coupled to vibrational modes involving a C=C stretch (see discussion above). Both χ(t) and σ (t) have a maximum at around 20 fs, with χ(t) exhibiting a peak at about 0.6 when the nuclei are in motion, in contrast to the constant value of 0.2 assumed after the laser illumination when the nuclei are clamped (see Fig. 7). Also in the case of thiophene, the nuclear motion enhances the value of χ(t) [Fig. 8b)]. The difference between fixed and moving nuclei is however less pronounced than in ethylene, which we attribute to the larger mass of S in comparison to H and C. The most pronounced discrepancy between fixed and moving nuclei is observed in benzene, for which χ(t) assumes a maximum value of about 1.5 when nuclei are moving, compared to 0.4 with fixed nuclei [Fig. 8c)]. The fact that χ(t) exceeds 1 in the former case is likely related to the fact that the pumped excitation is degenerate (see also Table S3).
To rationalize the increase of χ(t) upon the inclusion of the nuclear motion, we first note that the vibrational energies are insufficient to cause interband transitions. Thus, the ES population cannot be increased by the nuclear motion. We instead attribute the enhancement effect to the change of the GS electron density associated with a change of the nuclear configuration. The electronic excitation leads to a partial population of the excited-state PES, in a point which is not a minimum, such that forces begin to drive the nuclear configuration away from the GS equilibrium one (R). The new configuration (R + ∆R) has a different GS density, which is similar to that of the vertically excited state, i.e. ρ g (R + ∆R) ∼ ρ e (R). This relaxationrelated charge redistribution acts on top of the one stemming directly from the electronic excitation. As we can see in Fig.  8, the ratio of these two contributions depend on the system: In thiophene, the electronic excitation causes the most pronounced change in the density; in ethylene and particularly in benzene, this role is played by the nuclear motion.

IV. SUMMARY AND CONCLUSIONS
We have investigated the laser-induced dynamics of three small organic molecules (ethylene, benzene, and thiophene) by means of RT-TDDFT in the adiabatic local density approximation. The scope of our study was to introduce simple models and indicators to interpret the quantities computed with this method, and to better understand the physics of ultrafast dynamical processes. By exciting each molecule with a laser pulse in resonance with a low-lying electronic excitation, we have monitored the time evolution of the induced dipole moment and of the electronic population. The comparison with the results obtained from an exactly-solvable two-level model has provided a diagnostic tool to analyze the computed quantities and to interpret the corresponding results in a physically transparent manner. The analysis of the transient absorption spectra, computed with and without including the nuclear motion, has revealed the complex dynamics of laser-excited systems, even of such small size as those examined in this work. The appearance of excited-state absorption, induced by the laser pulse and probed at subsequent time delays, has been demonstrated by an ad hoc simulation triggered by a train of monochromatic pulses. Finally, the analysis of the coupled electron-nuclear motion has allowed us to disclose the effects of the external field on the vibrational degrees of freedom, and to reveal the crucial role of vibronic coupling in stabilizing the excitation already in the first few fs after the laser illumination. These effects have been shown by a visual representation of the induced charge density distribution as well as through general indicators applicable also to more complex molecular systems.
This work and the analysis reported therein offer useful tools to understand the results of RT-TDDFT in describing ultrafast dynamical processes within the given approximations. We expect that this contribution will further consolidate the role of RT-TDDFT, coupled to the Ehrenfest molecular dynamics scheme, as state-of-the-art methodology capable to address the fascinating and yet complex phenomena of lightmatter interaction in the femtosecond regime.

SUPPORTING INFORMATION
The Supporting Information includes further details about the linear-response calculations of the relevant excitations, a methodological section related to the excited-state absorption calculations of Sec. III B, additional computational details about the two-level model, and derivations of different formulas for the excited-state population from RT-ALDA.