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.

## I. INTRODUCTION

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 Furche^{11} 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 electrons^{15–20} and core electrons^{21–24} as well as other types of non-equilibrium dynamic phenomena, such as electronic stopping^{25–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.

## II. DYNAMICAL TRANSITION ORBITALS

### A. Formulation

We introduce dynamical transition orbitals, $\psi 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,

where $\psi nh(r,t)$ and $\psi 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

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

where {*ϕ*_{n}(**r**, *t*)} are the TD-KS wavefunctions. Physically observable $\u27e8\xc2(t)\u27e9=Tr[\xc2\Gamma ^(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 *N*_{elec}/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,

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

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

Similar to the static natural transition orbitals,^{13} we obtain the dynamical transition orbitals and the hole orbitals via the SVD of the $\u015c$ operator,

where *a*_{n}(*t*) are the singular values. Using the SVD property Eq. (A9) in the Appendix and Eq. (6), we have

#### 2. Particle orbitals

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

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

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

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 $|\varphi nDT(t)\u2009$ notation instead of $|\psi nDT(t)\u2009$. Analogous to the hole orbital case [Eqs. (8) and (9)], we have

and

#### 3. Equivalence between $|\varphi nDT(t)\u2009$ and $|\psi nDT(t)\u2009$

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

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

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

Comparing this equation to Eq. (13), one recognizes that both $|\varphi nDT(t)\u2009$ and $|\psi nDT(t)\u2009$ are the same set of eigenvectors (up to an arbitrary phase) for the operator Γ(*t*)$Pv$Γ(*t*). The corresponding eigenvalues equate *b*_{n}(*t*)^{2} and 1 − *a*_{n}(*t*)^{2}, and Eq. (2) is proved. Note that *a*_{n}(*t*) and *b*_{n}(*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

The left-hand side can be expanded as

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

Hence, Eq. (18) becomes

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

### B. Density matrix in particle–hole orbital basis

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

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

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

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 2*N*, 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

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

### C. Equation of motion

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

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

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

Multiplying both sides by $\u2009\psi mh(t)|$, we have

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

Equation (30) can be written as

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,

Similarly, for the particle population, we have

When *m* ≠ *n*, because of the orthogonality, $\sigma \u0307nh(t)\u2009\psi mh(t)|\psi nh(t)\u2009$ vanishes in Eq. (30), and the time derivative of the hole orbital can be written as an expansion of other occupied hole orbitals,

Similarly, for the particle orbitals, we have

where the summation index *m*′ is used to distinguish from m. While *N*_{occ} is the number of occupied orbitals, $Nvir$ is the number of unoccupied (virtual) orbitals, which is finite in practical calculations. Although the $\psi 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 $\sigma nh/p(t)=\sigma 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., $|\u2009\varphi a|\psi np(t)\u2009|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).

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., $\sigma nh(t)=an(t)2$], yield the expansion coefficients with 0 ≤ *a*_{n}(*t*), *b*_{n}(*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., $\sigma nh(t)+\sigma 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.

## III. LINEAR RESPONSE APPROXIMATION

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 *b*_{n}(*t*)^{2} ≈ 0. Due to Eq. (2), we have 0 ≤ *a*_{n}(*t*) ≤ *a*_{n}(*t*)^{2} with *a*_{n}(*t*) ≈ 1. Under this approximation of the weak perturbation, the electron density [Eq. (26)] and the one-particle density operator [Eq. (25)] yield

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,

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,

We now connect this dynamical transition orbital expression [Eq. (39)] to the many-body expression and show how the particle occupation number *b*_{n}(*ω*) 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

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

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

where $\Psi n$ is the many-body state in the n-th excited state and the operator is $\rho ^(r)=\u2211iNelec.\delta (r\u2212ri\u2032)$.^{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

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

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

Therefore, the poles in *b*_{n}(*ω*) indicate the system’s excitation energy *E*_{n} − *E*_{0} of the many-body expression. This perspective on the particle occupation number is demonstrated in the following in a numerical example.

## IV. APPLICATION

### A. Numerical implementation

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 $\psi i$ are calculated, while the unoccupied KS eigenstates $\psi a$ are not necessary.

$\varphi n(t)$ or $\Gamma ^(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 functions

^{45}and the parallel transport gauge orbitals.^{42}At each time step,

*a*_{n}(*t*) and $|\psi nDT(t)\u2009$ are obtained by solving the eigenvalue problem

$|\psi nh(t)\u2009$ are obtained by acting the projection operator for the occupied subspace,

The particle orbitals are obtained via

where *b*_{n}(*t*) can be obtained from *b*_{n}(*t*)^{2} = 1 − *a*_{n}(*t*)^{2} [Eq. (2)].In some cases, particularly when *b*_{n}(*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,

where *N*_{eff} defines the smaller subspace spanned by unoccupied KS eigenstates $\psi 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

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,

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}

### B. Application: Optical excitation

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,

where *ω*_{0} is the specific frequency the excitation is centered on. *γ* controls the Gaussian width of the external potential applied. The amplitude *E*_{0} 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 bohrs^{3} 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, *t*_{0} = 75 a.u., *γ* = 2(*eV*/*ℏ*)^{2}, and *E*_{0} = 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 *D*_{i,n}(*t*) and *D*_{a,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.

Molecule . | Mode . | i . | ⟨D_{i,n}(t)⟩%
. | a . | ⟨D_{a,n}(t)⟩%
. |
---|---|---|---|---|---|

H_{2}O | 1 | HOMO | 99.8 | LUMO | 94.4 |

C_{6}H_{6} | 1 | HOMO | 97.3 | LUMO+1 | 96.8 |

2 | HOMO-1 | 99.1 | LUMO | 99.3 |

Molecule . | Mode . | i . | ⟨D_{i,n}(t)⟩%
. | a . | ⟨D_{a,n}(t)⟩%
. |
---|---|---|---|---|---|

H_{2}O | 1 | HOMO | 99.8 | LUMO | 94.4 |

C_{6}H_{6} | 1 | HOMO | 97.3 | LUMO+1 | 96.8 |

2 | HOMO-1 | 99.1 | LUMO | 99.3 |

#### 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 bohrs^{3}) 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, *t*_{0} = 200 a.u., *γ* = 1(eV/*ℏ*)^{2}, and *E*_{0} = 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.

### C. Application: Electronic stopping

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) integrator^{58} 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 bohrs^{3} 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., $\u2212\u2211\sigma 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 $\sigma 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.

Mode . | i . | $\u27e8Di,n(t)\u27e9$% . | Var(D_{i,n}(t))%
. | a . | $\u27e8Da,n(t)\u27e9$% . |
---|---|---|---|---|---|

1 | HOMO | 99.4 | 0.00 | LUMO | 11.3 |

2 | HOMO-1 | 98.4 | 0.03 | LUMO+1 | 56.3 |

3 | HOMO-4 | 90.9 | 0.06 | NA | NA |

Mode . | i . | $\u27e8Di,n(t)\u27e9$% . | Var(D_{i,n}(t))%
. | a . | $\u27e8Da,n(t)\u27e9$% . |
---|---|---|---|---|---|

1 | HOMO | 99.4 | 0.00 | LUMO | 11.3 |

2 | HOMO-1 | 98.4 | 0.03 | LUMO+1 | 56.3 |

3 | HOMO-4 | 90.9 | 0.06 | NA | NA |

## V. CONCLUSION

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.

## SUPPLEMENTARY MATERIAL

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

## ACKNOWLEDGMENTS

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.

## DATA AVAILABILITY

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

### APPENDIX: SINGULAR VALUE DECOMPOSITION AND NATURAL TRANSITION ORBITALS

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

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

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

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

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

where $\xce$ is the identity operator. These projection operators $P^o$ and $P^v$ are Hermitian (i.e., $P^o\u2020=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 $\xd4$, which may be any real or complex matrix, can be expressed as

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

and

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

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,

#### 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 equation^{3,5} in linear-response TDDFT (LR-TDDFT), and the transition density is given by^{59}

where *ψ*_{a} and *ψ*_{i} are the unoccupied and occupied Kohn–Sham eigenstates, respectively. The tensor *X*^{n} (with the dimension of *N*_{vir} × *N*_{occ}) are obtained from solving the Casida equation. Note that this tensor obeys the orthogonality property,

By defining an operator

the singular value decomposition (SVD) yields

where $\lambda in$ are the singular values. Due to the normalization condition, the sum of the singular value squared is 1, i.e., $\u2211i\lambda i2=1$. In most cases, for each *n*, a single pair of particle and hole orbitals can be identified (from ${\varphi ih,\varphi 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 excitation^{60,61} such that

where $\varphi np(r)$ and $\varphi 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,