Time-dependent coupled-cluster theory

Recent years have witnessed an increasing interest in time-dependent coupled-cluster (TDCC) theory for simulating laser-driven electronic dynamics in atoms and molecules, and for simulating molecular vibrational dynamics. Starting from the time-dependent bivariational principle, we review different flavors of single-reference TDCC theory with either orthonormal static, orthonormal time-dependent, or biorthonormal time-dependent spin orbitals. The time-dependent extension of equation-of-motion coupled-cluster theory is also discussed, along with the applications of TDCC methods to the calculation of linear absorption spectra, linear and low-order nonlinear response functions, highly nonlinear high harmonic generation spectra and ionization dynamics. In addition, the role of TDCC theory in finite-temperature many-body quantum mechanics is briefly described along with a few other application areas. This article is categorized under:


| INTRODUCTION
The objective of time-dependent molecular electronic structure theory 1,2 is to solve the time-dependent Schrödinger equation (TDSE), where the Hamiltonian, b , contains the time-independent, clamped-nuclei Born-Oppenheimer electronic Hamiltonian, 3,4 b H 0 , and an explicitly time-dependent operator, b V t ð Þ, representing the interaction of the electrons with external driving forces, usually electromagnetic fields.Atomic units will be used throughout this review unless explicitly stated otherwise.With the appropriate initial state, j Ψ 0 i, the TDSE thus allows simulation of electronic dynamics directly corresponding to experimental setups.
However, within the clamped-nuclei approximation, time-dependent electronic structure theory can only be expected to yield reliable dynamics on time scales short enough that nuclear motion can be neglected-typically, up to a few femtoseconds.An important example is the charge migration across the nuclear framework initiated by ionization. 5For longer time scales, the nonadiabatic coupling of electronic and nuclear motion must be taken into account using, for example, semiclassical methods based on the Ehrenfest theorem 6,7 or, if feasible, fully quantum-mechanical methods like the time-dependent Feshbach close-coupling method. 8An obvious application of time-dependent electronic structure theory thus is the simulation of processes induced by attosecond laser pulses, 9 including attosecond transient absorption spectroscopy. 10However, the broad spectral range of attosecond laser pulses combined with relatively high intensities almost invariably induce ionization processes, placing heavy demands on the electronic structure method and basis sets, which must accurately capture both bound states and the electronic continuum.The most accurate approach is to solve the TDSE using mesh-based methods in both space and time such as, for example, the finite element discrete variable representation [FEDVR]. 11Although not without limitations, a less computationally demanding approach is based on density-functional theory in combination with scattering states expanded in B-splines. 9,12irect simulation of experiments is not the only valuable application of the TDSE; for example, energies and in principle also the associated stationary-state wave functions of b H 0 can be extracted from simulations without external driving forces simply by starting in a nonstationary state. 13With judicious but artificial choices of the interaction operator b V t ð Þ, the TDSE can also be used as an alternative to perturbation theory for the calculation of molecular optical properties and spectra.A weak electric-field kick applied to the electronic ground state, for example, yields linear absorption spectra from the induced electric-dipole moment.The spectra automatically include all electric-dipole allowed transitions, both valence and core excitations, from a single simulation or a few simulations, depending on symmetry and whether the spectrum is simulated for an aligned or randomly oriented sample. 2 This approach avoids diagonalization of large matrices and is potentially advantageous for systems with a high density of states where a large number of eigenvalues would be needed in traditional time-independent methods.Linear and low-order nonlinear optical propertiespolarizabilities and hyperpolarizabilities-can be extracted from induced electric and magnetic moments using either a ramped continuous wave 14 or a pulsed wave 15 to perturb the ground-state wave function.Only a few, relatively short simulations are required.
The main computational obstacles for methods aimed at solving the TDSE are the long simulation times required to achieve sufficient resolution in Fourier analyses of the recorded signals and the small time steps required to capture high-frequency components in the wave function.Since electron correlation effects must be accounted for, it is no surprise that the most widely used electronic-structure method for electronic dynamics simulations is time-dependent density-functional theory (TDDFT), [16][17][18] often called real-time (RT) TDDFT to clearly distinguish it from perturbationbased density response theory in the frequency domain.Similar to ground-state calculations, the TDDFT approach often strikes a reasonable balance between computational effort and accuracy.
For higher accuracy, one must turn to methods that parameterize the wave function explicitly.Since electronic excited states are often multi-configurational, high-accuracy simulations of electronic dynamics have been dominated by the multi-configurational time-dependent Hartree-Fock (MCTDHF) method [19][20][21][22][23] or the closely related complete, 24 restricted, 25,26 and generalized 27 active space self-consistent field methods.These methods suffer from the curse of dimensionality, making coupled-cluster (CC) [28][29][30][31] approximations attractive alternatives with their more benign polynomial scaling, at least for simulating electronic processes where the time-dependent wave function is dominated by a single, generally time-dependent, electron configuration.
Time-dependent CC (TDCC) theory was first formulated by Monkhorst 32 in 1977, albeit not with the purpose of studying electronic dynamics.Rather, Monkhorst 32 and later Dalgaard and Monkhorst 33 applied perturbation theory to the TDCC equations and derived expressions for linear response properties such as the frequency-dependent electricdipole polarizability, identifying excitation energies-the poles of the linear response function-as the eigenvalues of a non-hermitian matrix.The non-hermiticity implies that the eigenvalues are not necessarily real and, indeed, Takahashi and Paldus 34 observed complex eigenvalues in their orthogonally spin-adapted TDCC approach to excitation energies.However, they only obtained complex excitation energies with a (quoting Takahashi and Paldus 34 ) "very poor (in fact almost meaningless) CC representation for the ground state" of strongly correlated systems.5][46] While Monkhorst's TDCC formulation was based on a fixed reference determinant, Hoodbhoy and Negele 44,45 allowed the underlying spin orbitals to be time-dependent, determined by time-dependent Hartree-Fock theory.(We will use the terms spin orbital and orbital interchangeably.)Note, however, that Pigg et al. 46 used static orbitals in the most recent application of TDCC theory to nucleon dynamics.8][49] Shortly after the publication of Hoodbhoy and Negele's first article, 44 Schönhammer and Gunnarsson 50 applied TDCC theory to compute the spectral weight function from the phase factor of the TDCC wave function for prediction of core-level spectra of atomic and molecular adsorbates.The next application of TDCC theory was also to a solid-state problem.In 1985, Sebastian 51 used TDCC theory to simulate scattering of high-energy cations from surfaces, computing the probability of neutralization through a one-electron charge-transfer process.Sebastian proposed a conventional expression for the TDCC expectationvalue functional which, for truncated cluster operators, does not fulfill the Hellmann-Feynman theorem. 52,53 CC expectation-value functional that fulfills the Hellmann-Feynman theorem is a key result of Arponen's bivariational formulation 54 and of the equivalent constrained optimization approach of Helgaker and Jørgensen.55,56 This expectation-value functional quickly became the standard choice and was used by Bishop and Emary in their TDCC study of a two-level system in a quantized electromagnetic field.35 The remainder of this review is organized as follows.We start in Section 2 with a brief summary of traditional perturbative approaches to TDCC theory with emphasis on concepts that are important for the time-dependent bivariational theory presented in Section 3. In Section 4, we describe the time-dependent extension of equation-ofmotion CC theory, which is equivalent to TDCC theory in the limit of untruncated cluster operators, while Section 5 reviews TDCC theory applied to molecular vibrational dynamics.Finally, Section 6 briefly reviews a few other application areas where TDCC theory plays a role, and Section 7 contains our concluding remarks.

| PERTURBATION-BASED APPROACHES
Besides the early applications to dynamical systems mentioned above, TDCC theory has mainly been used as a starting point for the calculation of frequency-dependent response properties in molecular electronic-structure theory. 57nspired by Helgaker and Jørgensen's 55,56 constrained optimization (Lagrangian) approach to static properties as energy derivatives and by Olsen and Jørgensen's 58 time-dependent variational formulation of response theory, Koch and Jørgensen 59 generalized the perturbation-based TDCC approach of Monkhorst 32 and Dalgaard and Monkhorst 33 to formulate a general CC response theory.Equivalent to the bivariational theory of the "normal" TDCC method of Arponen, 54 the starting point is independent Ansätze for the wave function and its conjugate where the normalized reference Slater determinant, j Φ 0 i, is time-independent and usually taken to be the Hartree-Fock ground-state determinant.The cluster operators b where the summations are over all possible excitations out of the reference determinant: b X and b Y denote excitation and de-excitation operators such that Note that the excitation and de-excitation operators separately commute: b While the phase amplitude, τ 0 t ð Þ, plays no role in CC response theory as formulated by Koch and Jørgensen, it becomes important in the interpretation of TDCC dynamics. 60,61runcation of the cluster operators after single excitations yields the TDCC singles (TDCCS) model, after singles and doubles yields the TDCC singles and doubles (TDCCSD) model, and so on.If the cluster operators are not truncated, TDCC theory becomes equivalent to the formally exact time-dependent full configuration interaction (TDFCI) theory, albeit with a wave function which is not normalized.Regardless of truncation, the TDCC wave functions instead satisfy the (intermediate) normalization conditions The equations of motion for the amplitudes are derived by inserting j Ψ t ð Þi and h e Ψ t ð Þ j into the TDSE and its conjugate, respectively, followed by projection onto the excited determinants, yielding 59 i_ where, as usual, the "dot" denotes the time derivative.The phase amplitude is determined by Assuming that the electronic system is initially in its ground state and that the driving forces are adiabatically switched-on and weak, the amplitudes are expanded in orders of the perturbation, leading to equations that must be solved order-by-order.Rather than solving the equations in the time domain, it is assumed that the Fourier transforms of the interaction operator, b V t ð Þ, and of the nth order, n ≥ 1 amplitudes exist, such that the amplitude equations can be transformed to the frequency domain and solved at the frequencies of interest.
The expectation value of an operator b A takes the form originally suggested by Arponen 54 b A 3][64] Note that the expectation value does not depend on the phase amplitude τ 0 t ð Þ. Expanding the expectation value in orders of the perturbation then leads to an identification of linear, quadratic, cubic, etc. response functions.If, for example, b V t ð Þ represents the interaction of the electrons with an electric field in the electric-dipole approximation, the response functions are (apart from a sign) the electricdipole polarizability and hyperpolarizabilities.Excitation energies and transition moments are identified from the poles and residues of the response functions.If the interaction operator is further assumed periodic in time such that its Fourier series converges, equivalent expressions for the CC response functions can be conveniently derived as derivatives of the cycle-averaged quasienergy Lagrangian. 57,65olving the TDSE with the ground state as the initial condition is thus the starting point for CC response theory.In principle, the perturbative solutions obtained in the frequency domain can be transformed back to the time domain but this would require a dense grid of frequencies and is never done in practice.Hence, in effect, CC response theory is a time-independent method.Another time-independent approach is equation-of-motion CC (EOM-CC) theory, 30,[66][67][68][69][70] where excited states are written explicitly as linear excitation operators acting on the CC ground-state wave function.This is a fundamentally different approach from CC response theory where explicit expressions for excited-state wave functions are neither needed nor assumed.Still, the excitation energies obtained from CC response theory are identical to those obtained from EOM-CC theory: The same non-hermitian eigenvalue problem arises in both theories.Transition moments and response functions differ, however, and only those obtained from (truncated) CC response theory are properly size consistent. 71,72Both CC response theory and EOM-CC theory converge to the full configuration interaction (FCI) limit when no truncation of the cluster (and linear excitation) operators are introduced.A formulation of response theory starting from a time-dependent EOM-CC Ansatz has been given by Coriani et al. 73 For truncated cluster operators, the expectation-value functional (10) leads to broken symmetry properties of the CC response functions under complex conjugation.This issue, which results in spurious origin-dependence of some optical properties, is rather easily fixed by using only the real part of Equation (10), as suggested by Pedersen and Koch. 74Alternatively, one can enforce the proper symmetries on the response functions a posteriori. 655][76][77] In order to resolve this, the orbitals must be dynamical variables and Pedersen, Koch and coworkers proposed the time-dependent orbital-optimized CC (TDOCC) model 47 and later the time-dependent nonorthogonal orbitaloptimized CC (TDNOCC) model, 48 which were formulated as response theories.Never implemented in a production-level code, the TDOCC and TDNOCC models have not been extensively used or tested, but they have gained importance for studying laser-driven electronic dynamics through Kvaal's formulation of orbital-adaptive time-dependent coupled-cluster (OATDCC) theory, 49 which is an adaptation of the basic idea of MCTDHF theory.The OATDCC model is equivalent to TDNOCC theory if the underlying time-dependent orbital space is not split into active and external subspaces.For dynamics, it is advantageous to start from the bivariational approach of Arponen. 54

| ELECTRONIC DYNAMICS WITH BIVARIATIONAL CC THEORIES
While computationally demanding, RT TDCC theory offers clear advantages over perturbation-based approaches.Explicitly time-dependent simulations contain responses to all orders and, therefore, are able to describe highly nonlinear optical phenomena in a time-resolved manner.Furthermore, experimental parameters like pulse shape and pulse duration can be embedded directly into the simulation.The TDCC methods can accordingly serve as a theoretical complement to the increasingly topical field of experimental attosecond science.For certain use cases, such as the calculation of near-edge x-ray absorption spectra, TDCC simulations may even prove to be computationally competitive.
The first application of TDCC theory to simulate laser-driven molecular electronic dynamics was presented in 2011 by Huber and Klamroth 78 who used the semiclassical electric-dipole approximation, truncated the cluster operator b T t ð Þ after double excitations to obtain the TDCCSD model, and propagated only the τ amplitudes according to Equation (7).Consequently, induced electric-dipole moments could not be computed from Equation (10) and Huber and Klamroth instead resorted to an approximation using the configuration interaction singles-and-doubles expression, leading to relatively large errors in excitation energies.More detrimental to the prospects of TDCC theory of laser-driven electronic dynamics, Huber and Klamroth found that the TDCCSD method became numerically unstable in strong external fields and with increasing basis set quality.

| The bivariational formalism
Arponen's bivariational principle 54 naturally leads to the expectation-value functional (10) which satisfies the Hellmann-Feynman theorem, both in the time-independent and in the time-dependent case.In the FCI limit, Arponen's extended CC formulation is equivalent to the Lagrangian approach of Helgaker and Jørgensen, 55,56 as discussed in detail by Kvaal. 79he starting point is the bivariational action functional 54 S e Ψ, Ψ with the Lagrangian where the Hamilton function is The ket and bra, j Ψ t ð Þi and h e Ψ t ð Þ j, are independent approximations to the exact wave function and its conjugate, respectively.Under the condition that the bivariational action functional is complex analytic, requiring S be stationary with respect to variations in the bra and in the ket leads to the TDSE and its conjugate: These equations guarantee that the normalization condition is conserved.Evidently, Equations ( 14) and ( 15) are complex symplectic generalizations of the classical Hamiltonian equations 80 with the ket and bra functions as canonical variables.This connection with classical Hamiltonian mechanics has been extensively explored by Arponen and coworkers [81][82][83][84] (see also Ref. [64]) and suggests that j Ψ t ð Þi and h e Ψ t ð Þ j together represent the quantum state of the system.This observation led Pedersen and Kvaal 60 to propose the two-component state vector and the indefinite inner product This form still satisfies the Hellmann-Feynman theorem and naturally leads to the correct symmetries in response functions. 64,74For hermitian operators, the expectation-value functional (19) equals the real part of Equation (10).
Using z to denote the vector of all wave function variables, the Lagrangian becomes a function of z, _ z, and time t, Þ, and the stationarity condition becomes the Euler-Lagrange equations which may offer a simpler derivation of the equations of motion for the chosen approximate parameterization.

| Formulations of TDCC theory with static or time-dependent orbitals
The Ansätze for the ket and bra states take the general form such that, with λ 0 ¼ 1, the state vector ( 17) is normalized with respect to the indefinite inner product (18), The phase amplitude τ 0 is canonically conjugate to λ 0 .In the conventional TDCC theory for an N-electron system, j Φ 0 t ð Þi ¼j Φ 0 i is chosen to be the time-independent field-free Hartree-Fock ground-state determinant and h e Φ 0 t ð Þ j¼ hΦ 0 j.The underlying spin orbitals are orthonormal and the cluster operator b ) contains from single to n-tuple, 1 ≤ n ≤ N, excitation (de-excitation) operators with respect to the Hartree-Fock determinant.The cluster operators are parameterized by the amplitudes τ t ð Þ and λ t ð Þ, one amplitude per excitation and de-excitation, as in Equation ( 4).If n ¼ N, conventional TDCC theory is equivalent to the formally exact TDFCI theory, although the intermediate normalization exp i¼ 1 may cause severe numerical instabilities (see below for details).The amplitude equations of motion are given by Equations ( 7) and (8).
In time-dependent nonorthogonal orbital-optimized CC (TDNOCC) theory, 48 all determinants contributing to j Ψ t ð Þi and h e Ψ t ð Þ j are time-dependent variational parameters.The underlying spin orbitals constitute a biorthonormal set, e Single excitations (de-excitations) are redundant when the orbitals are timedependent 48,49 and, consequently, they are removed from b ) in TDNOCC theory.Inspired by MCTDHF theory, orbital-adaptive TDCC (OATDCC) 49 theory adds the concept of active orbital space to TDNOCC theory.The cluster operators are restricted to a subset of the orbital space, which is optimized along with its orthogonal complement throughout the dynamics.This approach is very important for describing ionization dynamics.
With time-dependent orbitals, the amplitude equations of motion become 49 i_ Since _ λ 0 ¼ 0, normalization with respect to the indefinite inner product is conserved and λ 0 ¼ 1 is a natural choice.The time-dependence of the orbitals gives rise to the operator where the creation and annihilation operators, which satisfy the usual fermionic anticommutator relations, refer to the biorthonormal orbitals.With time-dependent orbitals, the τ and λ amplitude equations become coupled and must be solved simultaneously.If all orbitals are chosen active, we may write 49 The nonvanishing components of η p q are determined by the linear equations i Ài where the right-hand sides depend on the correlated one-and two-electron effective density matrices, thus coupling the orbital evolution to the correlating amplitudes.For explicit expressions for the right-hand sides, we refer to Equations (30a) and (30b) of Ref. [49].The matrix elements A ib aj are given by where ρ is the one-electron effective density matrix.In order for Equations ( 27) and ( 28) to be well-determined, the matrix A ¼ A ib aj h i must remain nonsingular at any time t.While this cannot be mathematically guaranteed, the singularity has not been reported in any publication to date.
Constraining the spin orbitals to be orthonormal throughout the dynamics, h e ϕ p t ð Þ j¼ hϕ p t ð Þ j, TDNOCC theory turns into time-dependent orbital-optimized CC (TDOCC) theory. 47,85The Lagrangian is forced to be real, ℒ Re ℒ ð Þ, and the two orbital equations of motion ( 26) are then related by complex conjugation, with the right-hand side of the orbital equation of motion given in Equation ( 23) of Ref. [85].As in OATDCC theory, the orbital space can be split into active and inactive subspaces in TDOCC theory, 85 facilitating simulations of highly nonlinear optical phenomena such as ionization.
The TDOCC, TDNOCC, and OATDCC theories thus are very closely related and, unlike TDCC theory based on static orbitals, they all provide gauge invariant results regardless of the truncation level of the cluster operators.As demonstrated by Köhn and Olsen, 86 however, TDOCC theory does not reproduce TDFCI results in the limit of untruncated cluster operators when N > 2. While this unfortunate feature was long believed to apply to TDNOCC theory, too, Myhre 87 recently showed that the correct limit can be obtained in TDNOCC theory for any N. On the other hand, the benchmark studies of Sato and coworkers 85,[88][89][90] indicate that the deviation of TDOCC results from TDFCI results is often negligible, at least for strong-field dynamics, see Section 3.5 below.

| Numerical integration
Collecting the bivariational parameters in a single vector z, the equations of motion can be written in the form of a complex ordinary differential equation (ODE), A wealth of numerical integrators for ODEs have been developed, addressing numerical issues such as conservation of symplectic structure and stiffness, see, for example, the authorative treatise by Hairer, Lubich, and Wanner. 91tarting with the work of Huber and Klamroth, 78 a popular integrator for TDCC theory has been the explicit fourthorder Runge-Kutta (RK4) algorithm.Its popularity can most likely be traced to two properties: It is very simple to implement and requires exactly four evaluations of the function f per time step, making computational time easily predictable.The RK4 integrator, however, is not symplectic and may thus break physically important conservation laws.Pedersen and Kvaal 60 instead proposed to use the symplectic s-stage Gauss-Legendre integrator 91 (which is of even order 2s, s ¼ 1,2,3, …) and showed that this yields long-time conservation of energy close to machine precision and that, depending on the time-step size and on the initial guess employed for the iterative solution of the implicit equations, may yield fewer f-evaluations per time step than the RK4 integrator for s ≤ 3.
More recently, Wang, Peyton, and Crawford 92 investigated modified Runge-Kutta integrators with adaptive time step, which increases stability when the parameters oscillate rapidly and allows larger time steps when they do not.Remarkably, stability and accuracy is maintained also in conjunction with single-precision arithmetic, which allows highly efficient calculations on graphical processing units.The combination of larger time steps when possible and single-precision arithmetic leads to significant acceleration (more than an order of magnitude).Note, however, that Wang, Peyton, and Crawford used the frozen core approximation, thus excluding the highest frequencies from the amplitude oscillations.The high energies associated with core excitations and with ionization dynamics, in particular, introduce stiffness in the TDCC equations and Sato et al. 85 proposed to use an exponential Runge-Kutta integrator to handle this.

| Bivariational interpretation
The definitions of the state vector (17) and of the indefinite inner product (18) provide the foundation for analysis of the electronic dynamics in close analogy with conventional quantum mechanics.Autocorrelation functions (ACFs)overlaps of the quantum state with itself at different times-contain important information about the dynamics.The early application of TDCC theory to core-level excitations of adsorbates by Schönhammer and Gunnarsson 50 is an interesting example that does not involve an external driving force.
Pedersen and Kvaal 60 defined the ACF as and demonstrated that the total energies of the stationary states contributing to the electronic dynamics can be extracted from it by Fourier transformation when t 0 is taken to be the switch-off time of an external laser pulse.It should be stressed that the phase amplitude τ 0 t ð Þ is important for a correct calculation of the ACF.The work of Pedersen and Kvaal 60 was restricted to conventional TDCC theory with static orbitals, since the calculation of overlaps between time-dependent Slater determinants exhibits factorial scaling, hampering the practical application of ACFs in TDCC theories with time-dependent orbitals.
Not only the energies of participating stationary states, but also their populations during the dynamics can be computed using the indefinite inner product.Pedersen et al. 61 introduced the (two-component) operator, projecting onto stationary state I. Here, j Ψ I i and h e Ψ I j are the right and left wave functions of state I from EOM-CC theory (see below for the precise definition).The projection operator (32) is hermitian with respect to the indefinite inner product (18) and the population of state I at time t thus becomes While inherently real, the populations are neither bounded below by 0 nor above by 1.Like ACFs, stationary-state populations are considerably more challenging to compute with time-dependent orbitals than with conventional static orbitals, and Pedersen et al. 61 only presented results for the latter.Pedersen et al. also proposed a projection operator based on CC linear response theory, which led to populations that are practically indistinguishable from those obtained with the EOM-CC projector in most cases.It was found, however, that the linear-response projector may lead to spurious high-frequency oscillations in the populations and it was recommended to mainly use the EOM-CC projector.
As an example, Figure 1 shows the final TDCCSD populations of stationary states below the ionization energy for the LiH molecule after interaction with a short, chirped laser pulse with carrier frequency resonant with the transition from the 1 Σ þ ground state to the lowest-lying electric-dipole allowed 1 Π state.The final populations are plotted as functions of the laser chirp rate b.It is interesting to note that the greatest population of the resonant 1 Π state is achieved by a slight up-chirp, whereas a slightly larger down-chirp leads to virtually no population of the same state.These effects are caused by transitions among excited states that are nonlinear optical processes from the viewpoint of response theory.As illustrated in Ref. [61], such nonlinear processes, including quenched Rabi oscillations between excited sates, can be tracked by recording populations.In the FCI limit, the populations are strictly conserved in the absence of external driving forces but may show slight drifts and low-amplitude oscillations with truncated cluster operators.The TDCCSD populations agree well with TDFCI populations provided all states participating in the dynamics are well described at the EOM-CCSD level of theory.

| Strong-field and ionization dynamics
The TDCC theory with static orbitals has an inherent instability.If an intense, resonant laser pulse is applied to the TDCC ground state of a system dominated by the Hartree-Fock ground-state determinant, the ground state is rapidly depleted and the state of the system thus becomes essentially orthogonal to the Hartree-Fock determinant.Yet, by construction, the intermediate normalization condition exp i¼ 1 must hold at any time t.This causes violent behavior of the amplitudes.Indeed, Pedersen and Kvaal 60 found that even with untruncated cluster operators, TDCC theory fails when the resonant laser pulse is strong enough.This implies, for example, that Rabi oscillations between the ground state and an excited state are practically impossible to describe within static-reference TDCC theory.
Choosing the time-dependent Brueckner determinant-the single determinant with the greatest overlap with the TDFCI wave function at any time t-as reference determinant, one would expect the intermediate normalization condition to place much less severe demands on the amplitudes, at least when the Brueckner weight is sufficiently close to 1. Kristiansen et al. 93 showed that the time-dependent reference determinant in TDNOCC theory (equal to OATDCC theory without splitting of the orbital space) is, in fact, an excellent approximation to the Brueckner determinant and that, therefore, TDNOCC theory shows improved numerical stability compared with staticreference TDCC theory.An example is given in Figure 2, which shows the weight of the reference determinant in TDNOCCD and TDCCSD simulations of the Be atom exposed to an intense near-resonant laser pulse, comparing with the weight of the time-dependent Brueckner determinant and of the Hartree-Fock ground-state determinant in TDFCI simulations.Also shown are the norms of the doubles amplitudes.The weight of the reference determinant is defined by, 93 where is the two-component state vector representing the reference determinant, either time-dependent for TDNOCC, OATDCC, and TDOCC theory or the static Hartree-Fock determinant for conventional TDCC theory.The extreme behavior of the amplitudes in TDCCSD theory is clearly correlated with low Hartree-Fock weights and enhanced stability is obtained in TDNOCCD theory (labeled OATDCCD theory in Figure 2) where the reference determinant is an excellent approximation to the Brueckner determinant.Instabilities may still occur in TDNOCCD theory, 93 however, and an absolutely stable TDCC theory likely requires a multireference Ansatz with time-dependent orbitals.For a physically correct description of electronic dynamics with nonvanishing ionization probability, including high-harmonic generation (HHG), the underlying basis set used to expand the time-dependent spin orbitals must support the electronic continuum.In the TDOCC approach of Sato et al., 85 a FEDVR is used along with absorbing F I G U R E 1 Controlling the ratio of CCSD energy level populations for LiH by altering the chirp rate of a laser pulse.The squares mark reference populations from TDFCI simulations.The aug-cc-pVDZ basis set was used.(Reprinted from Ref. [61].Copyright ©2020.T. B. Pedersen, H. E. Kristiansen, T. Bodenstein, S. Kvaal, and Ø. S. Schøyen.Published by American Chemical Society.)boundary conditions and splitting of the spin-orbital space into active and inactive subspaces.The main goal of Sato et al. 85 was to simulate HHG processes induced by a few-cycles, near-infrared (800 nm) laser pulse for atoms, including estimation of one-and two-electron ionization probabilities.To keep the computational cost reasonably low, the core electrons were frozen.
The HHG spectra were computed from the absolute square of the Fourier transform of the induced dipole acceleration, while Sato et al. estimated one-and two-electron ionization probabilities as the probabilities of finding one or two electrons outside a sphere of radius 20 a:u: around the nucleus.Examples are given in Figures 3 and 4 for ionization probabilities and HHG spectra, respectively, where results from TDOCCD and TDOCCDT simulations are compared with time-dependent complete active space self-consistent field (TDCASSCF) results which, with the same orbital-space splitting, can be regarded as TDFCI benchmark results.Both HHG spectra and ionization probabilities show considerable electron correlation effects and triple excitations must be included in the TDOCC treatment to get results close to the TDFCI limit.Note, however, that even without triple excitations, the TDOCCD model yields very significant improvements over the uncorrelated Hartree-Fock method.
Sato and coworkers have since 2018 published several studies of ionization and HHG processes with less computationally demanding CC-like approximations based on the FEDVR.The simplest time-dependent orbital-optimized coupled-electron pair (TDOCEPA0) 88 approximation, which can be viewed as a linearization of the TDOCCD method, was shown to yield results roughly on par with the parent TDOCCD theory at low and intermediate laser intensities, while the effects of electron correlation are somewhat overestimated at higher intensities.Although the TDOCEPA0 approximation carries the same formal computational complexity, O N 6 À Á , as the parent TDOCCD method, the linearization leads to significant computational simplifications, including halving the number of equations to be propagated due to the symmetry Further reductions in computational time was obtained with the time-dependent orbital-optimized second-order Møller-Plesset (TDOMP2) 88 model, where only terms through second order in the fluctuation potential are retained in the TDOCCD Hamilton function ℋ, leading to an approximation scaling as O N 5 ð Þ.The TDOMP2 model is related to the time-dependent second-order CC model, TDCC2, 94 replacing single-excitation amplitudes by orbital rotations.The TDOMP2 method tends to overestimate electron-correlation effects at a wide range of intensities, but strikes a reasonable balance between accuracy and efficiency.Pathak, Sato, and Ishikawa 95 also computed TDOMP2 results with those obtained with a variant of TDCC2 theory where the single-excitation amplitudes are included alongside orbital optimization.It was found that the TDOMP2 method yields superior HHG spectra compared with the TDCC2-like method proposed by Pathak, Sato, and Ishikawa. 95Compared with higher-level theories, it was concluded that the accuracy of the TDOMP2 model is moderate.
Recently, Pathak, Sato, and Ishikawa 90 also introduced the TDOCCDT(4) method, which includes tripleexcitation corrections through fourth order in the Hamilton function and shows a computational complexity of O N 7 ð Þ, intermediate between the TDOCCD (O N 6 À Á ) and full TDOCCDT (O N 8 À Á ) models.Test calculations on the Ne and Ar atoms indicate that the TDOCCDT(4) method gives results almost indistinguishable from the full TDOCCDT model and the TDCASSCF model for HHG spectra and one-and two-electron ionization probabilities.
F I G U R E 2 TDCCSD, OATDCCD, and TDFCI simulations of Be with the cc-pVDZ basis exposed to a laser pulse with peak electric-field strength 1 a:u: and carrier frequency ω ¼ 0:2068175 a:u: (Reprinted from Ref. [93] with the permission of AIP Publishing.)

| Linear and low-order nonlinear optical properties
Transient absorption spectroscopy 10 is an important application of time-dependent electronic-structure theory, yielding time-resolved spectra containing much richer information than conventional steady-state spectroscopy.Skeidsvoll, Balbi, and Koch 96 adapted the theory of transient absorption by Wu et al. 10 to bivariational TDCC theory and presented simulations of transient core-level spectra of the LiH and LiF molecules at the TDCCSD level of The probabilities, as a function of time, of finding one (a) and two (b) electrons outside a sphere of radius R 0 = 20 a.u.Comparison of the results of TDHF, TD-OCCD, TD-OCCDT, and TD-CASSCF methods.(Reprinted from Ref. [85] with the permission of AIP Publishing).

F I G U R E 4
The HHG spectra of Ar exposed to a laser pulse with a wavelength of 800 nm and an intensity of 6 Â 10 14 W=cm 2 .Comparison of the results of TDHF, TD-OCCD, TD-OCCDT, and TD-CASSCF methods.The inset shows a close-up of the spectra from 50th to 80th harmonic order.(Reprinted from Ref. [85] with the permission of AIP Publishing.) theory.Using a resonant valence-exciting pump pulse followed by a core-exciting probe pulse, Skeidsvoll, Balbi, and Koch observed oscillations of intensities with the pump-probe delay caused by interference in the wavepacket generated by the pump pulse.Using a static reference determinant, ionization dynamics is out of reach and, therefore, only weak lasers were studied.In fact, the wavepacket generated by the pump laser is overwhelmingly dominated by the electronic ground state: the LiF ground-state population was found to be roughly 99:5% in Ref. [61].Consequently, the pump-probe spectra were dominated by features ascribable to ground-to-excited state transitions, that is, essentially linear absorption spectra.
Conventional linear absorption spectra can be extracted from simulations, too.A molecule, initially in its ground state, is exposed to a weak electric-field kick-a weak delta-function shaped laser pulse-which induces transitions from the ground state to all excited states that can be populated within the electric-dipole selection rules.For sufficiently weak fields, multiphoton transitions are virtually absent.Moreover, the excited-state populations are so small that transitions between them are virtually absent, too, yielding a linear absorption spectrum.The absorption cross section may then be obtained from the Fourier transform of the induced electric-dipole moment.The great advantage of this approach to linear absorption spectra is that the entire spectrum, including the high-frequency core-valence transitions, is obtained from 1 to 3 simulations, one for each polarization direction to emulate random orientation of the sample relative to the propagation direction of the laser.The main challenge is the very long simulation times required to achieve sufficient resolution of the simulated spectrum.The simulation time can be reduced by about a factor 5 using Padé approximants for the Fourier transformation. 97ristiansen et al. 98 used linear absorption spectra, including core excitations, generated from the dipole moment induced by an electric-field kick to validate their implementation of the TDCC2 model 94 by comparing with results from linear response theory.They also presented a derivation of the equations of motion for the TDOMP2 and TDNOMP2 (where the spin orbitals are required to be biorthonormal instead of orthonormal) models based on exponentially parameterized orbital rotations and the Euler-Lagrange Equation ( 20), and found that, despite the full orbital relaxation included in the TDOMP2 model, no significant improvement over TDCC2 spectra was obtained in the core region of the spectrum.
However, for frequency-dependent polarizabilities and hyperpolarizabilities extracted from TDCC2, TDCCSD, and TDOMP2 simulations with ramped monochromatic continuous-wave lasers as suggested by Ding et al., 14 Kristiansen et al. 98 found that TDOMP2 theory outperforms the TDCC2 model, producing linear and nonlinear response functions much closer to the full TDCCSD results.While the linear absorption spectra were in perfect agreement with results from linear response theory (to within the resolution of the Fourier transformation), the TDCC2 and TDCCSD polarizabilities and hyperpolarizabilities were found to deviate slightly.These deviations are most likely caused by the nonperturbative nature of TDCC simulations and by nonadiabatic effects not being entirely removed by the single-cycle ramping.

| ELECTRONIC DYNAMICS WITH EOM-CC THEORY
With a finite basis, the eigenfunctions of the time-independent Hamiltonian b H 0 of the many-electron systemparticlesystem-the stationary states-can in principle be used to expand the time-dependent wave function because the continuum (for which normalizable eigenfunctions do not exist) becomes discretized.Starting from the EOM-CC Ansätze 30,66-70 for the left and right eigenfunctions of the similarity-transformed Hamiltonian, one may formulate an alternative theory, based on the eigenstate superposition approach, which converges to the correct FCI limit but provides different results than TDCC theory with truncated cluster operators.This idea was used by Kjønstad and Koch 99 to recast the Born-Huang approximation 4 to the fully coupled time-dependent electronic-nuclear wavefunction of a molecular system in the framework of CC theory.

| Equation-of-motion coupled-cluster theory
In EOM-CC theory, 30,66-70 the time-independent excited states are parameterized on top of the CC ground state as where b T is the time-independent cluster operator of the ground-state CC wavefunction and j Φ 0 i is a static reference determinant, typically the Hartree-Fock ground-state determinant.The linear excitation and de-excitation operators b R I and b L I , respectively, are defined by The coefficients of b L I and b R I are collected in the left and right eigenvectors of the field-free similarity transformed Hamiltonian b where E ¼ diag E I ð Þ contains the energies of the ground and excited states, and the columns of R and the rows of L define the right and left eigenvectors, respectively, such that the biorthonormality condition is fulfilled (L I is the Ith row of L and R J the Jth column of R).

| Time-dependent equation-of-motion coupled-cluster theory
The most direct formulation of time-dependent EOM-CC (TD-EOM-CC) theory is to expand the time-dependent bra and ket wavefunctions in the basis of field-free EOM-CC states Equations of motion for the time-dependent expansion coefficients e C I t ð Þ and C I t ð Þ can either be obtained from the time-dependent bivariational principle or by requiring that the time-dependent Schrödinger equation holds for the left and right wavefunctions, that is, Insertion of the Ansätze (44) and ( 45) into the TDSE and projecting onto the EOM-CC ket (bra) state yields where H IJ t ð Þ are the matrix elements of the Hamiltonian b H t ð Þ in the EOM-CC basis, When the time-dependent expansion coefficients have been determined, the expectation value of an arbitrary operator b Ω is given by in accordance with the Hellmann-Feynman theorem.The matrix elements Ω IJ are defined as in Equation ( 50).Note that since h e Ψ t ð Þ j and j Ψ t ð Þi are not Hermitian conjugates, the expectation value of b Ω is generally complex valued.Of course, the indefinite inner product, Equation ( 18), could also be used here to obtain real values for hermitian operators.
For a time-dependent Hamiltonian on the form b , the matrix elements of the Hamiltonian can be written as where the matrix elements V IJ t ð Þ constitute a non-hermitian matrix with elements defined as in Equation ( 50).This approach has been used by different groups to study explicitly time-dependent optical processes within the EOM-CC framework with b , where b μ is the electric-dipole operator and F t ð Þ is the spatially uniform electric field of the laser.
For example, Sonk, Caricato, and Schlegel 100 studied the optical response of butadiene to short, intense laser pulses while Luppi and Head-Gordon 101 used the method to compute HHG spectra of H 2 and N 2 .In both articles, the transition dipole matrix D IJ ¼ e Ψ I jb μjΨ J D E was symmetrized according to because the authors were worried that the non-hermitian tranisiton dipole matrix in Equation ( 52) could lead to dynamics that do not conserve the norm of the wavefunction.However, if the left and right TDSE are satisfied, it follows that Hence, the biorthonormality of the left and right wavefunctions is conserved, provided that e Note that with symmetrization (and real energies), the Hamiltonian matrix (52) becomes hermitian and the left and right expansion coefficients become complex conjugates, e With the symmetrization (53), which was also used in Ref. [61] for test purposes (albeit with transition dipole moments from CC response theory rather than EOM-CC theory), the left and right expansion coefficients are related by complex conjugation.In a recent publication, 102 Skeidsvoll et al. used the TD-EOM-CCSD method without symmetrization to simulate weak-field attosecond processes in small molecules-that is, distinct left and right expansion coefficients were retained.Core-level pump-probe spectra of LiH and LiF were compared with TDCCSD results from Ref. [96] and found to be in good agreement.
The TD-EOM-CCSD approach can be much more efficient than the TDCCSD method if the participating states can be limited by a careful selection procedure, 102 but the main computational drawback of the approach is the full diagonalization of the Hamiltonian required for completely general quantum dynamics simulations.It is generally hard to predict a priori how many and which states are required to describe a given dynamical process and a prohibitively large number of states may be needed.Skeidsvoll et al. 102 used an asymmetric band Lanczos algorithm combined with state selection criteria based on transition strengths and frequencies to successfully limit the number of participating states, including high-lying core-excited states using core-valence separation. 103,104nother alternative is to avoid the diagonalization problem entirely by propagating directly in a determinant basis.While this effectively removes all issues related to diagonalization of large matrices-such as stability, selection of states, and large memory usage-the computational complexity is moved to the time propagation, exactly as in TDCC theory albeit with the benefit of linear parameterization.This would be especially disadvantageous in simulating transient absorption spectroscopy, since the ability to analytically propagate the wavefunction after the pulses are switched off is removed, exactly as in TDCC theory.Still, this may be an attractive approach in cases where one is interested in broad frequency ranges and large number of states.

| Time-dependent formulation of linear absorption spectra
6][107][108][109][110][111] The starting point is the line-shape function obtained from Fermi's Golden Rule, 112 which results from solving the TDSE to first order in time-dependent perturbation theory with the ground-state wavefunction as initial condition and using the rotating wave approximation (i.e., only ground-to-excited state one-photon transitions are included).The line-shape function thus implicitly assumes weak-field perturbations.Adapting to EOM-CC states and assuming the electric-dipole approximation for the semiclassical matter-field interaction, Fermi's Golden Rule for the line-shape function can be written as 107 where α x,y,z f g, ω I ¼ E I À E 0 , and we have assumed a Lorentzian line shape with fixed lifetime 1=γ for all excited states: Here, we follow Nascimento and DePrince 107 and include the ground state in the summation over states I.The ground state was not included in the work of Park, Perera, and Bartlett. 109Note that the line-shape function may be complex in EOM-CC theory; had we used the indefinite inner product instead, only the real part of this expression would be used.
Using that ℒ ω; ω I , γ ð Þis the Fourier transform of exp iω I t À γjtj ð Þ , the line-shape function can be rewritten as where the right and left dipole functions are defined as The line-shape function (57) thus may be evaluated by Fourier transformation of either of the dipole ACFs after multiplication by the damping factor exp Àγjtj ð Þ.Here, Only one of these propagations need to be performed to compute the dipole ACF.While DePrince and coworkers [105][106][107]111 used a projection procedure for the TDSE for either j M α Àt ð Þi or h e M α t ð Þ j, Park, Perera, and Bartlett 109 used the evolution of the dipole operators in the Heisenberg picture. Forally, at least, the two approaches are equivalent.The absorption spectrum is then obtained from Since this approach does not require explicit diagonalization of the similarity transformed Hamiltonian, it provides a direct route to the study of core excitation spectra 106,109 and has been generalized to include scalar relativistic effects at the exact two-component (X2C) level with a fifth-order Douglas-Kroll-Hess Hamiltonian. 109Scalar relativistic effects including spin-orbit coupling at the X2C level was explicitly addressed in Ref. [108], where Koulias et al. illustrated both frequency shift, activation of spin-forbidden transitions, and energy splitting of the 2 P 1=2 and 2 P 3=2 states in atoms and cations of the alkali and alkaline earth metal groups.
It should be noted that the methodology outlined above is not restricted to the electric-dipole approximation and can be generalized to account for beyond-dipole effects. 107,110Examples are shown in Figures 5 and 6, where the method was applied to compute linear absorption and isotropic electronic circular dichroism spectra in substituted oxiranes.Park, Perera, and Bartlett 110 reported an implementation which included higher-order multipole functions corresponding to full second-order oscillator strengths, simulating the 3p !4d quadrupole-allowed transition in the pre K-edge region of Ti 4+ and TiCl 4 .This study also includes scalar relativistic effects.
Since the computational cost of these methods is dominated by the time integration and since long simulation times are required to achieve sufficient resolution in the Fourier transform of the ACF, some effort has been spent on algorithms designed to accelerate the simulations.Nascimento and DePrince 106 used Padé approximants of the Fourier F I G U R E 5 Linear absorption spectra for substituted oxiranes computed at the TD-EOM-CC2/aug-cc-pVDZ level of theory.The solid black lines correspond to artificially broadened stick spectra obtained from standard, frequency-domain EOM-CC2 computations.The labels NH 2 , F, OH, and CH 3 represent the oxirane substituents.(Reprinted from Ref. [107] with the permission of AIP Publishing.) transform, achieving an order of magnitude speedup with small or no errors in the linear absorption spectrum.Later, however, Cooper et al. 111 reported that the Padé approximants may give faulty results in dense spectral regions.The same group later reported an implementation using a short iterative Lanczos (SIL) integration scheme. 111This method utilized the fact that the main computational cost in the time propagation is the evaluation of the right hand side of the TDSE, which amount to a matrix vector product HC.In the SIL approach, a tridiagonal approximation H k to the Hamiltonian matrix is constructed and used to propagate a moment vector within a Krylov subspace of dimension k.Due to the simple Hamiltonian form, they used the matrix representation of exp iH k dt À Á directly in the time propagation.The approximate Hamiltonian and the corresponding Krylov subspace must be regenerated regularly, but still the authors report up to an order of magnitude speedup relative to the RK4 integrator.Although they used an algorithm designed for hermitian matrices, it was shown that the SIL method generates frequency spectra with mostly negligible differences from those generated with the RK4 integrator.

| TIME-DEPENDENT VIBRATIONAL COUPLED-CLUSTER THEORY
Vibrational coupled-cluster (VCC) theory refers to the application of CC theory to the nuclear Schrödinger equation in the adiabatic Born-Oppenheimer approximation.There are two distinct flavors of VCC theory: First, a basis-free method based on bosonic CC theory 54,[113][114][115] was developed into sophisticated VCC theory by Banik, Pal, and Prasad 116 and Faucheaux and Hirata, 117 who coined the acronym XVCC.The second approach is the modal approach to VCC F I G U R E 6 Electronic circular dichroism spectra for the substituted oxiranes computed at the EOM-CC2/aug-cc-pVDZ level of theory in the length (left panels) and velocity (right panels) gauges.Colored lines correspond to TD-EOM-CC2-derived data, while the solid black lines are frequency-domain spectra obtained by artificially broadening the corresponding stick spectra.(Reprinted from Ref. [107] with the permission of AIP Publishing.) theory, simply termed VCC, developed by Christiansen. 118,119RT propagation has been developed in both flavors.Indeed, in the XVCC case, this is where it started.

| Bosonic VCC theory: XVCC
The bosonic VCC theory starts with the harmonic approximation, writing the nuclear Hamiltonian as an Mdimensional harmonic oscillator plus anharmonic perturbations.Each of the M modes have associated harmonic-oscillator ladder operators b a with joint vacuum j 0i, that is, the ground state of the M-dimensional harmonic oscillator.The XVCC nuclear wavefunction is given by with the untruncated cluster operator b S defined by The shown terms up to second order define the SUB2 approximation, while higher-order approximations, denoted SUBN, truncate b S at the Nth order.The theory now proceeds as in traditional CC theory.It is to be remarked, that the SUB2 approximation generates a wavefunction which is a squeezed state, that is, a general complex-valued Gaussian.
In Ref. [120], Prasad introduced TDCC theory for the study of Franck-Condon spectra, that is, the nuclear transition probabilities to the various vibrational eigenstates upon (instantaneous) excitation from the nuclear ground state at the equilibrium geometry to an excited electronic surface.This seems to be the first application of CC theory to the vibrational Schrödinger equation, and at the same time it is an early application of TDCC theory.Notable here is that only the ket ðj ΨiÞ is propagated in time, since the transition probability is obtained in terms of an ACF only dependent on this ket.Prasad presented an application to a two-dimensional (2D) model system previously studied by Heller using frozen Gaussians, 121,122 pointing out similarities and differences.For example, the TDCC treatment in the SUB2 approximation is roughly equivalent to the thawed Gaussian approximation.In a follow-up study, Sastry and Prasad 123 applied the above methodology to the Beswick-Jortner model 124 of photodissociation of the form ABC !AB þ C.
In Ref. [125], Latha and Prasad studied the possibility of using TDCC theory to describe nonadiabatic dynamics on conically intersecting potential-energy surfaces.In addition to treating multiple electronic surfaces, the basis is made time dependent using a time-dependent self-consistent field (TDSCF) procedure.Although a very simple model system was employed, the findings were encouraging.In particular, spurious peaks in the ACF arising from TDSCF theory alone were strongly alleviated by the CC couplings.

| Modal VCC theory
In the VCC theory developed by Christiansen, 118,119 the exponential Ansatz is applied to the vibrational Schrödinger equation for nuclear motion.Formally, the theory can be described as standard single-reference CC theory with multiple species of distinguishable particles, called modes.The starting point is the vibrational self-consistent field (VSCF) procedure, 119 which approximates the nuclear wavefunction as a Hartree product of M functions ϕ m 0 q m ð Þ called modals, where m ¼ 1, ÁÁÁ, M. Similar to the Hartree-Fock procedure, VSCF also produces excited Hartree products, where s ¼ s 1 , ÁÁÁ, s M ð Þenumerates the modals, and where q ¼ q 1 , ÁÁÁ, q m ð Þare the nuclear coordinates.A certain arbitrariness exists in the choice of the mode coordinates q m .The VSCF ground state Hartree product Φ 0 q ð Þ and the excited products form the many-mode Hilbert space basis.A second-quantization formalism can be set up which in a natural manner defines the exponential ansatz for the vibrational Schrödinger equation.Each modal ϕ m s m is associated with a creation operator b a †,m A general cluster operator reads b T ¼ P μ τ μ b X μ , with μ being a generic index enumerating the various m-mode simultaneous excitations, up to M simultaneous excitations.For example, a single-mode cluster operator reads where A m is the number of modals for mode m, while a two-mode cluster operator reads The cluster operator is usually truncated at a number n of simultaneous mode excitations, denoted the VCC[n] approximation. 118part from the qualitative differences arising from having multiple distinguishable particle species, the Hamiltonian in VCC theory is also qualitatively different, since the particles are now modes, and there are, in principle, arbitrarily many such particles interacting at the same time.In particular, the Baker-Campbell-Hausdorff expansion does not truncate for VCC for the exact Born-Oppenheimer potential-energy surface.The Hamiltonian is typically truncated at a maximum number of interacting modals (e.g., using a sum-of-products approximation).
The Christiansen group has also developed RT time-dependent VCC (TDVCC) theory using the bivariational framework in Section 3. In Ref. [126], Hansen et al. introduced the TDVCC formalism, along with an analysis of the separability of the bra and ket wavefunctions as well as the corresponding extensivity of expectation values.The authors also considered imaginary-time propagation for locating the ground-state solution.These theoretical results are also highly relevant for electronic-structure TDCC theory.The authors present an implementation of the TDVCC [2] method, using the Dormand-Prince 8(5,3) explicit Runge-Kutta method with adaptive step size control, see Ref. [127], Section II. 5.
In order to verify their implementation and analysis, the authors tested their TDVCC [2] code on the 2D Hénon-Heiles potential, as well as calculations on the water and formaldehyde molecules, using approximate potentials on a sum-over-products format coupling at most two modes per term, generated using Gaussian process regression. 128One interesting finding is that round-off errors in the asymptotic region of imaginary-time propagation can, even for such an accurate Runge-Kutta integrator, lead to slightly incorrect exponential decay rates, thereby predicting slightly wrong excitation energies.The authors also studied the integration of a system driven by an explicitly time-dependent laser pulse.
In a follow-up study, 129 an automated implementation of the full hierarchy of TCVCC[n] approximations was presented.Again, the Dormand-Prince 8(5,3) method was the chosen integrator for the numerical studies.The authors note that even though the integrator is not symplectic, it is sufficiently accurate so that there are no stability or nonconservation problems often associated with explicit Runge-Kutta methods. 91he authors performed several detailed numerical experiments, demonstrating the convergence to the vibrational FCI (VFCI) limit for a 5-mode system (formaldehyde).An extended truncation scheme inspired by single-reference based multireference theory 130 is also introduced, called VCC[kextn], where a single mode is chosen to have n À k more excitations than the VCC[k] truncation scheme.The convergence toward the VFCI limit was demonstrated to be enhanced when a single mode dominated the dynamics.
The authors also discussed ACFs in detail with respect to their separability properties for separated noninteracting systems, singling out the ACF defined by A t i , valid in standard hermitian dynamics for a real Ψ 0 ð Þ and a real Hamiltonian.However, the authors found severely unphysical behavior of this ACF.Finally, we mention that the authors applied their implementation to a larger system, studying the intramolecular vibrationalenergy redistribution of the imidazole molecule (with 21 modes) using an accurate many-term potential-energy surface.
In a related study, Hansen, Madsen, and Christiansen 131 implemented the full time-dependent extended CC method of Arponen, 54 that is, the TDEVCC method.Although unfeasible for larger systems than, say, M ¼ 6 modes, the TDEVCC method utilizes a double exponential Ansatz for the bra and ket vectors, implying full multiplicative separability and corresponding separability of expectation values.The authors observed that for ground-state energy calculations, EVCC theory does not offer a significant improvement over "plain" VCC theory, especially taking the computational cost into consideration.This is a finding consistent with the conventional wisdom in electronic-structure theory. 132However, the authors noted that for time-dependent calculations, TDEVCC[k] performs in general much better than TDVCC[k] with regards to the closeness of both the bra and the ket to the TDVFCI limit and accuracy of expectation values.Both ACFs of type A and B are correctly separable with TDEVCC[k].
The Christiansen group has also developed and implemented a VCC analogue of orbital-adaptive TDCC 49 (cf., Section 3), called time-dependent modal VCC (TDMVCC) theory.In Ref. [133], Madsen et al. presented an advanced implementation of orbital-adaptive theory with very promising results, including the apparent cure of exploding amplitude norms in TDVCC[n] calculations for the water potential-energy surface.This corroborates findings by Kristiansen et al., 93 who applied TDNOCC theory to electronic systems.

| Finite-temperature theory
Development of viable computational tools for the study of many-body quantum systems, especially in the condensed phase, at finite temperature is an active research area, see, for example, Refs.[134-137] for recent work within density-functional theory.The earliest efforts to cast quantum mechanics in a thermodynamical framework were based on the close resemblance of the statistical partition function and the quantum-mechanical evolution operator with imaginary time.The Matsubara formalism 138 for Green's functions in quantum-field theory is a notable example of this connection.However, when the time has been rotated to the imaginary-time axis, the dynamics are lost.To circumvent this problem there are several techniques to re-introduce real time, for example, thermofield dynamics and the Keldysh formalism. 139,140he first to explore CC theory in a finite-temperature setting were Altenbokum and coworkers. 115,141They used a density-matrix formulation with the Bloch equation replacing the Schrödinger equation.This, however, requires the knowledge of the full spectrum of the Hamiltonian thus quickly making the solution prohibitive.Several years later the thermal cluster cumulant (TCC) method was developed by Sanyal, Mandal, and Mukherjee as an extension of the thermofield dynamics with the CC method, requiring a thermal formulation of Wick's theorem. 142,143This resulted in a method resembling the cumulant expansion from statistical mechanics, and bypassed the need to know the full spectrum of the Hamiltonian.][144][145][146] In 2018, White and Chan 147 started from an explicitly time-dependent formulation of CC theory and replaced time by an imaginary time.This replaced the time-dependent formulation with a temperature-dependent formalism, dubbed the finite-temperature coupled-cluster (FTCC) method, and was further studied in 2020. 148Their formulation of FTCC theory is slightly different from the TCC model by Sanyal, Mandal, and Mukherjee, 142 but the methods are equivalent and lead to the same set of equations.
At the same time, Hummel 149 developed an imaginary-time time-dependent truncated CC method for the application to systems at finite temperature.The truncation scheme is the direct-ring coupled-cluster doubles (drCCD) method (i.e., the direct random-phase approximation) and leads to much more cost-effective equations than the FTCC formalism at the price of reduced accuracy.
Following shortly after these publications, Harsha, Henderson, and Scuseria 150 developed an imaginary-time TDCC method for finite-temperature systems from the thermofield dynamics formalism.Results from this method were published already in an earlier work comparing with a thermal configuration-interaction method. 151y including dynamics the study of non-equilibirium systems becomes possible.White and Chan 152 utilized the Keldysh formalism to extend the imaginary-time formalism to include a real component for RT dynamics. 140This method was dubbed the Keldysh coupled-cluster (Keldysh-CC) method.Up to this point all time-dependence and temperature-dependence was kept in the cluster amplitudes.However, in 2021 Peng et al. 153 extended the Keldysh-CC method to include orbital rotations.They formulated the Keldysh-OCC method as an extension of the TDOCC method 85 to finite-temperature systems.

| Sub-system embedding
Kowalski and Bauman 154 have developed the sub-system embedding sub-algebra CC formalism, which is a generalization of the complete active space CC (CASCC) formalism of Piecuch, Adamowicz, and coworkers. 130,155The excitations that stay completely in the CAS is an example of a sub-system embedding sub-algebra.In Ref. [154], the formalism is extended to the TDSE.Moreover, a CASCC generalization of unitary CC (UCC) is studied, including variational formulations of the dynamical equations of motion.The UCC method is one of the main contenders for quantum advantage on noisy intermediate-scale quantum (NISQ) devices. 156

| Green's function methods
In coupled-cluster Green's function (CCGF) theory, 157 the goal is to compute a quantum mechanical Green's function by means of CC theory.For example, the one-body retarded Green's function defined by is a causal propagator that expresses the amplitude of electron/hole propagation from time t 1 to a later time t 2 .Here, b a p t ð Þ ¼ e i b Ht b a p e Ài b Ht is the Heisenberg representation of b a p and θ is the Heaviside step function.Intuitively, G À contains information about causal single-particle processes in a quantum system starting out in the ground state, and can be used to compute a host of properties. 158The Green's function is closely related to the ACFs described elsewhere in this review.
In CCGF theory, the bivariational approximations for j Ψ 0 i and hΨ 0 j constitute the starting point, see Section 3. Conventionally, a momentum-space representation of G À is sought and approximated using many-body perturbaton theory.In recent work, however, RT propagation methods for j Ψ 0 i in conjunction with a cumulant approximation 159 (i.e., exponential Ansatz) was explored.
In Refs.[160, 161], Rehr, Vila and coworkers used RT propagation within EOM-CC theory, obtaining a nonperturbative expression for the cumulant appearing in G À in terms of the solution to a set of coupled first-order, nonlinear differential equations.The primary aim was to study x-ray absorption spectra of molecular systems, and it was shown that the nonlinear terms of the cumulant expansion yields significant improvements over the traditional linear approximation.In Ref. [162], the approach is further extended and applied to the core-hole spectral function for small molecular systems.

| CONCLUDING REMARKS
While still in its infancy, TDCC theory is emerging as a versatile tool in computational molecular science.As one might perhaps expect from its unquestionable status as the high-accuracy method in time-independent quantum chemistry, it can provide high accuracy relative to TDFCI theory for electronic and vibrational quantum dynamics.With timedependent orbitals and sufficient flexibility in the basis set used to expand these orbitals, TDCC theory has the potential to make decisive contributions to the understanding of molecular processes on the timescale of the electron-the attosecond timescale-where not only bound states, but also the electronic continuum must be taken into account.With relatively little development effort (since no response equations need to be implemented), TDCC theory can be used to compute full linear absorption spectra, including core-level excitations.Linear and low-order nonlinear response functions may be extracted from relatively short simulations.With proper use of Arponen's bivariational formulation, essentially all information that can be extracted from hermitian quantum dynamics can also be extracted from TDCC simulations despite the non-hermitian formulation.
Several challenges need to be overcome, however.The computational cost of TDCC theories is very high compared with the most widely used method, TDDFT, and two issues need to be addressed.First, the computational complexity or at least the prefactor of the evaluation of the function f in Equation ( 30) must be reduced.A general algorithm for this is significantly more challenging to formulate than reduced-scaling algorithms for the ground state, since different external driving forces may produce dramatically different responses in the wavefunction.For weak driving fields, for example, it is (very) small changes in the cluster amplitudes that produce the oscillations of interest, making screening procedures difficult to implement with full controllability of the accuracy, whereas strong fields can produce wavefunctions with very wide spatial distribution.
Second, the number of f-evaluations must be kept at a minimum, requiring careful selection of the numerical integrator and improved signal processing.Another major challenge, which is shared by all electronic quantum dynamics methods, is the representation of the electronic continuum.Ideally, one would wish to have the accuracy of mesh-based approaches like FEDVR at the cost of Gaussian-based electronic-structure theory.
Third, nuclear motion needs to be included to reliably extend simulation times beyond a few femtoseconds.This can be approached with either classical or quantum nuclear motion but will eventually require a time-dependent multireference CC wavefunction to describe, for example, a photo-induced chemical reaction.

RELATED WIREs ARTICLES
Real-time time-dependent electronic structure theory It is interesting to compare this with the ACF used by Pedersen and Kvaal, 60 Equation (31), which is the average of A and B. In the work of Pedersen and Kvaal, no unphysical values were observed for B t 0 , t ð Þ unless the integration of the equations of motion failed due to ground-state depletion (in which case also A t 0 , t ð Þ becomes ill-behaved).The third ACF studied was based on the relation Ψ 0