Interpretation of Coupled-Cluster Many-Electron Dynamics in Terms of Stationary States

We demonstrate theoretically and numerically that laser-driven many-electron dynamics, as described by bivariational time-dependent coupled-cluster (CC) theory, may be analyzed in terms of stationary-state populations. Projectors heuristically defined from linear response theory and equation-of-motion CC theory are proposed for the calculation of stationary-state populations during interaction with laser pulses or other external forces, and conservation laws of the populations are discussed. Numerical tests of the proposed projectors, involving both linear and nonlinear optical processes for He and Be atoms and for LiH, CH+, and LiF molecules show that the laser-driven evolution of the stationary-state populations at the coupled-cluster singles-and-doubles (CCSD) level is very close to that obtained by full configuration interaction (FCI) theory, provided that all stationary states actively participating in the dynamics are sufficiently well approximated. When double-excited states are important for the dynamics, the quality of the CCSD results deteriorates. Observing that populations computed from the linear response projector may show spurious small-amplitude, high-frequency oscillations, the equation-of-motion projector emerges as the most promising approach to stationary-state populations.


INTRODUCTION
Providing unique time-resolved insights into electronic quantum dynamics and with the exciting prospect of detailed manipulation and control of chemical reactions, 1 increasing experimental and theoretical research efforts have been directed toward attosecond science in the past couple of decadessee, for example, ref 2 for a recent perspective. While the initial step usually involves ionization induced by extremeultraviolet or near-infrared laser pulses, Hassan et al. 3 have demonstrated that optical attosecond pulses may be used to observe and control the dynamics of bound electrons with little or no ionization probability. Whether ionization plays a role or not, the rapid development of experimental methodology creates a strong demand for explicitly time-dependent quantum chemical methods that can accurately simulate the ultrafast many-electron dynamics driven by ultrashort laser pulses.
While real-time time-dependent density functional theory 4−6 is a highly attractive option from the viewpoint of computational efficiency, it suffers from a number of deficiencies caused largely by the reliance on the adiabatic approximation in most practical applications. 7 Improved accuracy can be achieved with wave function-based methods at the expense of increased computational costs. 7 In a finite basis, the exact solution to the time-dependent electronic Schrodinger equation is the full configuration interaction (FCI) wave function whose computational complexity, unfortunately, increases exponentially with the number of electrons. We are thus forced to introduce approximations. The perhaps most widely used time-dependent wave function approximation for simulating many-electron dynamics is multiconfigurational time-dependent Hartree− Fock (MCTDHF) theory 8−11 and the related time-dependent complete active space self-consistent field and restricted active space methods. 11−14 Restricting the participating Slater determinants to those that can be generated from a fixed number of electrons and a carefully selected active space of (time-dependent) spin orbitals, these methods still have the FCI wave function at the heart, eventually facing the exponential scaling wall as the number of active electrons and orbitals is increased.
Coupled-cluster (CC) theory offers a gentler polynomialscaling hierarchy of approximations that converge to the FCI wave function. Besides the differences in computational complexity, the two methods differ in the sense that MCTDHF captures static (strong) correlation, whereas single-reference CC theory aims at dynamical correlation effects. Yielding energies, structures, and properties with excellent accuracy for both ground and excited states of weakly correlated systems, it has become one of the most trusted methods of molecular quantum chemistry. 15 Recent years have witnessed increasing interest in time-dependent CC (TDCC) theory 16−20 for numerical simulations of many-body quantum dynamics in nuclear 21 and atomic and molecular 22−38 systems. In addition, TDCC theory has played a key role in a recent work on finitetemperature CC theory for molecular 39,40 and extended 41 systems. While the papers by Christiansen and co-workers 35,36 are concerned with vibrational CC theory and those of Pigg et al. 21 with nucleon dynamics, the remaining papers are focused on the dynamics of atomic and molecular electrons exposed to electromagnetic fields such as ultrashort laser pulses.
In many cases, the main goal is to compute absorption (or emission) spectra [24][25][26][27][28]37,38 by Fourier transformation of the induced dipole moment. This requires the calculation of the induced dipole moment for extended periods of time after the perturbing field or laser pulse has been turned off. While decisive for the features observed in the final spectrum, the dynamics during the interaction with the laser pulse is rarely analyzed in detail. Processes that occur during the pulse, such as high harmonic generation and ionization, are studied using TDCC theory in refs. 31,34 . Since energy is the physical quantity associated with time translations, textbook analyses of such interactions are naturally performed in terms of the population of the energy eigenstatesthe stationary statesof the fieldfree particle system, see, for example, ref 42. However, manybody theories such as TDFCI, MCTDHF, and TDCC theories do not express the wave function as a superposition of stationary states, making the analysis difficult to perform in simulations. Moreover, when approximations are introduced (truncation of the many-body expansion), the stationary states are hard to define precisely for nonlinear parameterizations such as MCTDHF and TDCC theories. The problem is particularly pronounced for approximate methods where the orbitals are time-dependent such that a different subspace of the full configuration space is spanned in each time step of a simulation, leading to energies and eigenvectors of the Hamilton matrix that vary depending on the laser pulse applied to the system. 43 This implies, for example, that identification of stationary-state energies by Fourier transformation of the postpulse autocorrelation function leads to pulse-dependent results. Still, several reports of population transfer during interaction with laser pulses have been reported recently 44−46 using MCTDHF theory.
The natural approach would be to define the stationary states from the zero-field Hamiltonian and zero-field wave function using, for example, linear response theory 47 or orthogonality-constrained imaginary time propagation. 48 The latter approach was investigated recently within the framework of MCTDHF theory by Loẗstedt et al., 49 who found that the stationary-state populations oscillate even after the pulse is turned off unless a sufficiently large number of active orbitals is included in the wave function expansion. In this work, we use both CC linear response (CCLR) theory 19,50,51 and equationof-motion CC (EOMCC) theory 52−55 to propose projectors whose expectation values yield stationary-state populations. Test simulations are presented with different laser pulses, and the TDCC results are compared with the exact (TDFCI) results.
The paper is organized as follows. In Section 2, we briefly outline the exact quantum dynamics on the basis of energy eigenstates and use analogies to propose projectors whose expectation values can be interpreted as stationary-state populations within TDCC theory. Technical details of the numerical simulations are given in Section 3, and numerical results are presented and discussed in Section 4 for atoms and diatomic molecules in few-cycle laser pulses, including chirped pulses. Concluding remarks are given in Section 5.

THEORY
2.1. Recapitulation of Exact Quantum Dynamics. Laser-driven quantum dynamics of a particle system is usually interpreted in terms of stationary states |n⟩ defined as solutions of the time-independent Schrodinger equation where H 0 is the time-independent Hamiltonian of the particle system and E n is the energy of the stationary state |n⟩. The stationary states evolve in time according to and are assumed to form a complete orthonormal set such that P n n n t n t P ( ) ( ) , 1 where 1 is the identity operator. Note that the continuum is formally included in the summation over states. The time evolution of the particle system is determined by the time-dependent Schrodinger equation [using atomic units (a.u.) throughout], where |Ψ(t)⟩ is the normalized state of the system with a known initial value |Ψ 0 ⟩ and the dot denotes the time derivative. Within semiclassical radiation theory, the time dependence of the Hamiltonian, stems from the operator V(t) describing the interaction between the particle system and external electromagnetic fields.
Expressing the time-dependent state as a superposition of stationary states, where C n (t) = ⟨n(t)|Ψ(t)⟩, the time-dependent Scrodinger equation may be recast as an ordinary differential equation with the initial conditions C n (0) = ⟨n|Ψ 0 ⟩. The population of stationary state |n⟩ at any time t ≥ 0 may be determined as the expectation value of the projection operator P n Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article 8) which is real and non-negative. Since the state is assumed normalized, ⟨Ψ(t)|Ψ(t)⟩ = 1, the stationary-state populations sum up to one, n ∑ p n (t) = 1, and, hence, are bound from above as well as from below: 0 ≤ p n (t) ≤ 1.
It follows from the time-dependent Schrodinger equation that the expectation value of some operator, say B, evolves according to the Ehrenfest theorem, where the last term vanishes when B is a time-independent operator such as a stationary-state projector. Hence, stationarystate populations are conserved in the absence of external forces as the projectors commute with the time-independent Hamiltonian, [P n ,H 0 ] = 0. Stationary-state populations are required to identify transient phenomena such as Rabi oscillations in simulations and may be used to determine the composition of the quantum state resulting from the application of a short pump laser, facilitating interpretation of spectra recorded by means of a subsequent probe laser. The stationary-state populations can be controlled by varying the pump laser parameters such as peak intensity, shape, and duration. Predicting final populations with varying laser parameters thus becomes a central computational goal.
Unfortunately, computing all stationary states of a particle system followed by integration of the time-dependent Schrodinger equation presents an insurmountable challenge, even if the number of stationary states is kept finite through the use of a finite set of basis functions. In practice, the quantum state is parameterized in a finite-dimensional Hilbert space spanned by well-defined basis vectors (Slater determinants, in the case of electronic systems) rather than stationary states. Populations of a few selected (low-lying) stationary states can then be computed as a function of time with varying pump laser parameters. This procedure is easily implemented for any electronic-structure method with explicit parameterization of orthogonal ground-and excited-state wave functions and has been used recently by, for example, Peng et al. 56 within time-dependent configuration-interaction theory to predict populations of stationary electronic states of the rigid decacene molecule with varying laser parameters. In the following, we will present an approach to the calculation of stationary-state populations within the framework of TDCC theory.
2.2. TDCC State Vector. Our starting point is Arponen's time-dependent bivariational formulation of CC theory 20 within the clamped-nucleus Born−Oppenheimer approximation. This allows us to parameterize the CC ket and bra wave functions as independent approximations to the FCI wave function and its Hermitian conjugate. The quantum state of an atomic or molecular many-electron system at time t is then represented by the TDCC state vector 29 i k j j j j j j y where the component functions are defined by Although the reference determinant |Φ 0 ⟩ should be constructed from time-dependent orthonormal 16,31,57 or biorthonormal 23,30,58 orbitals to capture the main effects of interactions between the electrons and external fields, we shall in this work use the static Hartree−Fock (HF) ground-state determinant for simplicity. As long as the external field does not lead to nearly complete depletion of the ground state, we found in ref 30 that the results obtained with static and dynamic orbitals are virtually identical. In addition, using static HF reference orbitals allows us to exploit well-known CC theories for excited states, as discussed in more detail below, and we avoid the complexity of computing overlaps between determinants in different nonorthogonal orbital bases.
The cluster operators are defined as where μ > 0 labels excitations out of the reference determinant, i.e., If the cluster operators include all excited determinants, the CC state becomes equivalent to the exact wave function, the FCI wave function. Approximations are obtained by truncating the cluster operators after singles to give the CCS method, after singles and doubles to give the CC singles-and-doubles (CCSD) model, and so on. Since the HF reference determinant is static, the time dependence of the cluster operators is carried by the amplitudes τ μ (t) and λ μ (t) only. The amplitude τ 0 (t) is a phase parameter related to the so-called quasi-energy, 59,60 and λ 0 (t) determines the normalization of the state, as discussed in more detail below.
In the closed-shell spin-restricted CCSD model, the singlesand-doubles excitation and de-excitation operators are defined as 61 where i, j and a, b refer to occupied and virtual spatial HF orbitals, respectively, and is a unitary group generator expressed in terms of the elementary second-quantization spin−orbital (α and β here refer to the spin-up and spin-down states) creation and annihilation operators. The equations of motion for the amplitudes are derived from the time-dependent bivariational principle and are given by 20,29 i Ht  (19) where the dot denotes the time derivative and where While H 0 is the time-independent molecular electronic Hamiltonian in the clamped-nucleus Born−Oppenheimer approximation, V(t) describes the interaction of the electrons with explicitly time-dependent external fields in the semiclassical approximation. Note that the normalization amplitude λ 0 is constant. The equations of motion ( 18 ) and ( 19 ) must be integrated with suitable initial conditions. In this work, we use the CC ground-state i k j j j j j j j y where T X e , The ground-state amplitudes satisfy the stationary CC equations (25) and τ 0 (t = 0) = 0 such that, in the absence of external perturbations, the time-dependence of the TDCC state vector correctly becomes |S(t)⟩⟩ = |S 0 ⟩⟩ exp(−iE 0 t), where E 0 is the CC ground-state energy with respect to which the TDCC state vector is normalized, i.e., provided we choose Re(λ 0 ) = 1. In practice, we choose λ 0 = 1. The indefinite inner product induces the expectation value expression 29 where the two-component form of the quantum mechanical operator C reads While the expectation value of an anti-Hermitian operator C † = −C is imaginary, the expectation value of a Hermitian operator C † = C is real. This symmetrized form of the CC expectation value was first introduced by Pedersen and Koch 62 in order to ensure correct symmetries, including time-reversal symmetry, of CC response functions. Using the expectation value expression to compute the electric dipole moment induced by an external laser field, absorption spectra can be obtained by Fourier transformation.
Given a set of orthonormal excited-state vectors |E n ⟩⟩, which are orthogonal to the ground state with respect to the indefinite inner product, we may define the projection operator and compute the population p n (t) of excited state n at time t as the expectation value This would provide a time-resolved picture of the populations of excited states within TDCC theory. Unfortunately, a fully consistent set of CC excited-state vectors is not known.
There are two distinct approaches to excited states in common use within CC theory today. 63 One is the EOMCC 52−55 approach where the excited states are parameterized explicitly in terms of linear excitation and deexcitation operators, which generate the excited-state vectors from the ground state. While making it straightforward to express the projection operator in eq 31, the linear Ansatz of EOMCC theory leads to size-intensivity issues in transition moments and response properties such as polarizabilities. 64−66 The other approach is CC response theory 19,50,51 where the time-dependent Schrodinger equation is solved in the frequency domain using a combination of Fourier transformation and adiabatic perturbation theory. The amplitudes at each order are then used to express linear, quadratic, and higher-order response functions. This leads to the identification of excitation energies and one-and multiphoton transition moments by analogy with response theory for exact stationary states. 47 It does not, however, lead to explicit expressions for the excited-state vectors, making it impossible to construct projection operators of the form ( 31 ). While the excitation energies from EOMCC theory and CC response theory are identical and properly size-intensive, the transition moments differ and no size-extensivity issues are present in CC response theory. 64,65 We note in passing that Hansen et al. 36 offer a detailed discussion of size-extensivity (and, by extension, sizeintensivity) and size-consistency issues within CC theory using the more concise and accessible concepts of additive and multiplicative separability.
2.4. Projection Operators from EOMCC Theory. Exploiting that EOMCC theory provides an explicit parameterization of "left" (bra) and "right" (ket) excited states, we investigate the EOMCC projector defined as  (34) and the linear excitation and de-excitation operators are given by The EOMCC excitation and de-excitation operators are truncated at the same level as the underlying CC ground-state cluster operators T 0 and Λ 0 . The ground-state component reads The EOMCC amplitudes are determined from the non-Hermitian eigenvalue problem where 1 is the unit matrix and ΔE is a diagonal matrix with the excitation energies ΔE n = E n − E 0 as elements. The elements of the non-Hermitian Jacobian matrix are defined by Note that the components of the EOMCC excited-state vector are biorthonormal n m nm δ ⟨Ψ ∼ |Ψ ⟩ = (39) and orthogonal to the CC ground state in the sense of This implies that the EOMCC projectors in eq 33 are Hermitian with respect to the indefinite inner product, annihilate the ground state, Pn|S 0 ⟩⟩ = 0, and are idempotent and orthogonal, P P P n m nm n δ̂=̂ (41) In the limit of untruncated cluster operators, it is readily verifiedusing the orthonormality of the left and right eigenvectors in eq 37that the EOMCC projectors satisfy the completeness relation If the cluster operators are truncated, the right-hand side must be corrected for excited determinants beyond the CC truncation level (e.g., triples, quadruples, etc. for EOMCCSD).
The one-photon transition strength obtained from the EOMCC projector is given by where B and C are Hermitian operators representing electric or magnetic multipole moments and B̂and Ĉare their twocomponent forms defined in eq 30. The transition strength is properly symmetrized with respect to simultaneous permutation of the multipole operators and complex conjugation, ⟨⟨S 0 | ĈPnB̂|S 0 ⟩⟩* = ⟨⟨S 0 |B̂PnĈ|S 0 ⟩⟩ and agrees with the commonly used expression in EOMCC theory, which is based on a configuration-interaction-like interpretation of the bra and ket states. This expression yields the correct FCI limit but, as mentioned above, the transition strength is not properly sizeintensive when the cluster operators are truncated. 64−66 We may now use eq 33 to extract excited-state populations from the TDCC state vector according to eq 32, While the EOMCC excited-state populations are manifestly real, they are neither bounded above by 1 nor below by 0. Lack of proper bounds is common in CC theory, but problems are rarely experienced in practical calculations as long as the CC state vector is a sufficiently good approximation to the FCI wave function. This, in turn, is related to the distance in Hilbert space between the reference determinant and the FCI wave function, as discussed in more detail by Kristiansen et al. 30 for TDCC theory.
Consistent with this analysis, we compute the ground-state population via the ground-state projector i k j j j j j j j j y 2.5. Projection Operators from CC Linear Response Theory. In lieu of explicitly defined excited states in CC response theory, we investigate the CCLR projector i k j j j j j j j j y { z z z z z z z z P where we have introduced the notation (f and g are arbitrary functions) The functions |Ψ n ⟩ and ⟨Ψ̃n| are defined in eq 34, and n n n 0 0 The amplitudes n ̅ μ are determined by the linear equations 51 where the superscript T denotes matrix transposition, 1 is the unit matrix, A is the CC Jacobian matrix defined in eq 38, ΔE n is the eigenvalue (excitation energy) corresponding to the right eigenvector n (cf. eq 37), and the symmetric matrix F is defined by The main arguments in favor of the CCLR projector ( 47 ) are as follows: Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article 1. it becomes identical to the EOMCC projector in the FCI limit, and 2. it yields the correct size-intensive CCLR ground-toexcited state transition strength. We demonstrate in the Appendix that ⟨Ψ̑n| = 0 in the FCI limit, implying that the CCLR projector ( 47 ) becomes identical to the EOMCC projector ( 33 ) in this limit. The CCLR projector ( 47 ) gives the correct transition strength, i.e., To make the equivalence evident, we note that in the notation used by Christiansen et al. 51 for "right" and "left" This allows us to recast the transition strength as which is identical to the expression obtained as a residue of the CCLR function in eq 5.
The correction term only vanishes in the FCI limit. They do, however, project onto the orthogonal complement of the CC ground state as i k j j j j j j j j y Excited-state populations extracted by the CCLR projectors according to eq 32 become While real, the CCLR population is neither bound below by 0 nor bound above by 1. The ground-state population is given by eq 46.
2.6. Conservation Laws. As noted in Section 2.1, exact stationary-state populations are conserved in time intervals where no external forces act on the particle system. Conservation laws in the framework of TDCC theory have been discussed at various levels of detail previously. 21,29,35,38,67,68 To this end, we recast the TDCC equations of motion in the Hamiltonian form 29 where the Hamiltonian function t ( , , ) τ λ = is defined by Introducing the Poisson-like bracket for any analytic function f = f(τ,λ,t) and g = g(τ,λ,t), the relation is readily obtained from the Hamiltonian eq 61.
which shows that energy is conserved whenever the Hamiltonian operator is constant in time, including before and after the application of external forces such as laser pulses. Note that this is true regardless of the truncation of the cluster operators and regardless of the initial conditions of the amplitudes. This observation was used in ref 29 to propose symplectic numerical integration as a stable method for solving the TDCC equations of motion. The time evolution of the TDCC expectation value in eq 29 is obtained by choosing f = ⟨Ψ̃|C|Ψ⟩. Using a derivation analogous to that of Skeidsvoll et al., 38 we find where we have introduced the projected commutator with Here, the summation over excited determinants is truncated at the same level as the cluster operators. The time evolution of the TDCC expectation value thus is Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article This is not quite the form of a generalized Ehrenfest theorem since, in general, Consequently, constants of motion in exact quantum dynamics are not necessarily conserved in truncated TDCC theory. In the FCI limit, however, P ∥ = 1 and the Ehrenfest theorem is recovered from eq 69.
The EOMCC stationary-state population, eq 44, is of the form ( 29 ) with C = |Ψ n ⟩⟨Ψ̃n|, and the time evolution, therefore, is given by eq 69. The proposed EOMCC projector thus breaks the conservation law of stationary-state populations in truncated TDCC simulations, although we note that it is properly restored in the FCI limit since C = |Ψ n ⟩⟨Ψ̃n| commutes with H 0 according to eqs 81 and 82 in the Appendix. It seems reasonable to expect that the conservation law is approximately fulfilled whenever the many-electron dynamics predominantly involves stationary states that are well approximated within (truncated) EOMCC theory.
The CCLR stationary-state population, eq 60, is not of the form ( 29 ). Instead, the time evolution is given by eq 64 with which only depends on time through the amplitudes (∂f/∂t = 0). Hence, the CCLR projector also breaks the conservation law of stationary-state populations in truncated TDCC simulations. It is restored in the FCI limit where the CCLR and EOMCC projectors are identical. Again, it seems reasonable to expect that the conservation law is approximately fulfilled whenever CCLR theory provides a sufficiently good approximation to the FCI states.

COMPUTATIONAL DETAILS
Explicitly time-dependent simulations are performed with a closed-shell spin-restricted TDCCSD Python code generated using a locally modified version of the Drudge/Gristmill suite for symbolic tensor algebra developed by Zhao and Scuseria. 69 The static HF reference orbitals and Hamiltonian integrals are computed by the Dalton quantum chemistry program 70,71 along with the response vectors required for the EOMCC and CCLR projectors using the implementations described in refs. 72−75 Tight convergence criteria are employed: 10 −10 a.u. for the HF orbital gradient norm (implying machine precision for the HF ground-state energy) and 10 −8 a.u. for the CCSD residual norms (both ground state and response equations). Dunning's correlation-consistent basis sets, 76−78 downloaded from the Basis Set Exchange, 79 are used throughout. We also perform TDFCI simulations using the contraction routines implemented in the PySCF package. 80 The TDCCSD and TDFCI equations of motion are integrated using the symplectic Gauss−Legendre integrator 81 as described in ref 29.
To test the proposed EOMCC and CCLR projectors, TDCCSD and TDFCI simulations are carried out for the He and Be atoms placed at the coordinate origin. Further tests are performed for the LiH molecule placed on the z-axis with the Li atom at the origin and the H atom at z = 3.08 a.u. and for the CH + ion placed on the z-axis with the C atom at the origin and the H atom at z = 2.13713 a.u. Finally, we study the time evolution of stationary-state populations during the optical pump pulse applied by Skeidsvoll et al. 38 to investigate transient X-ray spectroscopy of LiF. As in ref 38, we place the LiF molecule on the z-axis with the F atom at the origin and the Li atom at z = −2.9552749018 a.u. All electrons are correlated, and point group symmetry is not exploited in these simulations.
We assume the systems are initially in the ground state and expose them to a laser pulse described by the semiclassical interaction operator in the electric dipole approximation where d is the electric dipole operator of the electrons, u is the real unit polarization vector of the electric field, and ε(t) is the time-dependent electric-field amplitude.
Two forms of the electric-field amplitude are used in this work. One is the sinusoidal pulse where ε 0 is the field strength, ω 0 is the carrier frequency of the pulse, and t 0 is the time at which the pulse is turned on. The time-dependent phase ϕ(t) may effectively alter the instantaneous carrier frequency, creating a chirped laser pulse, if it depends on time at least quadratically. 82 In this work, we use the quadratic form which, for b ≠ 0, creates a linearly chirped laser pulse with instantaneous frequency ω(t) = ω 0 + a + b(t − t 0 ). The envelope G(t) controls the shape and duration of the pulse and is defined by where t d is the duration of the pulse and Θ(t) is the Heaviside step function. The second form of the electric-field amplitude used in this work is the Gaussian pulse with the envelope where t 0 is the central time of the pulse and σ is the Gaussian root-mean-square (rms) width and N defines the start time (t 0 − Nσ) and end time (t 0 + Nσ) of the pulse through the Heaviside step functions. Note that N thus introduces discontinuities of the Gaussian pulse at each end. Unless a sudden disturbance of the system is intended, one must choose N large enough that the discontinuities are negligible.  29,30 Dynamic orbitals, such as those of orbital-adaptive TDCC theory, 23 are required for a numerically stable integration of the equations of motion. 30 Consequently, in this work, we will focus on Rabi oscillations between excited states. Rabi oscillations between two excited states can be achieved by the application of two consecutive laser pulses, the first of which is resonant with a dipole-allowed transition from the ground state, while the second is resonant with a dipoleallowed transition between the resulting excited state and another one. The intensity and duration of the first pulse must be such that the ground-state population is significantly reduced but not entirely depleted. Nonlinear optical processes are thus involved, making it an ideal test case for the CCLR and EOMCC projectors, which are constructed on the basis of first-order perturbative arguments (first-order perturbation theory in the case of CCLR and linearization of the cluster exponential in the case of EOMCC), which cannot necessarily be expected to correctly capture higher-order optical processes. In particular, transition moments between excited states are quadratic response properties, which cannot be expressed solely in terms of linear response parameters. 51 It is, therefore, important to test if the proposed projectors correctly capture the effects of nonlinear optical processes.

RESULTS
Results for the He atom with the aug-cc-pVTZ basis set are presented in Figure 1. The integration of the TDCCSD equations of motion was performed with time step Δt = 0.1 a.u. = 2.42 as using the eighth-order (s = 4) Gauss−Legendre integrator and a convergence threshold of 10 −10 (residual norm) for the fixed-point iterations. The ground-and excitedstate energy levels shown in the left panel of Figure 1 were computed using CCSD linear response theory. In total, 14 excited states were computed, several of which lie above the ionization energy, which is estimated to be 0.902 a.u. using the total ground-state energy difference between the neutral and ionized atoms at the (spin unrestricted) CCSD/aug-cc-pVTZ level of theory. 83 Although the states above the ionization energy are unphysical, we keep them for the purpose of comparing with regular TDFCI simulations with the same basis set.
The right panel of Figure 1 shows the total energy-level populations computed using the CCLR projector as a function of time during the application of two consecutive laser pulses. Energy levels that are never populated above 0.01 are excluded. The level populations are computed by summing up the populations of all states belonging to each energy level, thus avoiding ambiguities arising from the arbitrariness of the basis of a degenerate subspace. As expected, the first laser pulse causes significant population of the 2 1 P energy level. The highlying 4 1 D level, which is dominated by the 1s 1 3d 1 electron configuration and located 0.972 a.u. above the 2 1 P level, also becomes populated toward the end of the first pulse because of the dipole-allowed transition 2 1 P → 4 1 D. The length-gauge oscillator strength of this transition is f = 0.484 compared with f = 0.355 for the 0 1 S → 2 1 P transition. The population of the 4 1 D level is still modest, however, since the transition can only occur once the 2 1 P level is significantly populated.
The second pulse induces several cycles of Rabi oscillations between the 2 1 P and 1 1 S levels. The Rabi oscillations are slightly perturbed by weak transitions between the 2 1 P level and the 0 1 S ground state, as witnessed by the increasing oscillation of the ground-state population when the population of the 2 1 P level is close to its maximum value. We also observe an even weaker perturbation caused by the higher-lying 4 1 D level. The CCLR populations agree both with TDFCI populations and EOMCC populations: the rms deviation for the entire simulation is 10 −3 between the CCLR and TDFCI populations and 3 × 10 −7 between CCLR and EOMCC populations. We have previously demonstrated that discrepancies with TDFCI simulations for the He atom can be reduced by tightening the computational parameters such as convergence thresholds and, most importantly, by reducing the time step of the numerical integration. 30 We thus ascribe the small discrepancies between the CCLR and TDFCI populations to the rather coarse discretization (Δt = 0.1 a.u.) employed in the numerical integration scheme and conclude that the proposed CCLR and EOMCC projectors behave correctly in the FCI limit.
It is of interest to compare the TDCCSD simulation with a much simpler model based on an eigenstate expansion propagated according to eq 7. Letting |n⟩,n = 0,1,2,...,14 represent the stationary states computed with CCSD linear response theory, all we need to integrate eq 7 in the presence of external laser pulses of the form ( 72 ) is the dipole matrix in the energy eigenbasis. This model is essentially identical to that employed by Sonk et al., 84 except that we use CCSD linear and quadratic response theories rather than EOMCCSD theory to build the dipole matrix. (Note, however, that the CCSD response and EOMCCSD approaches yield identical results for He.) It is also similar to the EOMCCSD model employed by Luppi and Head-Gordon, 85 who propagated both bra and ket states.
The only obstacle is that the "right" and "left" transition moments from CC response theory, cf. eqs 54 and 55, are not related by complex conjugation, yielding a spurious non-Hermiticity of the dipole matrix. Sonk et al. 84 and Luppi and Figure 1. Energy-level populations computed with the aug-cc-pVTZ basis set through CCLR projectors plotted as functions of time for He exposed to two consecutive laser pulses, the first resonant with the 0 1 S → 2 1 P transition and the second resonant with the 2 1 P → 1 1 S transition, with peak intensity 87.7 TW/cm 2 . The black curves show the populations during the TDFCI simulation.
Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article Head-Gordon 85 circumvented this issue by only using the Hermitian part of the matrix. In the present case, due to symmetry, the CCSD electric-dipole transition moments are either zero or may be chosen parallel to one of the three Cartesian axes, making it easy to define the off-diagonal dipole matrix elements as the negative square root of the dipole transition strength (cf. eq 56 with B = C a Cartesian component of the position operator). With the dipole matrix thus constructed for the 15 states, we have integrated eq 7 with the initial condition C n (t = 0) = δ n0 and the same consecutive laser pulses as in the He/aug-cc-pVTZ TDCCSD simulations. The Gauss−Legendre integrator was used with the same parameters and time step. The rms population deviation between the model and the full TDCCSD simulations is about twice that between the TDFCI and TDCCSD simulations. The maximum absolute deviation is an order of magnitude greater (0.02 vs 0.002), however, indicating that a potentially large number of states, including states above the ionization energy, may be required for the simple model to agree quantitatively with full TDCCSD simulations in general.
Moving away from the FCI limit, we repeat the study of excited-state Rabi oscillations for the Be atom with the aug-cc-pVDZ basis set. The integration parameters were the same as those used for He above. The ground-and excited-state energy levels shown in the left panel of Figure 2 were computed using CCSD linear response theory. In total, 21 excited states were computed, and the CCSD excitation energies agree with those of FCI theory to within 3 × 10 −4 a.u. Several of the computed excited states lie above the ionization energy, which is estimated to be 0.341 a.u. using the total ground-state energy difference between the neutral and ionized atoms at the (spin unrestricted) CCSD/aug-cc-pVDZ level of theory. 83 The highlying excited states are retained in order to compare with regular TDFCI simulations with the same basis set.
The right panel of Figure 2 shows the total energy-level populations computed using the CCLR projector as a function of time during the application of the two consecutive laser pulses. Energy levels that are never populated above 0.01 are excluded. The first laser pulse causes significant population of the 1 1 P energy level by excitation from the ground state (oscillator strength f = 0.478), although the transition is quenched by further excitation to the high-lying 7 1 D level from the 1 1 P level (f = 0.685). The 7 1 D level, which is dominated by the 1s 2 2s 1 3d 1 electron configuration, is located 0.195 a.u. above the 1 1 P level, and hence, the 1 1 P → 7 1 D transition is nearly resonant with the first laser. Consequently, the 7 1 D level population increases once sufficient population of the 1 1 P level is achieved toward the end of the first pulse. The second pulse induces a single-cycle Rabi oscillation between the 1 1 P and 2 1 S levels (f = 0.118), which is quite significantly perturbed by the transition between the 2 1 S and 4 1 P levels (f = 0.211). The population of the 4 1 P level drops to zero as the Rabi oscillation enters the final stage where the population of the 2 1 S level decreases.
The CCLR populations are in close agreement with TDFCI populations and with EOMCC populations: the rms deviation for the entire simulation is 7.4 × 10 −3 between the CCLR and TDFCI populations and 9.1 × 10 −5 between the CCLR and EOMCC populations.
Increasing the basis set to aug-cc-pVTZ, the CCSD levels, the higher-lying ones in particular, move down in energy as seen by comparing the left panel of Figure 3 with that of Figure  2.
The right panel of Figure 3 shows the variation of the level populations as the Be atom is exposed to the same sinusoidal laser pulses as in Figure 2, albeit with the carrier frequencies adjusted to match the 0 1 S → 1 1 P and 1 1 P → 2 1 S transitions at ω 0 = 0.195 a.u. and ω 0 = 0.0539 a.u., respectively. The duration is 10 optical cycles for each pulse, as above. The populations obtained from the CCLR and EOMCC projectors are virtually identical, with an overall rms deviation of 1.6 × 10 −4 .
The lowering of the 7 1 D level, which is now 0.162 a.u. above the 1 1 P level, implies that the probability of the 1 1 P → 7 1 D transition (f = 0.648) diminishes, resulting in a very low population of the 7 1 D level an and increased population of the Figure 2. Energy-level populations computed with the aug-cc-pVDZ basis set through CCLR projectors plotted as functions of time for Be exposed to two consecutive laser pulses, the first resonant with the 0 1 S → 1 1 P transition and the second resonant with the 1 1 P → 2 1 S transition, with peak intensity 0.877 TW/cm 2 . Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article 1 1 P level (compared with the aug-cc-pVDZ simulation) during the first pulse. Although the populations of the states involved thus are different, the perturbed Rabi oscillation induced by the second pulse is essentially the same as in Figure 2. While a full TDFCI simulation is too costly with the aug-cc-pVTZ basis set, we compare with the much simpler model introduced for He above. We use the same consecutive laser pulses and the Gauss−Legendre integrator with the same parameters as in the Be/aug-cc-pVTZ TDCCSD simulation for the model simulation with 22 states included.
The resulting energy-level populations, plotted as dotted lines in Figure 3, are remarkably similar to the full TDCCSD results. The maximum absolute deviations between the model and TDCCSD populations are just 15% for the 1 1 P level, 2% for the 2 1 S and 7 1 D levels, and below 1% for the remaining levels, including the ground state. Such good results can only be expected from the simple model when all, or very nearly all, participating CCSD states are included. How many states are needed will in general be very hard to estimate a priori.

Control by Chirped Laser Pulses.
Control by shaped laser pulses is an important challenge to theoretical simulations and requires information about population of energy levels. For the LiH molecule described with the aug-cc-pVDZ basis set, which is small enough to allow TDFCI reference simulations, we use chirped sinusoidal laser pulses to further test the proposed CCLR and EOMCC projectors within TDCCSD simulations. The laser pulses are polarized along the x-axis, perpendicular to the molecular axis, and the duration is kept fixed at t d = 378.4 a.u. = 9.152 fs, corresponding to 10 optical cycles of radiation resonant with the lowest-lying electric-dipole allowed transition from the ground state, the x-polarized 1 Σ + → 1 Π transition at ω 0 = 0.166 a.u. The laser pulses are turned on at t 0 = 0 a.u. with peak intensity 3.51 TW/cm 2 (ε 0 = 0.01 a.u.). The oscillator strength of this transition is estimated to be f = 0.208 by CCSD linear response theory with the aug-cc-pVDZ basis set. The phase of the laser pulse is defined such that the instantaneous frequency is 78) and the chirp rate b is varied between −0.513 and +0.513 fs −2 .
A few such laser pulses with different chirp rates are shown in Figure 4. The 31 lowest-lying states, organized into 21 energy levels in the left panel of Figure 5, were computed with CCSD linear response theory (aug-cc-pVDZ basis) and used to construct CCLR and EOMCC projectors for the simulations. The highest-lying energy level is 0.603 a.u. above the ground state, well beyond the ionization energy estimated by CCSD/ aug-cc-pVDZ total-energy difference to be 0.281 a.u. 83 The most important energy levels in the simulations are marked by their term symbols. These states are all predominantly singleexcited states, with at least 90% contribution from singles in the EOMCC excitation amplitudes. While the 1 Π level at 0.166 a.u. is well below the estimated ionization energy, the 1 Δ and 1 Σ + levels at 0.291 and 0.312 a.u., respectively, are slightly above. With x-polarized laser pulses, one-photon transitions from the ground state to the 1 Δ and 1 Σ + levels are electricdipole forbidden, implying that these excited levels can only become populated by nonlinear optical processes.
The final populations of these levels, computed immediately after the interaction with the chirped sinusoidal laser pulse, are shown in the right panel of Figure 5 along with a few reference TDFCI results. The TDCCSD (and TDFCI) equations of motion were integrated using the sixth-order (s = 3) Gauss-Legendre integrator with time step Δt = 0.1 a.u. and convergence threshold 10 −6 . The sum of the populations of the remaining 27 energy levels is labelled "Rest" and is seen to be insignificant for all but the most up-or down-chirped pulses. At b = 0 fs −2 , the pulse is resonant with the ground-state 1 Σ + → 1 Π transition and, therefore, other levels are hardly populated. The maximum population of the 1 Π level is observed at a slightly up-chirped pulse (at b = 0.023 fs −2 ), which prevents further excitation from the 1 Π level to the higher-lying 1 Σ + and 1 Δ levels. As the chirp rate increases, the laser pulse becomes increasingly off-resonant, and the groundstate population increases.
The population of the excited 1 Σ + level and, in particular, of the 1 Δ level increases with moderately down-chirped pulses because of transitions from the 1 Π level, whose probability increases as the laser frequency decreases. At b = −0.102 fs −2 , the population of the 1 Π level is just 9.3 × 10 −4 , while the excited 1 Σ + and 1 Δ populations are close to their maximum values. These nonlinear optical processes can easily be understood by studying the populations during interaction with the laser pulse in Figure 6.
During the first half of the pulse, the instantaneous laser frequency is nearly resonant with the ground-state 1 Σ + → 1 Π transition at 0.166 a.u. (f = 0.207), causing the population of  Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article the 1 Π level to increase. As the instantaneous frequency decreases, it comes closer to the transition frequencies of the 1 Π → 1 Δ and 1 Π → 1 Σ + transitions at 0.138 and 0.146 a.u., respectively. Although the instantaneous frequency is closer to the 1 Π → 1 Σ + transition than to the 1 Π → 1 Δ transition, the latter level is considerably more populated. This is explained by the difference in the oscillator strength for these transitions: f = 0.363 for the 1 Π → 1 Δ transition compared with f = 0.076 for the 1 Π → 1 Σ + transition. The differences between the EOMCC and CCLR projectors are, again, utterly insignificant, the typical rms population deviation between them being approximately 10 −5 regardless of the chirp rate. As can be seen from Figure 5, the TDCCSD simulations are in excellent agreement with TDFCI results. This is not unexpected since all states participating in the dynamics are single-excited states. The maximum absolute deviation in the 30 excitation energies, including the excited states with significant double-excited character (down to 11% singles contribution in the EOMCC excitation amplitudes), between CCSD linear response theory and FCI theory is just 0.0013 a.u.
4.3. Dynamics Involving Double-Excited States. In general, CCSD linear response theory performs poorly for states dominated by double-excited determinants relative to the HF ground state. For such states, excitation-energy errors are typically an order of magnitude greater than for singleexcited states, 86,87 roughly 0.01 a.u. for double-excited states compared with 0.001 a.u. for single-excited states of small molecules where FCI results are available. 87 In all examples presented above, the states participating significantly in the dynamics are all single-excitation dominated, explaining the close agreement observed between TDFCI and TDCCSD simulations.
With a ground-state wave function dominated by the HF ground-state determinant, one-photon transitions to excited states dominated by double-excited determinants are either electric-dipole forbidden or only weakly allowed. Accordingly, we expect double-excited states to influence laser-driven manyelectron dynamics mainly through nonlinear optical processes. In order to test the influence on TDCCSD dynamics, we consider the CH + molecule, which is a classic example of the relatively poor performance of CCSD linear response and EOMCCSD theory for such states, see, for example, refs. 88−90 The 1 Σ + ground state of CH + is dominated by the 1σ 2 2σ 2 3σ 2 electron configuration with some nondynamical correlation contribution from the double-excited 1σ 2 2σ 2 1π 2 configuration.
The two lowest-lying excited states form the 1 Π energy level and are dominated by the single-excited 1σ 2 2σ 2 3σ 1 1π 1 configuration. The three subsequent states form two energy levels, 1 Δ and 1 Σ + , and are almost purely double-excited states stemming from the 1σ 2 2σ 2 1π 2 electron configuration. Transitions from the ground state to these levels are either electricdipole forbidden ( 1 Δ) or very weak with oscillator strengths on the order of 10 −3 . Significant population of these levels, therefore, can only be achieved through nonlinear optical processes, requiring rather intense laser pulses.
In order to make TDFCI simulations feasible for CH + , we use a reduced aug-cc-pVDZ basis set where the diffuse p functions on hydrogen and the diffuse d functions on carbon have been removed. While removing these diffuse functions has little effect on the 5 lowest CCSD linear response excitation energies (rms deviation 0.001 a.u.), the effect is significant on the following 25 excitation energies with an rms deviation of 0.021 a.u. The 31 lowest-lying states, forming 21 energy levels, computed with CCSD linear response and FCI theory are shown in the left panel of Figure 7.
The ionization energy of CH + is estimated to be 0.878 a.u. using ionization-potential EOMCCSD theory and 0.876 a.u. using FCI theory with the 6-31G* basis set. 91 Hence, all computed states are below the estimated ionization energies.
We expose the CH + ion to a sinusoidal laser pulse with intensity 2654 TW/cm 2 (ε 0 = 0.275 a.u.) and carrier frequency ω 0 = 0.212 a.u., which is resonant with the 1 Π → 1 Σ + transition between excited states at the CCSD level of theory. The pulse is polarized along the y-axis, perpendicular to the bond axis, with duration t d = 66.7 a.u. = 1.61 fs, corresponding to 2.05 optical cycles. The TDCCSD and TDFCI equations of motion were integrated with the sixth-order (s = 3) Gauss-Legendre integrator with time step Δt = 0.05 a.u. and convergence threshold 10 −6 . The resulting energy-level populations are presented in the right panel of Figure 7. Populations of the 4 lowest-lying levels, including the ground state, are shown along with the sum of the populations of the remaining 17 levels, labelled "Rest." The total population of all computed levels is labelled "Sum." The effect of poorly described double-excited states at the CCSD level of theory is evident, although we do observe a qualitative agreement with FCI theory. We first note that the 21 levels included in the analysis only account for about 80% of the norm of the FCI wave function at the end of the pulse, implying that a physically correct description must also take  Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article ionization processes into account. This, of course, is not surprising, considering the high intensity of the pulse. The TDFCI and TDCCSD populations agree reasonably well during the first 0.75 fs of the simulation, whereas the TDCCSD errors increase as the simulation progresses. Involving a large number of states, the dynamics is considerably more complex than the dynamics of the cases presented above. For simplicity, we use a classification of the 30 excited states based on their single-or double-excitation character. Excited states with more than 90% contribution from singles in the EOMCCSD amplitude norm are classified as single-excited states, while states with less than 10% singles contribution are classified as double-excited states. States with 10−90% singles contribution are mixed states, classified as singles-dominated (>50% singles contribution) or doublesdominated (<50% singles contribution). Thus, the 30 excited CCSD states can be grouped into 7 single-excited states, 4 singles-dominated states, 5 doubles-dominated states, and 14 double-excited states. The total population of each class of states is presented in Figure 8.
After about 0.4 fs, the laser pulse induces transitions from the ground state into singles-dominated states, followed by transitions (from both the ground state and the singlesdominated excited states) into single-excitation states. Doubleexcited or doubles-dominated states are barely populated at this stage, explaining the reasonable agreement with TDFCI populations. Roughly half-way through the simulation, doubleexcited states and, to a smaller extent, doubles-dominated states become populated, mainly due to transitions from singles-dominated states. As soon as these processes occur, the agreement between TDCCSD and TDFCI deteriorates.
4.4. Population Conservation. As discussed in Section 2.6, stationary-state populations in the absence of external forces are strictly conserved in the FCI limit but may vary when the cluster operators are truncated. In order to investigate the breaking of the population conservation law within TDCCSD theory, we have conducted simple numerical experiments with several of the systems presented above. We apply the same sinusoidal laser pulses as above but continue the propagation after the pulses have been turned off, recording stationary-state population using the CCLR and EOMCC projectors. Figure 9 shows the conservation of TDCCSD populations after the laser pulses have been turned off for the He atom with the aug-cc-pVTZ basis set.
A maximum absolute deviation in the populations of 1.5 × 10 −3 is observed for the 1 1 S and 2 1 P levels, whereas the deviations for the remaining levels are at least 1 order of magnitude smaller. These deviations from exact conservation are likely caused by the discretization of the numerical integration.
Slightly larger deviations from strict conservation are observed for the Be atom with the aug-cc-pVTZ basis set in Figure 10.
The maximum absolute deviation of 0.005 is observed for the 1 1 P level. Caused by the truncation of the cluster operators in conjunction with discretization, this deviation is thrice greater than that observed for He above. Only weak oscillatory behavior is observed, indicating that the states involved in the dynamics are very well approximated at the CCSD level of theory.
Energy-level populations for LiH during and after interaction with a chirped laser pulse are plotted in Figure 11.
The chirp rate b = 0.03078 fs −2 yields a final state dominated by the 1 Π level, and the populations remain constant to an excellent approximation after the interaction ceases. The    Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article maximum absolute deviation is observed for the 1 Π level and is on par with that observed for the He atom: 1.7 × 10 −3 . As we saw above, the TDCCSD method is a much poorer approximation to TDFCI theory in the case of CH + , where double-excited states participate significantly in the dynamics. On these grounds, we expect a much poorer conservation of energy-level populations after interaction with the laser pulse and, indeed, large deviations can be seen in Figure 12.
The depletion of the ground state appears to continue after the pulse, and irregular weakly oscillatory behavior is observed for the 1 Π level with a maximum absolute deviation of 0.072, which is an order magnitude greater than the deviations found for He, Be, and LiH above. The double-excited 1 Δ and 1 Σ levels show maximum absolute deviations of 0.009 and 0.019, respectively.
4.5. Pump Spectrum of LiF. Skeidsvoll et al. recently reported a theoretical TDCCSD study of transient X-ray spectroscopy of the LiF molecule. 38 They applied a pump− probe laser setup with an optical pump pulse resonant with the lowest-lying dipole-allowed transition from the ground state, followed, at various delays, by an X-ray probe pulse resonant with the first dipole-allowed core excitation. The resulting time-resolved spectra were interpreted by means of excitation energies from EOMCC theory and core−valence separated EOMCC theory. 92 The pump absorption spectrum reported in Figure 7 of ref 38 contains weak unassigned features, one weak absorption above the two low-lying valence excitations and two very weak features below, which the authors speculated were due to two-photon absorptions. We will now use the EOMCC and CCLR projectors to investigate what might cause these weak features of the pump absorption spectrum. We use the same basis set, denoted aug-cc-p(C)VDZ, as in ref 38: the augcc-pVDZ basis set for Li and the aug-cc-pCVDZ basis set for F. The closed-shell TDCCSD equations of motion were integrated using the sixth-order (s = 3) Gauss−Legendre integrator with time step Δt = 0.025 a.u. = 0.60 as and convergence threshold 10 −6 for the fix-point iterations.
Initially in the ground state, we expose the LiF molecule to a shortened but otherwise identical z-polarized Gaussian laser pulse to that in ref 38, with field strength ε 0 = 0.01 a.u. (peak intensity I = 3.51 TW/cm 2 ), carrier frequency ω 0 = 0.2536 a.u., and Gaussian rms width σ = 20 a.u. The shortening consists in choosing the central time t 0 = 80 a.u (compared with t 0 = 160 a.u. in ref 38) and N = 4 (compared with N = 8 in ref 38). This implies that the electric-field amplitude jumps from zero to 3.3 × 10 −6 a.u. at t = 0 a.u. and from 3.3 × 10 −6 a.u. to zero at t = 160 a.u., whereas virtually no discontinuities can be observed with the pulse parameters used in ref 38 (they are on the order of 10 −16 a.u.). Since the pump pulse is quite weak, the effects of these discontinuities on the populations are negligible, as can readily be verified using a simple eigenstate expansion analogous to the one used for He and Be in Section 4.1.
Our TDCCSD results are presented in Figure 13. The first 30 excited states (20 energy levels, left panel of Figure 13) are all single-excitation dominated states, with 94.9−95.4% singles contribution to the norm of the EOMCCSD amplitudes. The highest-lying states are somewhat above the first ionization energy of LiF, which we estimate to be about 0.4 a.u. (11 eV) based on data available in ref 83.
The modest intensity of the pump pulse results in fairly little excitation from the ground state, well within reach of a perturbation-theoretical (Fermi's golden rule) treatment. In agreement with ref 38, the projectors predict the absorption to be dominated by the B 1 Σ + and E 1 Σ + states. The final population of the latter is roughly 53% of the former, in good agreement with the relative intensities of the pump spectrum reported in Figure 7 of ref 38. The population of the E 1 Σ + state reaches its maximum value 0.0025 at t = 2.26 fs; at the same time, the ground-state population reaches its minimum value 0.9945, indicating that the ensuing decay of the E 1 Σ + population is caused by transition back to the ground state.
The weak feature at higher frequency (at about 9.1 eV) in the pump spectrum of ref 38 is seen to be consistent with onephoton transition from the ground state to the H 1 Σ + state, whose final population is about 5% of that of the B 1 Σ + state. As speculated by Skeidsvoll et al., 38 the two very weak features below the B 1 Σ + line in ref 38 are indeed seen to arise from direct two-photon absorptions from the ground state to the L 1 Σ + and R 1 Σ + states. The only alternative explanation would be excitations between excited states, but this mechanism can almost certainly be ruled out since the population of these states starts before other excited levels are significantly populated and since no other excited states are depleted as the populations of these states increase.
Although the CCLR and EOMCC populations largely agree with an overall rms deviation of 7 × 10 −6 , the CCLR populations show spurious high-frequency oscillations in Figure 13. The oscillations are caused by the off-diagonal contributions from ⟨Ψ̑n| (eq 49) to the CCLR projector (eq 47). While these contributions are required to ensure proper  Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article size-intensivity of one-photon transition moments, they also cause nonorthogonality of the CCLR excited-state representation as expressed by eq 57. Since the CCLR and EOMCC projectors provided virtually identical results in the cases above, this observation serves as a recommendation of the EOMCC projector for the calculation of stationary-state populations. Figure 14 depicts a normalized pump spectrum of LiF generated from the final EOMCC populations using where γ = 0.01 eV is an artifical Lorentzian broadening of the excited levels. This approach implicitly assumes that the excited states only become populated through one-photon absorption from the ground state, thus excluding all nonlinear optical processes. The population-based pump spectrum agrees remarkably well with that reported by Skeidsvoll et al., 38 which was properly generated from Fourier transformation of the induced dipole moment. This supports the conclusion that the low-frequency features are two-photon absorptions and further strengthens the confidence in the proposed EOMCC projector for the calculation of stationary-state populations in TDCC simulations.

CONCLUDING REMARKS
We have proposed projectors for the interpretation of manyelectron dynamics within TDCC theory in terms of the population of stationary states. Two conditions are used to define suitable projectors from CCLR theory and from EOMCC theory: (i) the projector must reproduce the correct form of one-photon transition strengths and (ii) the projectors must yield populations that converge to the FCI results in the limit of untruncated cluster operators.
The CCLR and EOMCC projectors are tested numerically at the TDCCSD level of theory for the laser-driven dynamics in the He and Be atoms and in the molecules LiH, CH + , and LiF. It is demonstrated that the populations provide valuable insight into the linear and nonlinear optical processes occurring during the interaction of the electrons with laser pulses. For the He atom, it is verified numerically that the populations computed from both CCLR and EOMCC projectors agree with those computed from TDFCI simulations. For Be and LiH, the CCLR and EOMCC populations show excellent agreement with TDFCI popula-tions since the excited stationary states involved in the dynamics are dominated by single-excited Slater determinants. Such states are generally well described by CCSD theory. For CH + , we deliberately design the laser pulse such that doubleexcitation dominated states become populated, which reduces the agreement between TDFCI and TDCCSD populations. This is also reflected in the studies of the conservation of populations after the laser pulses are turned off. Theoretically, the TDCC populations will only be strictly conserved in the FCI limit. Numerically, we find that they are nearly conserved as long as the participating stationary states are well approximated at the CCSD level of theory. Finally, for LiF, we use the CCLR and EOMCC projectors to explain unassigned weak features in a theoretical TDCCSD pump spectrum reported recently. 38 Overall, the CCLR and EOMCC projectors yield very similar excited-state populations with typical rms deviations on the order of 10 −5 . For LiF, however, we observe smallamplitude high-frequency oscillations of the excited-state populations computed with the CCLR projector. Originating from a contribution that vanishes in the FCI limit, we speculate that such oscillations may increase for larger and more complex systems where TDCCSD theory may be further from TDFCI theory. Not showing signs of such spurious behavior, the EOMCC projector appears more attractive than the CCLR projector. This has the added benefit that the additional response eq 51 need not be solved, thus making the EOMCC projectors less computationally demanding than the CCLR projectors.
These findings call for further research aimed at a fully consistent definition of excited states in CC theory, and work in this direction is in progress in our labs.

■ APPENDIX FCI Limit
The CCLR projector, eq 47, becomes identical to the EOMCC projector, eq 33, in the FCI limit. To show this, we will now demonstrate that in the FCI limit.
In the FCI limit, the EOMCC wave functions satisfy the time-independent Schrodinger equation and its Hermitian conjugate where the ground-and excited-state wave functions constitute a biorthonormal set according to eqs 39 and 40. The resolution-of-the-identity reads where 1 is to be understood as the identity operator. To verify the resolution-of-the-identity, one simply inserts the definitions of the wave functions and exploits the biorthonormality of the Jacobian eigenvectors, eq 37, along with the completeness of the underlying determinant basis. According to eq 51, the amplitudes n ̅ μ can be recast as where we have used eq 37. Expanding the nested commutator and using 81−83 and, 39 we find