Sequential double excitations from linear-response time-dependent density functional theory

Molecules and solids can absorb and emit not only one but multiple photons over a short time period, leading to the ability to explore the spectroscopic characteristics of excited states. Such multiphoton excitations/emissions can either be coherent or incoherent. But, even if they are incoherent (as considered here) it is still a significant challenge to use electronic structure theory to model experimental observables from many frontier non-linear optical spectroscopic measurements. Wave-function theory and its computational tools could provide us with all the information we need to understand and correlate multi-photon processes; however, computational demands tend to be high. New technological advancements and algorithms simplifying the computation of the tensor elements required in high-level calculations may eventually enable correlated-wave function methods to be applicable for large molecular systems, but this is not available now.

In parallel with theoretical developments, experiments nowadays can resolve time scales where coherent nuclear motions in molecules can be tracked using methods like transient UV/vis or X-ray absorption, or Raman-based methods. Careful monitoring of the states of molecules is crucial in any mechanism proposed to control the course of photochemical reactions. Moreover, the fast-paced increase in experimental techniques and devices motivates the quest for theoretical and computational methods that offer scientists the opportunity to simulate ultrafast spectroscopic and quantum control experiments.

Methods based on time-dependent (TD) density-functional theory (TDDFT) are often useful for calculating dynamical properties of molecules in a computationally affordable manner. The standard type of calculations within TDDFT, used by solid-state scientists, is a natural extension of Kohn-Sham (KS) theory, in which the electrons are subject to a local TD potential. In the molecular case, however, a mixture of TD Hartree-Fock (non-local) and KS potentials is employed for the computations. This mixture can be formulated using the generalized KS method,¹ in which the electrons interact through an auxiliary non-local potential. The use of the non-local potential produces more accurate HOMO-LUMO gaps and improved excitation energies.

The quadratic response TDDFT formalism^2,3 provides a robust method to describe the absorption/emission of two photons (sequentially and/or coherently). The implementation of this theory by Sałek et al.⁴ comprises calculation of amplitudes corresponding to excitations of occupied orbitals, determination of the tensor elements of the first- and second-order exchange-correlation kernels, evaluation of various commutators of the amplitude and XC kernels, iterative solution of the resulting eigenvalue problem to find the excitation energies, and computation of the response functions. This quadratic response method scales as N⁶, a considerable expense for applications to large molecular structures. If incoherent (sequential) double excitations are required, the scaling can be reduced to that of linear-response calculations, as we show in this work.

Alternative ideas are being proposed. Fischer, Cramer, and Govind⁵ developed an algorithm in which an excited initial reference is produced, then this new reference is propagated in real time and a new spectrum is sampled. The authors noted that the excitation frequencies obtained from the real time propagation are shifted with respect to those calculated from the original spectrum (determined by a regular linear-response TDDFT calculation). This peak-shifting problem is due to the lack of initial-state and memory dependencies of traditional XC kernels.⁶ 

In this manuscript, we present a scheme to calculate the spectrum starting from an excited state of a molecule or solid. The methodology proceeds as follows: (i) One solves the regular linear response (LR) problem and selects a target excited state. Using the excitation vectors of the target excited state, the occupied ground-state orbitals are perturbed by the virtual ones (and vice versa). The intensity of the perturbation is controlled by a parameter, denoted here as λ. (ii) We show that the residues of the polarizability tensor of the perturbed system give the transition dipole moments and the oscillator strengths corresponding to transitions from the target excited state. The estimation of the oscillator strengths requires resolving the Casida equations with reference to the perturbed orbitals and calculation of the excited-state transition dipoles. These are obtained as the difference between the transition dipoles of the perturbed and unperturbed systems divided by the strength of the perturbation. Finally, (iii) a formal property, in the limit λ → 0 the excitation frequencies from the target excited state tend to those derived from the traditional LR spectrum. Under this property, we consider small values of the perturbation parameter and use the excitation energies obtained from the original LR problem to calculate the frequencies for excitations from the target excited state. The total computational cost to calculate all the oscillator strengths in a given energy window is a multiple (as low as two) of the cost of a regular LR routine. To validate the theory, we test the method by comparing with experimental data for X-ray absorption near-edge spectroscopy (XANES) and transient optical (Vis/UV) absorption.

Let us assign $Hˆ0$ as the electronic Hamiltonian of the molecular system (within the Born-Oppenheimer picture) and assume that the spectrum and eigenstates of the system are known. Denoting the eigenkets of the system as {|Ψ_J〉}, we have $Hˆ|ΨJ〉=EJ|ΨJ〉$⁠; we assume the wavefunctions {Ψ_J} are real-valued. The oscillator strength for a transition from the ground state to a state J is $fJ=2/3ΩJ|〈ΨJ|μˆ|Ψ0〉|2$⁠, where Ω_J = E_J − E₀ is the excitation frequency, |Ψ₀〉 is the ground-state ket, and $μˆ=(xˆ,yˆ,zˆ)$ is the dipole vector operator. Suppose we choose the state I as the target excited state, the transition strength of the J←I excitation is $fJI=2/3ΩJI|〈ΨJ|μˆ|ΨI〉|2$⁠, where Ω_JI = Ω_J − Ω_I.

Now suppose that the ground state is slightly perturbed by the target excited state I and that the system is subject to an electric field, $Vˆ(t)=−μˆ⋅E(t)$⁠. The TD Hamiltonian of the molecule is $Hˆ(t)=Hˆ0+Vˆ(t)$⁠. The problem we need to solve now is $i∂t|Ψ(t;λ)〉=Hˆ(t)|Ψ(t;λ)〉$⁠, with initial condition Ψ(0) = Ψ₀ + λΨ_I, where λ is a real number. The TD wavefunction depends on the strength of the initial state perturbation; hence, we denote it as Ψ(t; λ).

We define the following polarizability tensor in frequency space as

We use l and m as indices for spatial degrees of freedom, (x, y, z). μ_l(t; λ) is the lth component of the net dipole, $〈Ψ(t;λ)|μˆ|Ψ(t;λ)〉$⁠. Henceforth we use the superindex “F” to denote quantities in frequency space. A more general form of the polarizability tensor (in real time) is δμ_l(t; λ)/δE_m(s) at E = 0, where s < t. For our purposes, however, Eq. (1) is sufficient to obtain the required information.

In  Appendix A we show that

Therefore the oscillator strength of the J←I transition can be extracted from the trace of the polarizability tensor, as follows:

The polarizability tensor is a second degree polynomial in terms of λ; the right hand side of the above equation is independent of λ. Alternatively, if we define the λ-dependent transition moment as

then

This form is more convenient than Eq. (3) because it only involves first derivatives. In Secs. III and IV, we elaborate the scheme to calculate the above quantity by perturbing the initial state of the system of auxiliary electrons.

This work is based on propagation of orbitals, which (for quantum chemical applications) are usually handled in computational packages in the form of molecular eigenvectors. Let us thus present a slightly different working formalism based on orbital propagation. We use the indices i, j to denote occupied orbitals and a, b for virtual ones. To refer to any orbital in general (occupied or virtual) we use the index p. The labels for the spin variables are σ and τ. Here we adopt a generalized Kohn-Sham (GKS) approach: The TD electronic density of the system is represented by an auxiliary system of electrons which interact via an orbital dependent, parametrized, electron-electron interaction. Let us denote this interaction operator as $wˆσθ$⁠. Consider, for example, the potential $wˆσθ=θ(jˆ−kˆσ)$⁠, where θ is a positive number, $jˆ=jˆ↑+jˆ↓$⁠, and $kˆσ$ are the Coulomb and σ-exchange operators, respectively. The GKS equations are of the form (with ħ = 1)

the 1-particle spin-polarized Hamiltonian, $hˆσ(t)$⁠, is $tˆ+vˆs,σ[nσ′](t)$⁠, where $vˆs,σ[nσ′](t)=wˆσθ(t)+vˆHXC,σθ[nσ′](t)+vˆ0$⁠, $tˆ$ is the kinetic energy operator, and $vˆ0$ is the static 1-body external potential of the electron. The operator $uˆ(t)$ is a driving field. The spin-polarized Hartree-exchange-correlation (HXC) potential $vˆHXC,σθ$ is a functional of the TD spin-densities, $nσ(r,t)=∑i|ψiσ(r,t)|2,σ=↑,↓$⁠. This potential accounts for 100% the correlation effects and (1 − θ) × 100% the exchange and Hartree-repulsion interactions.

Let us express the Hamiltonian as $hˆσ(t)=hˆ0,σ+Δˆσ(t)$⁠, where $Δˆσ(t)=hˆσ(t)−hˆ0,σ$⁠. The Hamiltonian $hˆ0,σ$ is determined by the ground-state density of the system (which we assume is given). The initial condition for each orbital is ψ_iσ(0) = φ_iσ. {φ_pσ} are eigenfunctions of the zero order Hamiltonian (⁠$hˆ0,σ$⁠). Hence, $hˆ0,σ|φpσ〉=ϵpσ|φpσ〉$⁠.

We write the perturbed ket |ψ_iσ(t)〉 as the sum of its zero order form, exp(−iϵ_iσt)|φ_iσ〉, and a perturbation, which is a linear combination of virtual orbitals,

This ket represents a perturbed orbital, not an orbital associated to a pure, excited state of auxiliary electrons (Section V). Now, we expand Q_aiσ(t) in terms of two response vectors, x(t) and y(t),

Here, c is a complex number. The perturbation operator is re-expressed as

where

The Hamiltonian $hˆσ$ can be expanded linearly in terms of the response coefficients ${Qaiσ,Qaiσ∗}$⁠. If all the terms are factorized according to their common c and c^∗ factors, and the functions multiplying these factors are each set to zero, we arrive at the TD LR equation,

where $vaiσ(t)=〈φaσ|vˆ(t)|φiσ〉$ and $vaiσ†(t)=〈φaσ|vˆ†(t)|φiσ〉$⁠. The matrix A depends on the orbital excitation energies and the hybrid XC kernel, while B only depends on the latter. The specific form of these matrices is shown in  Appendix B. The Casida equations can be derived from Eq. (11) by setting $vˆ(t)=0$⁠, x(t) = Xexp(−iωt), and y(t) = Yexp(+iωt), where X and Y are real vectors. The solution to the Casida equations is a set of eigenvectors and frequencies. We denote this set as {X^J, Y^J, Ω_J}, where Ω_J is the frequency of the excitation from the ground state to the state J.

The linear equations derived above are consistent with the traditional LR framework in frequency space if we write

The expansion for the virtual orbitals, {|ψ_aσ〉}, is obtained by swapping the labels i and a in Eq. (7). The following condition ensures orthogonality between perturbed virtual and occupied orbitals:

This relation implies that X_iaσ = − Y_aiσ and Y_iaσ = − X_aiσ. Furthermore, the eigenvectors {X^J, Y^J} satisfy the orthonormalization condition (X^I)^TX^J − (Y^I)^TY^J = δ^IJ. The superscript “T” refers to transpose of vectors.

We now suppose that the driving field is weak and let c = 1 (as usual in TD perturbation theory). The perturbed kets are written as |ψ_iσ(t)〉 = exp(−iϵ_iσt)|φ_iσ〉 + |ξ_iσ(t)〉. The response orbital ξ_iσ(t) is expanded in terms of the eigenvectors. First let |ξ_iσ(t)〉 = [x_aiσ(t) + y_aiσ(t)]exp(−iϵ_iσt)|φ_aσ〉, then the expansion to solve the TD LR equations is

We refer to the parameter λ_J as a response amplitude. The evolution equation for this amplitude is obtained by inserting the above expansion in Eq. (11). The result is (in Sec. IV we derive a similar equation for a more general case)

The response amplitudes are controlled by the external TD field. For example, assume that the external field is a perturbation of the form $uˆ(t)=2vˆ0cos(ωt)$⁠, then $vˆ(t)=vˆ0exp(−iωt)$⁠. By solving the above equation, with initial conditions λ_J(0) = 0, we get

Furthermore, by solving the evolution equation in frequency space for systems that are subject to TD electric fields, one can calculate, for instance, the residues of the polarizability tensor, and show that the oscillator strengths agree with those obtained by Casida,⁷ that is,

where $μaiσ=〈φaσ|μˆ|φiσ〉$⁠.

In this section we consider the propagation of orbitals describing auxiliary electrons in presence of weak scalar fields, where the initial state is perturbed, and show that the residues of the polarizability tensor (in frequency space) provide the necessary matrix elements to calculate f_JI. We assume that the kernel is memoryless and independent of the initial wavefunctions of the auxiliary and real systems of electrons (although the formalism can be adapted for more general kernels).

The perturbed orbitals are prepared following the results from Sec. III. Suppose a weak external field is turned on (at some t < 0), and that the only amplitude responding to such perturbation is λ_I. This quantity evolves until it reaches some value (and then the field is off at t = 0), suppose that value is a positive real number, λ. If no further driving field is applied then the orbitals describing the system would be {exp(−iϵ_iσt)|ζ_iσ(t; λ)〉} (for t > 0), where

We denote |ζ_0,pσ〉 = |ζ_pσ(0)〉.

Suppose a second weak driving field is turned on, we write the new set of TD orbitals to study the corresponding response of the system as {η_pσ(t; λ)}; the initial conditions are {η_pσ(0; λ) = ζ_0,pσ(λ)}. The equation describing the evolution of the perturbed system is

where $hˆσ′(t;λ)$ is a new GKS Hamiltonian, $hˆσ′(t;λ)=hˆσ(t;λ)+Δˆσ′(t;λ)+uˆ(t)$⁠, and $uˆ(t)=cvˆ(t)+c∗vˆ†(t)$ is the second driving field. $hˆσ(t;λ)$ is the Hamiltonian used in Sec. III, it reads $hˆσ(t;λ)=tˆ+vˆs,σ[nσ′I](t)$⁠, where

The TD interaction operator is $Δˆσ′(t;λ)=vˆs,σ[n̄σ′](t)−vˆs,σ[nσ′I](t)$⁠. The spin-densities for the system in presence of the second field are $n̄σ(t)=∑i|ηiσ(r,t)|2$⁠, σ = ↑, ↓.

Based on Sec. III, we expand the kets {|η_iσ(t; λ)〉} as follows:

here $Q̄aiσ(t;λ)$ is written as in Eq. (7) in terms of new, TD response vectors, $x̄(t;λ)$ and $ȳ(t;λ)$⁠,

The expansion shown in Eq. (21) ensures that |η_iσ(t; λ)〉 tends to exp(−iϵ_iσt)|ζ_iσ(t; λ)〉 as the intensity of the second field goes to zero.

The next step is the linearization of the Hamiltonian $hˆσ′$⁠. We expand both $hˆσ$ and $Δˆσ′$⁠. $hˆσ$ is expanded with respect to the ground state, $hˆσ(t;λ)=hˆ0,σ+Δˆσ(t;λ)$⁠, where $Δˆσ(t;λ)$ is expressed linearly in terms of $XaiσI$ and $YaiσI$⁠. And $Δˆσ′(t;λ)$ is expanded with respect to the perturbed state, which is determined by ${nσI(t)}$⁠. $Δˆσ′(t;λ)$ is thus approximated as a linear function with respect to $Q̄aiσ(t;λ)$⁠. For numerical convenience, we ignore the time-dependency of $(n↑I,n↓I)$ (set $nσI(t)=nσI(0)$⁠) in the expansion of $Δˆσ′(t;λ)$⁠.

To obtain the TD LR equations we apply the bra 〈ζ_0,aσ(λ)| to both sides of Eq. (19). After neglecting quadratic quantities ( Appendix B) and factorizing common terms, we get

where $v̄aiσ(t;λ)=〈ζ0,aσ|vˆ(t)|ζ0,iσ〉$⁠, $v̄aiσ†(t;λ)=〈ζ0,aσ|vˆ†(t)|ζ0,iσ〉$⁠, and the matrices coupling the orbital excitations are

To solve these equations, a new basis of eigenvectors is needed. If the ground-state orbitals, and the eigenvectors are real quantities, then the second Casida equations can be written as

The solution to the above eigenvalue problem is a set of eigenvectors and transition frequencies, ${X̄J(λ),ȲJ(λ),ΩJ}$⁠. In analogy with the first LR formalism, here the eigenvectors also satisfy the orthonormality rule $δJI=(X̄J)TX̄I−(ȲJ)TȲI$⁠.

The LR eigenvalue equations for the perturbed reference have a very similar structure to those for the ground-state reference. They differ in the matrix elements involving the perturbation $Δˆ′$⁠. Calculating the matrix elements only requires the perturbed orbitals {ζ_0,pσ(λ)} and the density $nσI(r,0;λ)$⁠. In  Appendix B we show the explicit form of the matrix elements for the second LR problem.

In the exact theory the new frequencies obtained from the solution of the second linear response problem are equal to those from the first LR problem.⁶ In the present case, the excitation frequencies from Eq. (25) vary as functions of λ. Nonetheless, in the limit λ → 0 these frequencies become those obtained from the first LR eigenvalue equations.

We now take c = 1. Using the new basis of excitation vectors we can expand the TD vectors as $x̄(t;λ)=∑Jλ̄J(t;λ)X̄J(λ)$ and $ȳ(t;λ)=∑Jλ̄J∗(t;λ)ȲJ(λ)$⁠. Let us examine the response of the system to an additional perturbation, a TD electric field, $vˆ(t)=−E(t)⋅μˆ$⁠. By insertion of the expansions of $x̄$ and $ȳ$ in Eq. (23) we obtain

where K is an index running over excited states. Multiply the first and second rows by the transposed vectors $(X̄J)T$ and $(ȲJ)T$⁠, respectively, then complex-conjugate the second row, and add the two rows together. By means of the orthonormality rule, these operations yield

If we write $λ̄J(t)=∫0+∞dω/2πλ̄JF(ω)e−iωt$ then $λ̄JF(ω)=[X̄J+ȲJ]TūF(ω)/(ω−ΩJ)$⁠, where $ūaiσF(ω)=〈ζ0,aσ|uˆF(ω)|ζ0,iσ〉$⁠.

After some algebraic operations ( Appendix A), it can be shown that the polarizability obeys

where $μ̄l(λ)$ is the lth (x, y, or z) component of the single-particle transition dipole,

Eq. (28) can be re-expressed as

The second term in the right hand side of the above equation is $2ΩJI〈ΨJ|μˆl|Ψ0〉〈ΨJ|μˆm|Ψ0+λΨI〉$ ( Appendix A) and (in principle) is linear with respect to λ. The third term is negligible in comparison with the other two, or zero in the Tamm-Dancoff approximation (the case where $ȲJ=0$⁠, for all J).

The counterpart of Eq. (30) in wavefunction theory is ( Appendix A)

From the last two equations we obtain that

This function allows us to calculate f_JI, Eq. (5). Furthermore, the calculation of d_JI,l does not require modification to the conventional LR TDDFT algorithm. It is sufficient to generate the perturbed orbitals {ζ_0,pσ} and rerun the LR calculations.

The calculation of f_JI can be carried out applying perturbation theory (i.e., by solving the second eigenvalue problem perturbatively, inserting the explicit expansion of $X̄J$⁠, $ȲJ$⁠, and $μ̄$ in the definition of d_JI,l, and applying the first derivative) or using finite differences. In this work we choose the latter. To illustrate the application of our method, Section VI, we employ the formula,

where λ is a small number, and μ_J0,l is the usual transition dipole, [X^J + Y^J]^Tμ_l.

In our preliminary calculations we noted that the quality of the eigenvalues of Eq. (25) degrades for large values of λ. We thus use small values of λ to calculate f_JI. For λ < 1 there are small shifts in the second frequencies. We discard these small shifts and take the frequencies of the J←I transitions from the first eigenvalue problem (λ = 0) and calculate Ω_JI as Ω_J − Ω_I. In the formalism we assumed that the order of the spectrum after solving the second eigenvalue problem is the same as the original one, e.g., $X̄J→XJ$ as λ → 0 and for every J. However, this order may not be preserved, especially for degenerate excited states. Because the value of λ is small, we can reorganize the spectrum by matching the character of the excitation vectors. For example, suppose that the LUMO is twofold degenerate (or nearly degenerate), there is a LUMO-A and a LUMO-B, and that the first root is a HOMO → LUMO-A transition and the second is HOMO → LUMO-B. The perturbation slightly breaks this degeneracy and flips the order of the roots so the first root after solving the second problem is dominated by a HOMO → LUMO-B transition, while the second root is mainly HOMO → LUMO-A. We revert the order of these roots to ensure that the reorganized roots and their vectors are very close to the original ones (only differing by a small perturbation in the excitation amplitudes). In general, the sorting algorithm to reorganize the roots is based on the condition: Two pairs of eigenvectors $(X̄K,ȲK)$ and (X^J, Y^J) correspond to one another if $|(X̄K)TXJ−(ȲK)TYJ|>1−Cr$⁠, where C_r > 0 is a number close to zero. The best value of C_r is the minimum for which all the roots can be reorganized using the mentioned criterion.

Let us consider the root I and define the operator $Tˆ1I=∑iaσ(XaiσI+YaiσI)aˆaσ†aˆiσ$⁠. It can be shown that the Slater determinant describing the initial state of the perturbed, auxiliary system of electrons satisfies

here Ψ^I(λ) is the determinant composed of occupied, perturbed orbitals, {ζ_0,iσ(λ)}. The determinant |Φ〉 is the antisymmetrized product of ground-state orbitals, {φ_iσ}. The many-electron wave function |Ψ^I〉 is thus a coupled-clusters expansion in terms of singles only.

Let us expand linearly the exponential operator, $exp[λTˆ1I]≈I+λTˆ1I$ (this form is exact if the excitations take place from the same single occupied orbital). It is plausible to assume that the ground-state determinant is perturbed by a determinant associated to an excited state, where the intensity of this perturbation is determined by λ. Hence, an approximation to the excited-state wavefunction is

As discussed by Tavernelli, Curchod, and Rothlisberger,⁸ this type of auxiliary many-body wavefunction is useful to combine LR TDDFT calculations with the surface hopping algorithm, and it also provides the correct oscillator strengths for excitations from the ground state.

The auxiliary wavefunctions can be used to estimate the strength of the J←I transition,

This quantity can be calculated using the function,

In the Tamm-Dancoff approximation, we have that

where the right hand side is $〈ΦJ|μˆl|ΦI〉$⁠. Therefore,

The right hand side of the above equation is an approximation to $2/3|〈ΦJ|μˆ|ΦI〉|2$ if the Y vectors are included. The missing terms are products of Y vectors. We neglect these terms and use the above equation to estimate the auxiliary oscillator strengths.

In this section we apply the second LR method to the excited-state absorption of a diplatinum complex, pyrene, and azobenzene. This includes an application to a pump/probe experiment (for the diplatinum complex) exploiting the power of XANES.^9–12 The calculations illustrated here are non-relativistic and zero-order, i.e., the solvent effects are neglected. Regarding relativity, deep core excitations should be treated using non-local Breit theory or quantum electrodynamics; the orbital-level shifts often obtained from a scalar relativistic approximation such as ZORA are unsatisfactory.^13,14

We implemented the present second LR method in the NWChem suite.¹⁵ The procedure we apply in this section consists in running a ground-state geometry optimization and a LR calculation using the reference ground state. Then the target is selected and the perturbed orbitals, {ζ_0,pσ}, are generated to solve the second eigenvalue problem (Eq. (25)). From these calculations one extracts the eigenvectors to estimate the excited-state absorption profile, using Eq. (33). For these calculations we used λ = 1/2.

Let us consider the X-ray absorption of an optically excited dinuclear platinum(II) complex,¹⁶ [Pt(ppy)(μ-^tBu₂pz)]₂ (ppy = 2-phenylpyridine, ^tBu₂pz = 3, 5-di-tert-butylpyrazolate). In the X-ray experiment the molecules are first subject to a transient visible absorption by a laser pulse at 527 nm, inducing the S₁←S₀ transition, dominated by a transition from the HOMO to the LUMO. The valence excited state is then probed by an X-ray pulse, with duration around 100 ps, right after the valence excitation. The HOMO in the S₀ state is a dσ^∗ anti-bonding orbital between the platinum atoms. The HOMO → LUMO transition results in a transfer of charge from the platinum atoms to the ligands. To determine the evolution of the excited molecules in the S₁ state, a different study was carried out (measuring optical transient absorption anisotropy) to investigate the first few picoseconds after the valence excitation.¹⁷ In the resulting proposed mechanism, the molecule is excited to the S₁ state, then undergoes internal conversion to a T_1b state (this process takes about 145 fs), which decays to a T_1a configuration in about 2.4 ps. This T_1a has a lifetime around 700 ns and its structure was investigated by XANES.¹⁶ 

The structural changes of the diplatinum complex caused by the laser perturbation are determined using the X-ray probe pulse, which promotes a 2p electron to the vacated dσ^∗ level of the T_1a state. Due to the depletion caused by the laser, the effective Pt–Pt bond order increases, resulting in a contraction of the Pt–Pt distance, by about 0.2 Å, in the T_1a state.¹⁶ The first sharp peak around 11 565 eV in the XANES difference spectrum (laser on minus laser off) corresponds to the orbital promotion and was not covered by the scattering theory employed in Ref. 16. Describing the filling of the T_1a dσ^∗ hole can be achieved using unrestricted ground-state and TDDFT calculations. In addition, here we show that if the excited S₁ state is probed by X-ray absorption, then the orbital promotion should be observed as well. The S₁ state corresponds to the “LUMO” occupied by two electrons, and the “HOMO” level is empty.

For the calculations we use the LDA0 XC functional,¹⁸ a combination of 25% HF exchange, 75% local exchange, and 100% LDA correlation. We chose this functional because of its low computational cost and improved description of the HOMO-LUMO gap. The basis set used, DZP, offers a good balance of accuracy and speed. The target excitation is that of the orbital promotion. Other nearby transitions are negligible compared to that of refilling the vacated dσ^∗ level. Thus the DZP set is convenient in this case.

In Fig. 1(b) we show the comparison between the calculated and experimental difference spectra. The latter is calculated as the difference between the XANES spectrum taken at 100 ps after the excitation and the one taken from the ground state (S₀). On the other hand, we obtained two theoretical difference spectra, calculated by comparing the absorption profiles of the target excited (T_1a or S₁) and ground states. We optimized the triplet T_1a geometry before running the unrestricted TDDFT calculations. The second LR calculations with respect to the S₁ state were performed using the ground-state geometry.

The second LR calculations capture the first large oscillation in the spectrum, in agreement with the unrestricted LR TDDFT calculations and with the chemical argument made in Ref. 16 (that the first peak is due to the filling of the vacated dσ∗ orbital). The calculations indicate that the molecule absorbs energy around 2 eV from the laser. Then the system, when probed, will absorb for the 2p → dσ^∗ transition. The excitation energy of this transition is slightly less than that associated to the transition from the 2p level to the LUMO of the ground state. This last transition has a similar oscillator strength to that of the 2p → dσ^∗ transition. The overall effect results in the observed oscillation in the difference spectrum. The remaining features, additional oscillations, of the experimental difference spectrum (including the second positive peak) are attributed to the local structure of the triplet state.¹⁶ The optimized singlet and triplet geometries differ by a small extent. This makes the agreement between the second and unrestricted LR calculations quite evident. It must be remarked, however, that the relaxation of the T_1a and S₁ geometries should induce changes in the absorption profiles of these excited states.

The next application is the transient absorption spectrum of pyrene.¹⁹ In Ref. 19, pyrene (in n-octane) is excited at 290 nm, promoting a transition to the second excited singlet, and first probed after 60 ps of the excitation; and in Ref. 20, pyrene (in 1,1,2,2-tetrachloroethane) is excited at 266 nm and then probed after 0.7 ps. The lifetime of the pyrene excited singlet is 960 ns. For this molecule we apply the Coulomb-attenuated method to the LDA0 functional, i.e., we add a non-local long-range electron-electron interaction to the GKS Hamiltonian.²¹ This method improves the estimation of excitation frequencies.²² The resultant functional, CAM-LDA0,²³ is slightly less expensive and performs well with respect to the CAM-B3LYP functional.¹⁸ We optimized the geometry by using the LDA0 XC functional and the basis set 6-31G*. The excitation calculations for the second excited singlet were performed at the cc-pVDZ level. Fig. 2(b) shows the experimental transient absorption spectra of pyrene and our calculations. The theoretical spectrum from the second LR calculations shows a defined peak at 480 nm, and the absorption around 370 nm is reproduced by the calculation. In addition, the spectrum shown in yellow, Fig. 2(b), was calculated using the auxiliary wavefunctions and Eq. (39). The main peak of this spectrum, around 400 nm, is further scaled for comparison. This peak is in fact about fifty times larger than the peak around 375 nm from the second LR calculations, while the peak around 480 nm is missing. The transition intensities in both cases, however, do not agree completely with experimental data. Inclusion of excited-state dynamics, non-adiabatic effects, and coupling with the bath could improve the results.

For the third example we apply the same methodology shown above for azobenzene. Fig. 3(b) shows the comparison between the experimental data²⁴ and second LR calculations. Azobenzene, in hexane, is pumped into the second excited singlet and then probed at 60 fs after the laser pump pulse. The measurements show a wide absorption band around 465 nm, while the calculated one is red-shifted by 35 nm. Additionally, the experimental band peaked at 400 nm differs from the computed one by 50 nm. In contrast, the spectrum calculated using $fJIaux$ (solid yellow line) shows a wide horizontal gap between the left and right bands, which is reduced after calculating f_JI. The relative heights between the two bands are improved by the second LR calculations.

We developed a theoretical method to study sequential double excitations. The approach consists in generating, from a LR-TDDFT calculation, occupied and virtual orbitals that are perturbed by a target excited state, and then using these orbitals to rerun the LR-TDDFT calculations and calculate the corresponding oscillator strengths and transition moments. The computational cost of the second LR calculations is the same of regular LR-TDDFT computations and lower than those of methods solving the full quadratic response problem. Despite the lengthy, formal justification, our method is simple and can be extended to other wave-function methods, such as coupled clusters, and include condensed phase dynamics.

We acknowledge support for this work from the Ultrafast Initiative of the U.S. Department of Energy (DOE), Office of Science, Office of Basic Energy Sciences, through Argonne National Laboratory under Contract No. DE-AC02-06CH11357. For vis-UV development we acknowledge the DOE, under Grant No. DOE DE-FG02-10ER16153. M.A.M. is grateful to Carlos Borca, Nicholas Jackson, and Dr. Inkoo Kim and Dr. Amrit Poudel, for productive discussions.

Up to first order, the TD wavefunction is

where $|Ψ(0)(t)〉=exp(−iHˆ0t)|Ψ(0)〉$ and

Recall that the initial state is |Ψ(0)〉 = |Ψ₀〉 + λ|Ψ_I〉. Let us write $Ψ(1)(t)=∑JcJ(t)ΨJ$⁠. The lth component of the averaged dipole is (up to first order)

Apply 〈Ψ_J|(id/dt) on both sides of Eq. (A2) and Fourier transform. This leads to