We expand the concept of natural transition orbitals in the context of real-time time-dependent density functional theory (RT-TDDFT) and show its application in practical calculations. Kohn–Sham single-particle wavefunctions are propagated in RT-TDDFT simulation, and physical properties remain invariant under their unitary transformation. In this work, we exploit this gauge freedom and expand the concept of natural transition orbitals, which is widely used in linear-response TDDFT, for obtaining a particle–hole description in RT-TDDFT simulation. While linear-response TDDFT is widely used to study electronic excitation, RT-TDDFT can be employed more generally to simulate non-equilibrium electron dynamics. Studying electron dynamics in terms of dynamic transitions of particle–hole pairs is, however, not straightforward in the RT-TDDFT simulation. By constructing natural transition orbitals through projecting time-dependent Kohn–Sham wave functions onto occupied/unoccupied eigenstate subspaces, we show that linear combinations of a pair of the resulting hole/particle orbitals form a new gauge, which we refer to as dynamical transition orbitals. We demonstrate the utility of this framework to analyze RT-TDDFT simulations of optical excitation and electronic stopping dynamics in the particle–hole description.

Over the past few decades, time-dependent density functional theory (TDDFT), based on the Runge–Gross theorem,1 has become a computationally attractive approach for studying excitation and non-equilibrium dynamics of an electronic system in molecules and extended materials.2 Linear response TDDFT is now widely employed to calculate optical excitation and the corresponding excited state properties,3–5 and the real-time propagation approach to TDDFT (RT-TDDFT) has become increasingly popular for studying non-linear excitation dynamics and charge transport properties.6–10 While the increasingly sophisticated exchange-correlation approximation and computational power are promising for the more accurate prediction of excited state and dynamical properties for increasingly more complex systems, developing a physically intuitive understanding from the first-principles simulation is also desirable. In the context of linear-response TDDFT in the frequency domain, the one-particle transition density matrix (TDM) by Furche11 is widely used for interpreting electronic excitation.

The more recent work by Li and Ullrich further expanded and introduced the idea of a particle–hole map using the particle–hole transition density matrix.12 The one-particle TDM also serves as the basis of constructing the natural transition orbitals,13 which is widely used to describe individual electronic excitations approximately as transition events between a pair of natural transition orbitals.4,14 This framework attempts to effectively reduce an individual excitation of a many-body electronic system to a transition from a hole orbital to a particle orbital. While the linear response approximation works well for studying optical excitation, the RT-TDDFT approach enables us to study a much broader aspect of the inherently dynamical phenomenon of optical excitation of valence electrons15–20 and core electrons21–24 as well as other types of non-equilibrium dynamic phenomena, such as electronic stopping25–35 and dielectric breakdown.36–39 At the same time, RT-TDDFT simulations have received much less attention in terms of interpreting the electron dynamics from the viewpoint of electronic transitions. Li and Ullrich have formulated the time-dependent extension of the TDM, and its application to TDDFT was discussed by noting its potential usefulness to the real-time propagation method.40 Motivated by the idea of natural transition orbitals in linear-response TDDFT, we introduce dynamical transition orbitals for RT-TDDFT dynamics in this work.

The electron density remains invariant under a unitary transformation of Kohn–Sham orbitals in DFT. This gauge freedom, which is inherent also in the time-dependent Kohn–Sham (TD-KS) equations, has been exploited in the RT-TDDFT simulation in recent years. Under arbitrary unitary gauge transformations, physical properties of the electron dynamics are unchanged, but certain gauge transformations allow for reformulations of the TD-KS equations that are advantageous for some specific purposes.41 Lin and co-workers, for example, have used the parallel transport gauge so that oscillations of the propagating orbitals are minimized, allowing for significantly larger propagation time steps to be used with implicit time integration schemes.42–44 Kanai and co-workers have used the maximally localized Wannier function gauge in the RT-TDDFT simulation for treating the homogeneous electric field with the periodic boundary conditions (PBC) and for deducing the chemically intuitive understanding at the molecular level.45 Similarly, new dynamical transition orbitals are constructed in this work as a unitary transformation of the occupied TD-KS orbitals, and each dynamical transition orbital is formulated as a time-dependent linear combination of a pair of hole and particle orbitals. The general idea is to divide the eigenstate space into two complementary subspaces, and the particle and hole orbitals are obtained by performing the singular value decomposition (SVD) in each subspace. Mathematically, the choice of the complementary subspaces is not unique, and we use the subspaces spanned by the occupied KS eigenstates and the unoccupied (virtual) KS eigenstates so that transitions between the hole and particle orbitals can be interpreted as single-particle excitations. Using the property that each dynamical transition orbital is composed of a single pair of hole and particle orbitals, the time-dependence of the hole/particle populations help us gain physical insight into electronic excitation phenomena from the RT-TDDFT dynamics. Because the dynamical transition orbitals are a unitary transformation of the underlying TD-KS orbitals, physical properties remain invariant in the transformation. We note that time-dependent density matrix functional theory exploits a related idea, for which time-dependent orbitals termed natural orbitals are generated by diagonalizing the one-particle density matrix.46–50 

The rest of this work is organized as follows. In Sec. II, we formulate dynamical transition orbitals as a unitary transformation of the TD-KS orbitals such that time-dependent linear combination of a pair of hole and particle orbitals yields individual transitions. The equations of motion are then derived as they show when discontinuities might arise and how the dynamic continuity of the dynamical transition orbitals can be enforced. In Sec. III, we discuss the linear response approximation to the present formalism. In Sec. IV, we show the numerical implementation and practical application of the dynamical transition orbitals in a few different types of electron dynamics in several real chemical systems.

We introduce dynamical transition orbitals, ψnDTr,tn=1..N, as a unitary transformation of N occupied time-dependent Kohn–Sham (TD-KS) wave functions. We will show in the following that the dynamical transition orbitals can be constructed as a linear combination of the hole and particle orbitals,

|ψnDT(t)=an(t)|ψnh(t)+bn(t)|ψnp(t),
(1)

where ψnh(r,t) and ψnp(r,t) are the hole and particle orbitals that derive from the singular value decomposition (SVD) of the occupied and unoccupied (virtual) eigenstate subspaces, respectively. The coefficients (which are referred to as the occupation number here) satisfy

an(t)0bn(t)0an(t)2+bn(t)2=1.
(2)

The dynamical transition orbitals are a unitary transformation of the TD-KS wavefunctions, and the one-particle density matrix (i.e., density operator) can be expressed as

Γ^(t)=nN|ϕn(t)ϕn(t)|=nN|ψnDT(t)ψnDT(t)|,
(3)

where {ϕn(r, t)} are the TD-KS wavefunctions. Physically observable Â(t)=Tr[ÂΓ^(t)] are invariant under the unitary transformation. For a concise notation, we omit the spin degrees of freedom in the following proof, and thus, N is equal to Nelec/2. Key properties of the projection operator and density operator we employ for the following proof are briefly summarized in Subsection 1 of the  Appendix.

1. Hole orbitals

As done in the work on the time-dependent transition density matrix by Li and Ullrich,40 the state-to-state transition amplitude matrix is defined by projecting time-dependent Kohn–Sham wavefunctions onto the ground-state Kohn–Sham eigenstates,

Sij(t)=ψi|ϕj(t),
(4)

where ψi(r) are the eigenstates and ϕj(r, t) are the time-dependent wavefunctions. Then, a new operator can be defined via

Ŝ(t)=ij|ψiSijϕj(t)|=i|ψiψi|j|ϕj(t)ϕj(t)|=P^oΓ^(t).
(5)

Noting the idempotent property of the projection operator P^o and that the density matrix is Hermitian, Eq. (5) yields

Ŝ(t)Ŝ(t)=Γ^(t)P^oP^oΓ^(t)=Γ^(t)P^oΓ^(t).
(6)

Similar to the static natural transition orbitals,13 we obtain the dynamical transition orbitals and the hole orbitals via the SVD of the Ŝ operator,

Ŝ(t)=n|ψnh(t)an(t)ψnDT(t)|,
(7)

where an(t) are the singular values. Using the SVD property Eq. (A9) in the  Appendix and Eq. (6), we have

Ŝ(t)Ŝ(t)|ψnDT(t)=Γ^(t)P^oΓ^(t)|ψnDT(t)=an(t)2|ψnDT(t).
(8)

Similarly, using the complementary property Eq. (A10) in the  Appendix yields the hole orbitals via

Ŝ(t)Ŝ(t)|ψnh(t)=P^oΓ^(t)P^o|ψnh(t)=an(t)2|ψnh(t).
(9)

2. Particle orbitals

For the particle orbitals, we project the time-dependent Kohn–Sham wavefunctions onto the unoccupied Kohn–Sham eigenstates,

Taj(t)=ψa|ϕj(t).
(10)

Analogous to Eq. (5), we can write the operator

T^(t)=aj|ψaTajϕj(t)|=P^vΓ^(t).
(11)

Then, the dynamical transition orbitals and the particle orbitals can be obtained via the SVD

T^(t)=n|ψnp(t)bn(t)ϕnDT(t)|.
(12)

Here, the dynamical transition orbitals corresponding to the particle orbitals in Eq. (12) are not necessarily the same as the dynamical transition orbitals corresponding to the hole orbitals [Eq. (7)], and therefore, we used the |ϕnDT(t) notation instead of |ψnDT(t). Analogous to the hole orbital case [Eqs. (8) and (9)], we have

Γ^(t)P^vΓ^(t)|ϕnDT(t)=bn(t)2|ϕnDT(t)
(13)

and

P^vΓ^(t)P^v|ϕnp(t)=bn(t)2|ϕnp(t).
(14)

3. Equivalence between |ϕnDT(t) and |ψnDT(t)

We show here that the hole orbital and the particle orbital share the same dynamical transition orbital. Using the relationship between the occupied and unoccupied projection operators [Eq. (A7) in the  Appendix], Eq. (8) can be written as

Γ^(t)(ÎP^v)Γ^(t)|ψnDT(t)=an(t)2|ψnDT(t).
(15)

The term with the identity operator on the left-hand side can be rearranged as

(Γ^(t)ÎΓ^(t))|ψiDT(t)=Γ^(t)|ψiDT(t)=m|ψmDT(t)ψmDT(t)|ψnDT(t)=mδm,n|ψmDT(t)=|ψnDT(t),
(16)

and thus, Eq. (15) can be written as

Γ^(t)P^vΓ^(t)|ψnDT(t)=(1an(t)2)|ψnDT(t).
(17)

Comparing this equation to Eq. (13), one recognizes that both |ϕnDT(t) and |ψnDT(t) are the same set of eigenvectors (up to an arbitrary phase) for the operator Γ(t)PvΓ(t). The corresponding eigenvalues equate bn(t)2 and 1 − an(t)2, and Eq. (2) is proved. Note that an(t) and bn(t) are real positive numbers since they are singular values of the SVD [see Eq. (A8) in the  Appendix].

4. Dynamical transition orbitals from particle and hole orbitals

Operating the density operator on the dynamical transition orbitals and using Eq. (A7) in the  Appendix, we have

Γ^(t)|ψnDT(t)=(P^o+P^v)Γ^(t)|ψnDT(t)=P^oΓ^(t)|ψnDT(t)+P^vΓ^(t)|ψnDT(t).
(18)

The left-hand side can be expanded as

m|ψmDT(t)ψmDT(t)|ψnDT(t)=mδm,n|ψmDT(t)=|ψnDT(t).
(19)

Using Eqs. (5) and (7) for the first term and Eqs. (11) and (12) for the second term, the right-hand side of Eq. (18) is expanded as

m|ψmh(t)an(t)ψmDT(t)|ψnDT(t)+   m|ψmp(t)bn(t)ψmDT(t)|ψnDT(t)  =mδm,nan(t)|ψmh(t)+mδm,nbn(t)|ψmp(t)  =an(t)|ψnh(t)+bn(t)|ψnp(t).
(20)

Hence, Eq. (18) becomes

|ψnDT(t)=an(t)|ψnh(t)+bn(t)|ψnp(t).
(21)

We now have established that the dynamical transition orbitals can be expressed as a linear combination of the hole and particle orbitals with the coefficients satisfying Eq. (2).

In the basis of the dynamical transition orbitals, the one-particle density matrix for a time-dependent Kohn–Sham system is

Γ^(t)=n|ψnDT(t)ψnDT(t)|.
(22)

Using Eq. (21), it can be expanded in terms of the particle and hole orbitals as

Γ^(t)=nan(t)2|ψnh(t)ψnh(t)|+bn(t)2|ψnp(t)ψnp(t)|+an(t)bn(t)|ψnp(t)ψnh(t)|+an(t)bn(t)|ψnh(t)ψnp(t)|.
(23)

By expressing the density matrix in the basis of |ψnh(t) and |ψnp(t), we may regard the resulting density matrix as being composed of four sub-matrices (hole–hole, particle–hole, hole–particle, and particle–particle),

ψmh(t)Γ^(t)|ψnh(t)=δm,nan(t)2,ψmp(t)Γ^(t)|ψnp(t)=δm,nbn(t)2,ψmh(t)Γ^(t)|ψnp(t)=δm,nan(t)bn(t),ψmp(t)Γ^(t)|ψnh(t)=δm,nan(t)bn(t).
(24)

The density matrix is composed of the four block matrices with each sub-matrix being diagonal. At each instance of time, in the particle–hole basis, the dimension of the density matrix is 2N, given by N particle orbitals and N hole orbitals. For later discussion, we note that the density operator is expressed in the continuous position basis as

Γr,r,t=r|Γ^t|r=nant2ψnh*r,tψnhr,t+bnt2ψnp*r,tψnpr,t+antbntψnp*r,t×ψnhr,t+antbntψnh*r,tψnpr,t.
(25)

The electron density is given by the diagonal element of the density matrix,

ρr,t=drδrrΓr,r,t=nant2|ψnh*r,t|2+bnt2|ψnp*r,t|2+antbnt×ψnp*r,tψnhr,t+antbntψnh*r,tψnpr,t=nant2|ψnh*r,t|2+bnt2|ψnp*r,t|2+2antbntReψnp*r,tψnhr,t.
(26)

Having introduced our formulation of dynamical transition orbitals within the framework of real-time TDDFT, it is instructive to derive the equation of motion for the particle and hole orbitals in our efforts to gain physical insight into the dynamics of these orbitals and their occupations. Given the definitions Eqs. (A1) and (A2) in the  Appendix, the projection operators do not change in time

P^̇o=0,P^̇v=0.
(27)

Recall that the occupation numbers can be given by the following eigenvalue equation:

P^oΓ^(t)P^o|ψnh(t)=an(t)2|ψnh(t).
(28)

For convenience, we define the population of hole orbital as σnh(t)an(t)2. Taking the first derivative with respect to time on both sides of Eq. (28) yields

P^oΓ^̇(t)P^o|ψnh(t)+P^oΓ^(t)P^o|ψ̇nh(t)=σ̇nh(t)|ψnh(t)+σnh(t)|ψ̇nh(t).
(29)

Multiplying both sides by ψmh(t)|, we have

ψmh(t)|P^oΓ^̇(t)P^o|ψnh(t)+ψmh(t)|P^oΓ^(t)P^o|ψ̇nh(t)  =σ̇nh(t)ψmh(t)|ψnh(t)+σnh(t)ψmh(t)|ψ̇nh(t).
(30)

Note that the second term on the left-hand side (being Hermitian operator) is

ψmh(t)|PoΓ(t)Po|ψ̇nh(t)=σmh(t)ψmh(t)|ψ̇nh(t).
(31)

Equation (30) can be written as

ψmh(t)|PoΓ̇(t)Po|ψnh(t)=σ̇nh(t)ψmh(t)|ψnh(t).
(32)

By setting m = n in Eq. (32) and interchanging the left-hand and right-hand sides, we arrive at the equation of motion for the hole population,

σ̇nh(t)=ψnh(t)|PoΓ̇(t)Po|ψnh(t).
(33)

Similarly, for the particle population, we have

σ̇np(t)=ψnp(t)|PvΓ̇(t)Pv|ψnp(t).
(34)

When mn, because of the orthogonality, σ̇nh(t)ψmh(t)|ψnh(t) vanishes in Eq. (30), and the time derivative of the hole orbital can be written as an expansion of other occupied hole orbitals,

|ψ̇nh(t)=mnNocc..ψmh(t)|P^oΓ^̇(t)P^o|ψnh(t)σnh(t)σmh(t)|ψmh(t).
(35)

Similarly, for the particle orbitals, we have

|ψ̇np(t)=mnNvir.ψmp(t)|P^oΓ^̇(t)P^o|ψnp(t)σnp(t)σmp(t)|ψmp(t),
(36)

where the summation index m′ is used to distinguish from m. While Nocc is the number of occupied orbitals, Nvir is the number of unoccupied (virtual) orbitals, which is finite in practical calculations. Although the ψnDT(t) remain as unitary transforms of the underlying TD-KS wavefunctions, the vector space spanned by the particle orbitals is not necessarily invariant in time. As shown in the above equations of motion [Eqs. (35) and (36)], one must pay special attention to a possibility of having discontinuity when σnh/p(t)=σmh/p(t) in time. For example, as shown in Fig. 1(a), the particle populations could become equal at a number of instances in time for the two dynamical transition orbitals. Near these individual crossing points, assigning a particular state index to a specific particle/hole orbitals becomes numerically rather arbitrary. However, since the particle/hole orbitals must change continuously in time, we can identify the correct state index assignment by requiring the populations of individual particle/hole orbitals to change smoothly. The state index assignment as indicated by two different colors in Fig. 1(a) was made possible by employing this additional requirement. The observation of such crossing points stems from having degenerate KS eigenstates in some systems such as polyacetylene we simulated here. When we calculate the projection of one of the two particle orbitals onto KS eigenstates [i.e., |ϕa|ψnp(t)|2], we find sharp peaks at these crossing points for the two degenerate KS eigenstates that give dominant contributions to the specific particle orbital, as shown in Fig. 1(b).

FIG. 1.

Dynamical transition orbitals in the RT-TDDFT simulation of a trans-polyacetylene with an applied external potential as an example of how the time derivative of the particle/hole orbitals and their populations change smoothly. (a) The particle population of two dominant transition modes. (b) The projection of the particle orbital for transition mode 1 [shown in red in (a)] onto the Kohn–Sham eigenstates with the two most dominant contributions.

FIG. 1.

Dynamical transition orbitals in the RT-TDDFT simulation of a trans-polyacetylene with an applied external potential as an example of how the time derivative of the particle/hole orbitals and their populations change smoothly. (a) The particle population of two dominant transition modes. (b) The projection of the particle orbital for transition mode 1 [shown in red in (a)] onto the Kohn–Sham eigenstates with the two most dominant contributions.

Close modal

The dynamical transition orbitals are a gauge transformation of the underlying time-dependent Kohn–Sham wavefunctions by design, and the time evolution of the electron density governs the system’s dynamics in the real-time TDDFT simulation. The one-particle density matrix is propagated in the time-dependent density matrix functional theory. The so-called natural orbitals are obtained by diagonalizing the density matrix as its eigenstates,46 and the eigenvalues (i.e., occupation numbers) are not one in general.49 In real-time TDDFT, on the other hand, the eigenvalues of the one-particle density matrix are one by construction. The basic idea behind our dynamical transition orbital approach is to divide the Hilbert space into the occupied and unoccupied subspaces, and we obtain the dynamical transition orbitals as the eigenstates of the one-particle density matrix with the projection onto the subspaces [Eqs. (8) and (13)]. Then, the eigenvalues, or equivalently the populations [i.e., σnh(t)=an(t)2], yield the expansion coefficients with 0 ≤ an(t), bn(t) ≤ 1 for the time-dependent hole and particle orbitals that form the dynamical transition orbitals as their linear combination [Eq. (21)]. Thus, there is no net gain or loss but only changes in the population between the hole orbital and particle orbital for individual dynamical transition orbitals in real-time dynamics [i.e., σnh(t)+σnp(t)=1]. If the hole and particle orbitals themselves are invariant in time (up to a phase), then the population changes give dynamical information on the nature of the electronic transition in terms of a convenient particle–hole picture. Such a situation is possible if a external perturbation is sufficiently weak, and we explore this aspect in the context of the linear response approximation in the following.

In Sec. II, the particle–hole orbital dynamics are formulated explicitly in time by introducing the dynamical transition orbitals. In order to connect to the linear response theory, which is often used to study optical excitation, we consider a small perturbation to the Hamiltonian. This entails to having a sufficiently small population increase for the particle orbital such that bn(t)2 ≈ 0. Due to Eq. (2), we have 0 ≤ an(t) ≤ an(t)2 with an(t) ≈ 1. Under this approximation of the weak perturbation, the electron density [Eq. (26)] and the one-particle density operator [Eq. (25)] yield

ρ(r,t)=n|ψnh(r,t)|2+bn(t)ψnp*(r,t)ψnh(r,t)+ψnh*(r,t)ψnp(r,t).
(37)

Then, the density change (with respect to the equilibrium state) is given by changes in the hole–particle and particle–hole submatrices of the density operator,

δρ(r,t)=ρ(r,t)ρeq(r)=nbn(t)ψnp*(r,t)ψnh(r,t)+ψnh*(r,t)ψnp(r,t).
(38)

If the particle and hole orbitals are invariant in time (assuming a weak perturbation), the Fourier transform allows us to write the density change as the product of two terms that are separately independent of time and space,

δρ(r,ω)=nbn(ω)ψnp*(r)ψnh(r)+ψnh*(r)ψnp(r).
(39)

We now connect this dynamical transition orbital expression [Eq. (39)] to the many-body expression and show how the particle occupation number bn(ω) contains information about the excitation energy. By truncating the evolution operator in the interaction picture at the first-order term in the Dyson series, the density change can be written as

δρr,t=dtdrχr,t;r,tvextr,t,
(40)

where χ is the density–density response function that connects the density response to an external perturbation (e.g., applied electric field). The frequency-dependent density response is obtained via Fourier transform as

δρr,ω=χr,r,ωvextr,ωdr.
(41)

The response function, with causality taken into account, in the Lehmann representation is given by

χ(r,r,ω)=limη0nΨ0|ρ^(r)|ΨnΨn|ρ^(r)|Ψ0ω(EnE0)+iηΨ0|ρ^(r)|ΨnΨn|ρ^(r)|Ψ0ω+(EnE0)+iη,
(42)

where Ψn is the many-body state in the n-th excited state and the operator is ρ^(r)=iNelec.δ(rri).2 Note that in the linear-response formulation of TDDFT (LR-TDDFT), χ(r, r′, ω) is related to the non-interacting Kohn–Sham response function via the exchange-correlation kernel. Equation (41) now can be written as

δρ(r,ω)=limη0nΨ0|ρ^(r)|ΨnΨn|vext(ω)|Ψ0ω(EnE0)+iηΨ0|vext(ω)|ΨnΨn|ρ^(r)|Ψ0ω+(EnE0)+iη.
(43)

By assuming that transition density Ψ0|ρ^(r)|Ψn is real and also that a single pair of particle and hole orbitals dominantly contributes to the transition density (see the  Appendix) such that

Ψ0|ρ^r|Ψn=Ψn|ρ^r|Ψ0=Reψnh*rψnpr=12ψnh*rψnpr+ψnhrψnp*r,
(44)

comparing Eq. (39) to Eq. (43) shows that the frequency-dependent term is equal to the occupation number of the particle orbital,

bn(ω)=2limη0Ψn|vext(ω)|Ψ0ω(EnE0)+iηΨo|vext(ω)|Ψnω+(EnE0)+iη.
(45)

Therefore, the poles in bn(ω) indicate the system’s excitation energy EnE0 of the many-body expression. This perspective on the particle occupation number is demonstrated in the following in a numerical example.

The dynamical transition orbital approach can be implemented on top of most RT-TDDFT formalisms based on the TD-KS equation. We employ the RT-TDDFT in the plane wave pseudopotential formalism as discussed in Refs. 51 and 52. The RT-TDDFT method has been implemented in the Qb@ll branch of Qbox code.53,54 While it is, in principle, possible to directly propagate the dynamical transition orbitals according to the equations of motion for the particle/hole orbitals [Eqs. (35) and (36)], Sec. II C makes it clear that such a procedure is not practical. Instead, the dynamical transition orbitals are obtained at each time step as a unitary transform of time-dependent KS wavefunctions in the following way:

  • The occupied KS eigenstates ψi are calculated, while the unoccupied KS eigenstates ψa are not necessary.

  • ϕn(t) or Γ^(t) are propagated in the RT-TDDFT simulation. Note that it is possible to use any gauge in the propagation, including the time-dependent maximally localized Wannier functions45 and the parallel transport gauge orbitals.42 

  • At each time step, an(t) and |ψnDT(t) are obtained by solving the eigenvalue problem

Γ^(t)P^oΓ^(t)|ψnDT(t)=an(t)2|ψnDT(t).
(46)
  • |ψnh(t) are obtained by acting the projection operator for the occupied subspace,

P^o|ψnDT(t)=an(t)|ψnh(t).
(47)
  • The particle orbitals are obtained via

|ψnp(t)=1bn(t)ψnDT(t)an(t)|ψnh(t),
(48)

where bn(t) can be obtained from bn(t)2 = 1 − an(t)2 [Eq. (2)].In some cases, particularly when bn(t) is very close to 0, Eq. (48) suffers from an ill-conditioned problem in numerical calculations as we have observed in some cases discussed in the following. We circumvent this numerical problem by defining an effective projection operator that spans smaller numbers of unoccupied KS eigenstates,

P^veff=aNeff|ψaψa|,
(49)

where Neff defines the smaller subspace spanned by unoccupied KS eigenstates ψa. In practice, this effective space can be chosen by gradually increasing the unoccupied Kohn–Sham eigenstates until a convergence is reached such that the particle orbital does not change noticeably. The particle orbitals can be then obtained via

|ψnp(t)P^veff|ψnDT(t)
(50)

and by enforcing the proper normalization. The procedure is concisely summarized in a flowchart, as shown in Fig. 2. In the following applications of the dynamical transition orbitals approach, we note it when this approximation is used. As we study electronic excitation in optical excitation and in electronic stopping, we analyze each dynamical transition orbital also by projecting constituting particle/hole orbitals onto the reference Kohn–Sham eigenstates,

Di,n(t)=|ψi|ψnh(t)|2,
(51)
Da,n(t)=|ψa|ψnp(t)|2.
(52)
FIG. 2.

Flowchart for the numerical procedure of obtaining dynamical transition orbitals, hole orbitals, and particle orbitals in RT-TDDFT dynamics.

FIG. 2.

Flowchart for the numerical procedure of obtaining dynamical transition orbitals, hole orbitals, and particle orbitals in RT-TDDFT dynamics.

Close modal

We present a few examples of how the dynamical transition orbitals are used for interpreting RT-TDDFT dynamics in the particle–hole transition viewpoint. We study the dynamics of valence electrons in both optical excitation and electronic stopping excitation. Although this approach is equally applicable to studying core-electron excitation in principle, our use of plane wave pseudopotential formulation makes studies of core electrons quite challenging in practice despite some recent progress.35 

Optical excitation can be simulated and the absorption spectrum can be obtained within the RT-TDDFT framework by applying external potential perturbation.15–18 Based on the time-dependent Maximally Localized Wannier Functions (TD-MLWFs) with the length gauge, the external electric field can be applied as an additional scalar potential in the Hamiltonian to both isolated molecule and extended systems on an equal footing.45 In order to study the specific region of the optical excitation spectrum, the Gaussian envelope potential is used,

E(t)=E0cos(ω0t)eγ(tt0)2,
(53)

where ω0 is the specific frequency the excitation is centered on. γ controls the Gaussian width of the external potential applied. The amplitude E0 is set sufficiently small such that the perturbative response of the electronic system is in the linear response regime.

1. Isolated molecules

For isolated molecule cases, the water molecule and benzene molecule are considered here to illustrate how the dynamical transition orbital formulation is used to study the electronic excitation in real-time dynamics. The simulation box is set to 80 × 80 × 80 bohrs3 such that the interaction between the molecule and its periodic image is negligible. The Gaussian envelope potential is applied along the x direction with ω0 = 6 eV, t0 = 75 a.u., γ = 2(eV/)2, and E0 = 0.01 a.u. All atoms are represented by Hamann–Schluter–Chiang–Vanderbilt (HSCV) norm-conserving pseudopotentials,55,56 and the local-density approximation (LDA) XC functional is employed with a plane wave cutoff energy of 30 Ry. The integration time step of 0.1 a.u. is used with the explicit fourth-order Runge–Kutta integrator. After the initial application of the external potential (for the time duration of 150 a.u.), the simulation is continued for another 600 a.u. For the particle orbitals, we employ Eq. (50) with the lowest ten unoccupied KS eigenstates in the projection.

For brevity, the term “transition mode” is used here to refer to the time-dependent change in a specific dynamical transition orbital indexed by n, which is a linear combination of a pair of hole/particle orbitals. For the water molecule, there is a well-defined excitation peak at 5.67 eV in the absorption spectrum, corresponding to the lowest excited state. For the benzene molecule, the first prominent excitation peak at 6.61 eV derives from two degenerate excited states. While it is conceptually not straightforward to probe this quasi-degeneracy in the RT-TDDFT simulation, the Casida equation clearly shows the two near-degenerate excited state solutions within the linear-response TDDFT.57 For the water molecule, the dynamical transition orbitals show that a single transition orbital is largely responsible for the observed electronic response, as shown in Fig. 3(a). For the benzene molecule, two transition orbitals dominate the excitation process. Other transition modes contribute negligibly to the electronic excitation process, as shown in Fig. 3(b). The particle and hole orbitals can be projected independently onto the Kohn–Sham eigenstates [see Eqs. (51) and (52)]. As shown in Table I, for the water molecule, the dynamical transition orbital is dominantly of HOMO → LUMO transition. For the benzene molecule, the two dominating dynamical transition orbitals derive from HOMO → LUMO +1 and HOMO −1 → LUMO transitions. HOMO/HOMO −1 and LUMO/LUMO −1 are essentially degenerate, with the energy difference being less than 0.001 eV. As the external perturbing field is applied to the systems, the particle populations start to increase with a characteristic oscillation, as shown in Fig. 3. According to Eq. (45), this oscillation could be directly linked to the many-body excitation energy within the linear response approximation. As shown in Table II, the peak positions in the absorption spectrum are indeed nearly identical to the peaks in the Fourier transform of the time-dependent particle occupation number (see the supplementary material for information on how the particle orbital occupation changes with the optical pulse electric field). In the derivation of Eq. (45), it was assumed that the hole/particle orbitals are invariant to time under a weak perturbation [Eq. (39)]. Figure 4 shows time-dependent changes in the projections Di,n(t) and Da,n(t) [Eqs. (51) and (52)] for examining the extent to which this assumption is valid. After the initial response of the equilibrium system to the external perturbation, the projections become essentially unchanged in time after 50 a.u. The variation of the projections afterward is indeed less than 0.1% for the rest of the entire simulation that lasts for more than 700 a.u., and the assumption of time-invariant particle/hole orbitals is valid here. While advantages of the dynamical transition orbitals may not be obvious for these simple molecular cases here, it shows that the formalism reduces to the natural transition orbital description when the particle/hole orbitals are invariant in time. We now consider more complicated excitation dynamics for an optically excited condensed matter system and for the electronic stopping in order to illustrate more obvious benefits of employing the dynamical transition orbital gauge in analyzing RT-TDDFT simulations.

FIG. 3.

Particle population [i.e., σnp(t)bn(t)2] change for (a) a water molecule and for (b) a benzene molecule. Four most significant transition modes are shown.

FIG. 3.

Particle population [i.e., σnp(t)bn(t)2] change for (a) a water molecule and for (b) a benzene molecule. Four most significant transition modes are shown.

Close modal
TABLE I.

Time-averaged projections of time-dependent particle/hole orbitals onto a Kohn–Sham eigenstate, Di,n(t) and Da,n(t). i refers to the occupied Kohn–Sham eigenstates, and a refers to unoccupied Kohn–Sham eigenstates. HOMO/HOMO −1 and LUMO/LUMO −1 are essentially degenerate, with the energy difference being less than 0.001 eV.

MoleculeModeiDi,n(t)⟩%aDa,n(t)⟩%
H2HOMO 99.8 LUMO 94.4 
C6H6 HOMO 97.3 LUMO+1 96.8 
 HOMO-1 99.1 LUMO 99.3 
MoleculeModeiDi,n(t)⟩%aDa,n(t)⟩%
H2HOMO 99.8 LUMO 94.4 
C6H6 HOMO 97.3 LUMO+1 96.8 
 HOMO-1 99.1 LUMO 99.3 
TABLE II.

The peak positions of the Fourier-transformed particle occupation and in the optical absorption spectrum corresponding to the lowest excitation energies.

MoleculeModebn(ω) (eV)Absorption spectrum(eV)
H25.73 5.67 
C6H6 1 and 2 6.61 6.61 
MoleculeModebn(ω) (eV)Absorption spectrum(eV)
H25.73 5.67 
C6H6 1 and 2 6.61 6.61 
FIG. 4.

Time evolution of the projection of the particle and hole orbitals on the Kohn–Sham eigenstates with the two most significant contributions, Di,n(t) and Da,n(t) [Eqs. (51) and (52)].

FIG. 4.

Time evolution of the projection of the particle and hole orbitals on the Kohn–Sham eigenstates with the two most significant contributions, Di,n(t) and Da,n(t) [Eqs. (51) and (52)].

Close modal

2. Extended system

The RT-TDDFT simulation is suitable also for studying extended systems with the periodic boundary conditions (PBC). Most condensed matter systems exhibit dense manifolds in the electronic density of states, and correspondingly, the absorption spectrum shows continuous excitation energy. In this subsection, we demonstrate how the dynamical transition orbitals can be used to describe optical excitation with a minimal number of single-particle orbitals unlike for the description given in terms of Kohn–Sham eigenstates. As an example, we excite crystalline silicon at the first prominent peak around 2.5 eV in the absorption spectrum (according to TDDFT based on the LDA XC functional).45 We employ a 32-atom supercell (41.0525 × 10.2631 × 10.2631 bohrs3) with the periodic boundary conditions (PBC) with the Γ-point approximation as an example for illustrating the approach. The supercell is elongated with four silicon unit cells in the x direction, along which the optical pulse [using the Gaussian envelope potential, Eq. (53)] is applied with the parameters ω0 = 2.5 eV, t0 = 200 a.u., γ = 1(eV/)2, and E0 = 0.003 a.u. The length gauge is used for modeling the corresponding spatially-homogeneous electric field with the PBC by employing the TD-MLWF approach.45 All atoms are represented by HSCV norm-conserving pseudopotentials, and the LDA XC functional is employed with a plane wave cutoff energy of 40 Ry. The integration time step of 0.05 a.u. is used with the explicit fourth-order Runge–Kutta integrator. The total simulation time is 400 a.u., and then, the system is propagated for another 150 a.u. for reaching a steady state.

TD-MLWFs are projected onto the Kohn–Sham eigenstates to generate projection matrices at instances of time. By taking the difference of the projection matrix amplitudes at time t and the initial time, it shows how the TD-MLWF changes in response to the external perturbation in the excitation process. Figure 5 shows the changes at t = 200 a.u. and at t = 400 a.u., well after the perturbing electric field ceases. They show that all TD-MLWFs respond to the external potential rather similarly as perhaps expected, yielding the hole population near HOMO and the excited electron population at energies much higher than LUMO. By adapting the dynamical transition orbitals, one can achieve a more compact and illustrative representation of the electronic excitation in terms of the particle–hole description. Figure 6 shows that the same excitation dynamics is described by changes in only three dynamical transition orbitals (transition modes 1–3). Two other transition modes also show a weak oscillation, but the corresponding particle population decay back to zero. Unlike the isolated molecule cases, the corresponding projections of the particle/hole orbitals on the KS eigenstates vary quite substantially in time, as shown in Fig. 7. Individual hole orbitals initially consist of several Kohn–Sham eigenstates, while particle orbitals are somewhat more singular in character. As the optical pulse field ceases (i.e., t ≈ 400 a.u.), particle/hole orbitals tend to be dominated by a single Kohn–Sham eigenstate. For all the three transition modes, the hole is generated among the Kohn–Sham eigenstates near HOMO (i.e., valence band maximum), while the excited electron resides at the energies much higher than LUMO (i.e., conduction band minimum). HOMO-1/HOMO-2 and LUMO+23/LUMO+24 are essentially degenerate between them, and the transition modes 1 and 2 yield a very similar response to the applied field (see Fig. 6). Unlike in the case of the molecular systems, the particle/hole orbitals are far from being time-invariant in this case. At the same time, using the dynamical transition orbital gauge, the optical excitation dynamics is described by the particle–hole dynamics of only three transition modes. There is a clear advantage over the TD-MLWF for representing the excitation dynamics in this particle-hole description.

FIG. 5.

Changes in the projection matrix of TD-MLWFs onto Kohn–Sham eigenstates at t = 200 a.u. and t = 400 a.u., with respect to the equilibrium ground state, for the crystalline silicon. The y-axis represents the Kohn–Sham index and the x-axis represents the TD-MLWF index.

FIG. 5.

Changes in the projection matrix of TD-MLWFs onto Kohn–Sham eigenstates at t = 200 a.u. and t = 400 a.u., with respect to the equilibrium ground state, for the crystalline silicon. The y-axis represents the Kohn–Sham index and the x-axis represents the TD-MLWF index.

Close modal
FIG. 6.

Particle population [i.e., σnp(t)bn(t)2] changes in the five most dominant transition modes for crystalline silicon excited at 2.5 eV.

FIG. 6.

Particle population [i.e., σnp(t)bn(t)2] changes in the five most dominant transition modes for crystalline silicon excited at 2.5 eV.

Close modal
FIG. 7.

Projection of the time-dependent particle/hole orbitals onto the Kohn–Sham eigenstates for the three most dominant transition modes for crystalline silicon. The relative magnitudes of the projection matrix are shown in the color map, ranging from 0% to 100%. The x-axis is time, and the y-axis is the KS (occupied/unoccupied) eigenstates index.

FIG. 7.

Projection of the time-dependent particle/hole orbitals onto the Kohn–Sham eigenstates for the three most dominant transition modes for crystalline silicon. The relative magnitudes of the projection matrix are shown in the color map, ranging from 0% to 100%. The x-axis is time, and the y-axis is the KS (occupied/unoccupied) eigenstates index.

Close modal

When a charged particle moves within or near a system at a very high speed (e.g., the proton kinetic energy of keV ∼ MeV), its kinetic energy can be partly transferred to the electronic subsystem of the matter through excitation and ionization of electrons. This non-equilibrium dynamics is responsible for the so-called electronic stopping. The RT-TDDFT simulation has been successfully employed to study various electronic stopping phenomena for ion irradiation in recent years.33 At the same time, interpreting the RT-TDDFT simulation of the non-equilibrium process is not straightforward. Recently, the new TD-MLWF approach has made it possible to interpret the dynamics in terms of chemical moieties.34,45 The dynamical transition orbital approach here can be also used effectively to study the non-equilibrium processes involved in the electronic stopping in terms of the electronic excitation. Following the earlier work of the TD-MLWF methodology,45 we here study the electronic stopping process of a proton projectile with a benzene molecule. The impact parameter (distance between the proton path and benzene) is set to be 8 a.u. We considered the near-peak stopping power velocity of 0.625 a.u., which induces a very large energy transfer, causing strong electronic excitation in benzene. The simulation is performed using the enforced time-reversal symmetry (ETRS) integrator58 with the time step of 0.1 a.u. Note that the dynamical transition orbital approach is independent of the particular integration scheme used. The LDA XC functional was used, and HSCV pseudopotentials are used with the plane wave cutoff of 30 Ry. The simulation cell of 50 × 30 × 30 bohrs3 is used. Note that the proton is completely ionized at t = 0 before the energy transfer takes place.

Figure 8(a) shows that the total hole population loss [i.e., σnh(t)] as a function of time, together with the particle orbital populations for the three dominant transition modes. Other transition modes contribute negligibly to the excitation. Figures 9(b)9(d) shows the corresponding hole and particle orbitals for these three distinct transition modes after the projectile proton passes by. Qualitatively speaking, the electron is excited from the hole orbital to the particle orbital in individual transition modes (i.e., a → b, c → d, and e → f). The extent of the excitation in each transition mode is quantified by the particle population σnp(t) in Fig. 8(a). The electronic stopping process often involves both electronic excitation within the target matter (i.e., benzene) and charge transfer between the target matter and the projectile ion (i.e., proton). Figure 9 shows that there are two types of transition modes. In transition modes 1 and 2, both electron excitation within the benzene molecule and the charge transfer from the molecule to the projectile proton take place. However, transition mode 3 facilitates only the charge transfer. For more detailed analysis, we project the dynamical transition orbitals onto the Kohn–Sham eigenstates of the benzene molecule for each transition mode, as shown in Figs. 8(b)8(d) and in Table III. Figures 8(b)8(d) show that each hole orbital is mainly composed of a single Kohn–Sham eigenstate of the benzene molecule (HOMO, HOMO-1, or HOMO-4) after the projectile proton passes by (t > 20 a.u.). These hole orbitals are largely invariant in time and the variance of which is around 0.1% (Table III). For the particle orbitals, the situation is quite different. The particle orbitals in the transition modes 1 and 2 derive noticeably from LUMO and LUMO+1, which are localized on the benzene molecule. The rest of the contribution comes from the eigenstates localized on the projectile proton [see Figs. 9(b) and 9(d)]. It is interesting to note that transition mode 2 involves the charge transfer to the projectile proton in more of the p-type orbital unlike for transition mode 1. Unlike the hole orbitals, the contributions from LUMO and LUMO+1 to these two particle orbitals do not stay constant in time, as shown in Figs. 8(b) and 8(c). This is quite different from the optical excitation case where the hole and particle orbitals remain largely constant in time. The particle orbital for the third transition mode has negligible contribution from the eigenstates localized on the benzene molecule, and it is mostly of the projectile proton [see Fig. 9(f)]. To summarize, three dominant transition modes are observed in the electronic stopping process: (1) HOMO → LUMO transition and HOMO → proton charge transfer, (2) HOMO-1 → LUMO+1 transition and HOMO-1 → proton charge transfer, and (3) HOMO-4 → proton charge transfer.

FIG. 8.

(a) Total hole population loss (red dashed line) and time-dependent particle population [i.e., σnp(t)bn(t)2] of the three most dominant transition modes (three solid lines). [(b)–(d)] The relative magnitudes of the projection matrix onto the KS eigenstates are shown using the color map, ranging from 0% to 100%. The x-axis is time, and the y-axis is the KS (occupied/unoccupied) eigenstates index.

FIG. 8.

(a) Total hole population loss (red dashed line) and time-dependent particle population [i.e., σnp(t)bn(t)2] of the three most dominant transition modes (three solid lines). [(b)–(d)] The relative magnitudes of the projection matrix onto the KS eigenstates are shown using the color map, ranging from 0% to 100%. The x-axis is time, and the y-axis is the KS (occupied/unoccupied) eigenstates index.

Close modal
FIG. 9.

Particle and hole orbitals for the three most dominant transition modes at t = 68 a.u., after the projectile proton has passed by the benzene molecule.

FIG. 9.

Particle and hole orbitals for the three most dominant transition modes at t = 68 a.u., after the projectile proton has passed by the benzene molecule.

Close modal
TABLE III.

Time-averaged projection of time-dependent particle/hole orbitals onto a Kohn–Sham eigenstate for the three dominant transition modes, Di,n(t) and Da,n(t), after the projectile proton moves past the benzene molecule (after t = 40 a.u.). i refers to the occupied Kohn–Sham eigenstates. a refers to the unoccupied (virtual) Kohn–Sham eigenstates. NA indicates that no major Kohn–Sham eigenstates contribute to the particle orbital. Var(Di,n(t)) is the variance for Di,n(t).

ModeiDi,n(t)%Var(Di,n(t))%aDa,n(t)%
HOMO 99.4 0.00 LUMO 11.3 
HOMO-1 98.4 0.03 LUMO+1 56.3 
HOMO-4 90.9 0.06 NA NA 
ModeiDi,n(t)%Var(Di,n(t))%aDa,n(t)%
HOMO 99.4 0.00 LUMO 11.3 
HOMO-1 98.4 0.03 LUMO+1 56.3 
HOMO-4 90.9 0.06 NA NA 

We formulated dynamical transition orbitals as a unitary transformation of time-dependent Kohn–Sham orbitals in real-time time-dependent density functional theory (RT-TDDFT). Each dynamical transition orbital consists of a single pair of time-dependent particle and hole orbitals, and the orbital populations yield time-dependent transitions of the particle–hole pair. Therefore, electron dynamics can be interpreted in terms of a set of single-particle excitations from a hole orbital to a particle orbital within individual dynamical transition orbitals. The gauge invariant nature of the dynamical transition orbital formulation makes the approach particularly appealing to RT-TDDFT dynamics. The equations of motion for the particle and hole orbitals are derived to show that the particle/hole orbitals and populations evolve continuously in time. Within the linear response approximation, we show how the particle–hole transitions within individual dynamical transition orbitals reduce to the transition between the well-known static natural transition orbitals.13 We demonstrated the utility of this new dynamical transition orbital formulation through its application to the RT-TDDFT simulation of isolated molecules and extend systems. In particular, we studied the optical excitation of a water molecule and a benzene molecule and also of crystalline silicon with the periodic boundary conditions. Then, the approach was applied to a difficult case of electronic stopping excitation in a benzene molecule with a proton projectile. The dynamical transition orbitals enabled clear interpretation of the excitation dynamics in terms of the particle–hole transitions within the benzene molecule and the charge transfer between the benzene molecule and the projectile proton.

See the supplementary material for the figures for the optical pulse electric field with the particle orbital occupation changes.

This work was supported by the National Science Foundation under Award Nos. CHE-1954894 and OAC-17402204. An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program. This research used resources of the Argonne Leadership Computing Facility (ALCF), which is a DOE Office of Science User Facility supported under Contract No. DE-AC02-06CH11357.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1. Brief review of mathematical properties

We here briefly review necessary mathematics used to discuss the dynamical transition orbital formulation in Sec. II.

a. Projection operator

Projection operators for the occupied space (P^o) and unoccupied (virtual) space (P^v) in the ground state are defined as

P^o=iocc.|ψiψi|,
(A1)
P^v=avir.|ψaψa|.
(A2)

Here, indexes i and a are used for the occupied and unoccupied KS eigenstates, respectively. For any wave function ϕ in the Hilbert space, there exist unique vectors v and w such that

|ϕ=ivi|ψi+awa|ψa,
(A3)

where the element of the vectors v and w can be determined simply by

|ψiψi|ϕ=vi|ψi,
(A4)
|ψaψa|ϕ=wa|ψa.
(A5)

Operating the P^o+P^v on ϕ, one finds

P^o+P^v|ϕ=P^o+P^vivi|ψi+awa|ψa=iviP^o|ψi+awaP^v|ψa=|ϕ.
(A6)

Therefore, the projection operators for the occupied space (P^o) and virtual space (P^v) in the ground state are complementary,

P^o+P^v=Î,
(A7)

where Î is the identity operator. These projection operators P^o and P^v are Hermitian (i.e., P^o=P^o), and they are also idempotent (i.e., P^oP^o=P^o).

b. Singular value decomposition

According to the singular value decomposition (SVD), a general operator Ô, which may be any real or complex matrix, can be expressed as

Ô=n|AnλnBn|,
(A8)

where λn is real and positive and An, Bn are orthonormal vectors that can be obtained from the eigen-equations,

ÔÔ|Bn=λn2|Bn
(A9)

and

ÔÔ|An=λn2|An.
(A10)

Here, the singular value is required to be real and positive,

λn=λn2.

Meanwhile, there might be infinite choices for phases of Bn and An from Eqs. (A9) and (A10). In practice, only one eigen-equation needs to be solved, whereas another corresponding vector is solved by projection,

Ô|Bm=n|AnλnBn|Bm=n|Anλnδm,n=λm|Am.
(A11)

2. Natural transition orbitals in linear response TDDFT

Within the Tamm–Dancoff approximation (TDA), the coupling of the excitation to the de-excitation is neglected. The TDA is widely used to solve the Casida equation3,5 in linear-response TDDFT (LR-TDDFT), and the transition density is given by59 

Ψ0|ρ^(r)|Ψn=iaXainψi*(r)ψa(r),
(A12)

where ψa and ψi are the unoccupied and occupied Kohn–Sham eigenstates, respectively. The tensor Xn (with the dimension of Nvir × Nocc) are obtained from solving the Casida equation. Note that this tensor obeys the orthogonality property,

aiXaim*Xain=δm,n.
(A13)

By defining an operator

X^n=ai|ψaXainψi|,
(A14)

the singular value decomposition (SVD) yields

X^n=i|ϕihλinϕip|,
(A15)

where λin are the singular values. Due to the normalization condition, the sum of the singular value squared is 1, i.e., iλi2=1. In most cases, for each n, a single pair of particle and hole orbitals can be identified (from {ϕih,ϕip}) for having the largest singular value close to 1, whereas other singular values are close to 0. Then, we can approximate that the single pair of particle and hole orbitals (a particle–hole excitation indexed with n) dominates the many-body excitation60,61 such that

Xn|ϕnhϕnp|,
(A16)

where ϕnp(r) and ϕnh(r) are the particle and hole orbitals that are identified as having the singular value close to 1. Then, this single particle–hole excitation approximation allows us to write Eq. (A12) simply as a product of a pair of particle and hole orbitals,

Ψ0|ρ^(r)|Ψnϕnh*(r)ϕnp(r).
(A17)
1.
E.
Runge
and
E. K. U.
Gross
, “
Density-functional theory for time-dependent systems
,”
Phys. Rev. Lett.
52
(
12
),
997
(
1984
).
2.
C. A.
Ullrich
,
Time-Dependent Density-Functional Theory: Concepts and Applications
(
Oxford University Press
,
Oxford
,
2011
), pp.
124
132
.
3.
M. E.
Casida
, “
Time-dependent density functional response theory for molecules
,”
Recent Advances In Density Functional Methods: (Part I)
(World Scientific,
1995
) pp.
155
192
.
4.
A.
Dreuw
and
M.
Head-Gordon
, “
Single-reference ab initio methods for the calculation of excited states of large molecules
,”
Chem. Rev.
105
(
11
),
4009
4037
(
2005
).
5.
M. E.
Casida
, “
Time-dependent density-functional theory for molecules and molecular solids
,”
J. Mol. Struct.: THEOCHEM
914
(
1-3
),
3
18
(
2009
).
6.
Y.
Miyamoto
,
A.
Rubio
, and
D.
Tománek
, “
Real-time ab initio simulations of excited carrier dynamics in carbon nanotubes
,”
Phys. Rev. Lett.
97
(
12
),
126104
(
2006
).
7.
Y.
Tateyama
 et al, “
Real-time propagation time-dependent density functional theory study on the ring-opening transformation of the photoexcited crystalline benzene
,”
J. Chem. Phys.
124
(
12
),
124507
(
2006
).
8.
X.
Li
and
J. C.
Tully
, “
Ab initio time-resolved density functional theory for lifetimes of excited adsorbate states at metal surfaces
,”
Chem. Phys. Lett.
439
(
1-3
),
199
203
(
2007
).
9.
M. E.
Casida
, “
Time-dependent density functional response theory for molecules
,” in
Recent Advances In Density Functional Methods: (Part I)
(World Scientific,
1995
), pp.
155
192
.
10.
X.
Li
 et al, “
Real-time time-dependent electronic structure theory
,”
Chem. Rev.
120
,
9951
(
2020
).
11.
F.
Furche
, “
On the density matrix based approach to time-dependent density functional response theory
,”
J. Chem. Phys.
114
(
14
),
5982
5992
(
2001
).
12.
Y.
Li
and
C. A.
Ullrich
, “
The particle-hole map: Formal derivation and numerical implementation
,”
J. Chem. Phys.
145
(
16
),
164107
(
2016
).
13.
R. L.
Martin
, “
Natural transition orbitals
,”
J. Chem. Phys.
118
(
11
),
4775
4777
(
2003
).
14.
C.
Adamo
and
D.
Jacquemin
, “
The calculations of excited-state properties with time-dependent density functional theory
,”
Chem. Soc. Rev.
42
(
3
),
845
856
(
2013
).
15.
K.
Yabana
and
G. F.
Bertsch
, “
Time-dependent local-density approximation in real time
,”
Phys. Rev. B
54
(
7
),
4484
(
1996
).
16.
K.
Lopata
and
N.
Govind
, “
Modeling fast electron dynamics with real-time time-dependent density functional theory: Application to small molecules and chromophores
,”
J. Chem. Theory Comput.
7
(
5
),
1344
1355
(
2011
).
17.
K.
Yabana
and
G. F.
Bertsch
, “
Time-dependent local-density approximation in real time: Application to conjugated molecules
,”
Int. J. Quantum Chem.
75
(
1
),
55
66
(
1999
).
18.
M. A. L.
Marques
 et al, “
Time-dependent density-functional approach for biological chromophores: The case of the green fluorescent protein
,”
Phys. Rev. Lett.
90
(
25
),
258101
(
2003
).
19.
F.
Ding
 et al, “
Quantum coherent plasmon in silver nanowires: A real-time TDDFT study
,”
J. Chem. Phys.
140
(
24
),
244705
(
2014
).
20.
D. S.
Ravithree
 et al, “
Real-time TDDFT investigation of optical absorption in gold nanowires
,”
J. Phys. Chem. C
123
(
23
),
14734
14745
(
2019
).
21.
K.
Lopata
 et al, “
Linear-response and real-time time-dependent density functional theory studies of core-level near-edge x-ray absorption
,”
J. Chem. Theory Comput.
8
(
9
),
3284
3292
(
2012
).
22.
C. D.
Pemmaraju
 et al, “
Velocity-gauge real-time TDDFT within a numerical atomic orbital basis set
,”
Comput. Phys. Commun.
226
,
30
38
(
2018
).
23.
A. R.
Attar
 et al, “
Femtosecond x-ray spectroscopy of an electrocyclic ring-opening reaction
,”
Science
356
(
6333
),
54
59
(
2017
).
24.
C. D.
Pemmaraju
, “
Valence and core excitons in solids from velocity-gauge real-time TDDFT with range-separated hybrid functionals: An LCAO approach
,”
Comput. Condens. Matter
18
,
e00348
(
2019
).
25.
A.
Ojanperä
,
A. V.
Krasheninnikov
, and
M.
Puska
, “
Electronic stopping power from first-principles calculations with account for core electron excitations and projectile ionization
,”
Phys. Rev. B
89
(
3
),
035120
(
2014
).
26.
R.
Ullah
,
E.
Artacho
, and
A. A.
Correa
, “
Core electrons in the electronic stopping of heavy ions
,”
Phys. Rev. Lett.
121
(
11
),
116401
(
2018
).
27.
J. M.
Pruneda
 et al, “
Electronic stopping power in LiF from first principles
,”
Phys. Rev. Lett.
99
(
23
),
235501
(
2007
).
28.
I.
Maliyov
,
J.-P.
Crocombette
, and
F.
Bruneval
, “
Electronic stopping power from time-dependent density-functional theory in Gaussian basis
,”
Eur. Phys. J. B
91
(
8
),
172
(
2018
).
29.
A.
Schleife
,
Y.
Kanai
, and
A. A.
Correa
, “
Accurate atomistic first-principles calculations of electronic stopping
,”
Phys. Rev. B
91
(
1
),
014306
(
2015
).
30.
K. G.
Reeves
,
Y.
Yao
, and
Y.
Kanai
, “
Electronic stopping power in liquid water for protons and α particles from first principles
,”
Phys. Rev. B
94
(
4
),
041108
(
2016
).
31.
D. C.
Yost
and
Y.
Kanai
, “
Electronic stopping for protons and α particles from first-principles electron dynamics: The case of silicon carbide
,”
Phys. Rev. B
94
(
11
),
115107
(
2016
).
32.
D. C.
Yost
,
Y.
Yao
, and
Y.
Kanai
, “
Examining real-time time-dependent density functional theory nonequilibrium simulations for the calculation of electronic stopping power
,”
Phys. Rev. B
96
(
11
),
115134
(
2017
).
33.
D. C.
Yost
,
Y.
Yao
, and
Y.
Kanai
, “
First-principles modeling of electronic stopping in complex matter under ion irradiation
,”
J. Phys. Chem. Lett.
11
(
1
),
229
237
(
2020
).
34.
D. C.
Yost
and
Y.
Kanai
, “
Electronic excitation dynamics in DNA under proton and α-particle irradiation
,”
J. Am. Chem. Soc.
141
(
13
),
5241
5251
(
2019
).
35.
Y.
Yao
,
D. C.
Yost
, and
Y.
Kanai
, “
K-shell core-electron excitations in electronic stopping of protons in water from first principles
,”
Phys. Rev. Lett.
123
(
6
),
066401
(
2019
).
36.
T.
Otobe
 et al, “
First-principles electron dynamics simulation for optical breakdown of dielectrics under an intense laser field
,”
Phys. Rev. B
77
(
16
),
165104
(
2008
).
37.
C.
Wang
 et al, “
First-principles calculations of the electron dynamics during femtosecond laser pulse train material interactions
,”
Phys. Lett. A
375
(
36
),
3200
3204
(
2011
).
38.
K.
Yabana
 et al, “
Time-dependent density functional theory for strong electromagnetic fields in crystalline solids
,”
Phys. Rev. B
85
(
4
),
045134
(
2012
).
39.
A.
Yamada
and
K.
Yabana
, “
Energy transfer from intense laser pulse to dielectrics in time-dependent density functional theory
,”
Eur. Phys. J. D
73
(
5
),
87
(
2019
).
40.
Y.
Li
and
C. A.
Ullrich
, “
Time-dependent transition density matrix
,”
Chem. Phys.
391
(
1
),
157
163
(
2011
).
41.
G. F.
Bertsch
 et al, “
Real-space, real-time method for the dielectric function
,”
Phys. Rev. B
62
(
12
),
7998
(
2000
).
42.
W.
Jia
 et al, “
Fast real-time time-dependent density functional theory calculations with the parallel transport gauge
,”
J. Chem. Theory Comput.
14
(
11
),
5645
5652
(
2018
).
43.
W.
Jia
and
L.
Lin
, “
Fast real-time time-dependent hybrid functional calculations with the parallel transport gauge and the adaptively compressed exchange formulation
,”
Comput. Phys. Commun.
240
,
21
29
(
2019
).
44.
A.
Dong
and
L.
Lin
, “
Quantum dynamics with the parallel transport gauge
,”
Multiscale Model. Simul.
18
(
2
),
612
645
(
2020
).
45.
D. C.
Yost
,
Y.
Yao
, and
Y.
Kanai
, “
Propagation of maximally localized Wannier functions in real-time TDDFT
,”
J. Chem. Phys.
150
(
19
),
194113
(
2019
).
46.
K.
Pernal
,
O.
Gritsenko
, and
E.
Jan Baerends
, “
Time-dependent density-matrix-functional theory
,”
Phys. Rev. A
75
(
1
),
012506
(
2007
).
47.
K.
Pernal
 et al, “
Adiabatic approximation of time-dependent density matrix functional response theory
,”
J. Chem. Phys.
127
(
21
),
214101
(
2007
).
48.
K. J. H.
Giesbertz
,
E. J.
Baerends
, and
O. V.
Gritsenko
, “
Charge transfer, double and bond-breaking excitations with time-dependent density matrix functional theory
,”
Phys. Rev. Lett.
101
(
3
),
033004
(
2008
).
49.
H.
Appel
and
E. K. U.
Gross
, “
Time-dependent natural orbitals and occupation numbers
,”
Europhys. Lett.
92
(
2
),
23001
(
2010
).
50.
R.
Van Meer
,
O. V.
Gritsenko
, and
E. J.
Baerends
, “
Natural excitation orbitals from linear response theories: Time-dependent density functional theory, time-dependent Hartree-Fock, and time-dependent natural orbital functional theory
,”
J. Chem. Phys.
146
(
4
),
044119
(
2017
).
51.
A.
Schleife
 et al, “
Plane-wave pseudopotential implementation of explicit integrators for time-dependent Kohn-Sham equations in large-scale simulations
,”
J. Chem. Phys.
137
(
22
),
22A546
(
2012
).
52.
A.
Schleife
 et al, “
Quantum dynamics simulation of electrons in materials on high-performance computers
,”
Comput. Sci. Eng.
16
(
05
),
54
60
(
2014
).
53.
F.
Gygi
, “
Architecture of Qbox: A scalable first-principles molecular dynamics code
,”
IBM J. Res. Dev.
52
(
1.2
),
137
144
(
2008
).
54.
E. W.
Draeger
 et al, “
Massively parallel first-principles simulation of electron dynamics in materials
,”
J. Parallel Distrib. Comput.
106
,
205
214
(
2017
).
55.
D. R.
Hamann
,
M.
Schlüter
, and
C.
Chiang
, “
Norm-conserving pseudopotentials
,”
Phys. Rev. Lett.
43
(
20
),
1494
(
1979
).
56.
D.
Vanderbilt
, “
Optimally smooth norm-conserving pseudopotentials
,”
Phys. Rev. B
32
(
12
),
8412
(
1985
).
57.
T. P.
Rossi
 et al, “
Kohn–Sham decomposition in real-time time-dependent density-functional theory: An efficient tool for analyzing plasmonic excitations
,”
J. Chem. Theory Comput.
13
(
10
),
4779
4790
(
2017
).
58.
A.
Castro
,
M. A. L.
Marques
, and
A.
Rubio
, “
Propagators for the time-dependent Kohn–Sham equations
,”
J. Chem. Phys.
121
(
8
),
3425
3433
(
2004
).
59.
T.
Etienne
, “
Transition matrices and orbitals from reduced density matrix theory
,”
J. Chem. Phys.
142
(
24
),
244103
(
2015
).
60.
P. R.
Surján
, “
Natural orbitals in CIS and singular-value decomposition
,”
Chem. Phys. Lett.
439
(
4-6
),
393
394
(
2007
).
61.
I.
Mayer
, “
Using singular value decomposition for a compact presentation and improved interpretation of the CIS wave functions
,”
Chem. Phys. Lett.
437
(
4-6
),
284
286
(
2007
).

Supplementary Material