# Three Lagrangians for the complete-active space coupled-cluster method
**Authors**: Simen Kvaal
> Hylleraas Centre for Quantum Molecular Sciences, Department of Chemistry, University of Oslo, P.O. Box 1033 Blindern, N-0315 Oslo, Norway
(July 13, 2023)
## Abstract
Three fully variational formulations of the complete-active space coupled-cluster (CASCC) method are derived. The formulations include the ability to approximate the model vectors by smooth manifolds, thereby opening up the possibility for overcoming the exponential wall of scaling for model spaces of CAS type. In particular, model vectors of matrix-product states are considered, and it is argued that the present variational formulation allows not only favorably-scaling multireference coupled-cluster calculations, but also systematic correction of tailored coupled-cluster calculation and of quantum chemical density-matrix renormalization group methods, which are fast and polynomial scaling, but lacks the ability to properly resolve dynamical correlation at chemical accuracy. The extension of the variational formulations to the time-domain is also discussed, with derivations of abstract evolution equations.
## I Introduction
In this article, we discuss some Lagrangian formulations of the complete-active space coupled-cluster (CASCC) method introduced in the 1990s by Oliphant, Adamowicz and Piecuch [1, 2, 3], see also the work by Stolarczyk [4]. For a comprehensive review with literature, see Ref. [5]. More specifically, we will consider the semilinear wavefunction formulation of Adamowicz and coworkers [6, 7] but linked energy-independent working equations to the best of our knowledge first proposed rather recently by Kowalski [8].
The motivation for the present work is threefold: Firstly, any multireference approach employing a complete model space (i.e., full CASCI model vectors) face an exponential wall of scaling for large systems. Systems that are truly strongly correlated requires CASs with the number of orbitals proportional to the number of particles, which is clearly not feasible. A formulation of CASCC that allows approximations of the CASCI model vector in terms of differentiable manifolds such as matrix-product states (MPS) and more generally using the density-matrix renormalization group algorithm (DMRG) [9, 10, 11, 12, 13] could enable affordable and accurate multireference calculations for quite large systems. The MPS ansatz is excellent at resolving static correlation at polynomial cost, and there is a growing interest in its application in quantum chemistry. Key to achieving a successful combination of coupled-cluster methods with such approximations is a fully variational framework. Secondly, literature searches for fully variational formulations of the semilinear CASCC method similar to the formulation which has now become standard in traditional single-reference coupled-cluster (SRCC) theory have been fruitless. The Lagrangian formulation of SRCC was the culmination of several independent researchers’ efforts in the 1970–1980s to provide easily computable molecular properties and response to external perturbations [14, 15, 16, 17, 18, 19, 20].
Although the standard Lagrangian can be used for CASCC when the linear CAS amplitudes are mapped to an exponential ansatz, see e.g. [5], a Lagrangian for the semilinear case is desired, since the mapping from linear to exponential parameterizations can be cumbersome and numerically ill-conditioned, especially when excited states are targeted. Moreover, the numerical task of solving a linear (non-Hermitian) eigenvalue problem is different from solving a polynomial root equation.
Thirdly, a fully variational formulation of CASCC would allow its use in real-time propagation of coupled-cluster quantum states, a topic which has garnered recent interest in the attophysics community [21, 22, 23, 24, 25, 26, 27]. As the MPS ansatz is also being employed for real-time propagation of quantum chemical systems [28, 29, 30], one can speculate that a time-dependent CASCC method with a MPS-based model space description can be far superior to existing approaches.
The CASCC method is well-known to suffer from some issues that originate with its single-reference formalism: The use of a formal reference in the CAS gives a bias towards this reference. This in turn leads to undesirable potential-surface discontinuities. On the other hand, CASCC is relatively straight-forward to implement and to analyze compared to true multireference CC approaches. Moreover, when absolute energies are compared, CASCC seems to have similar accuracy as the genuine MRCC approches [31]. A closely related method is the tailored coupled-cluster method (TCC) [32], an externally corrected CC approach that is actively being studied, including the DMRG-TCC method where the CAS model function is replaced by a matrix-product state and optimized with the DMRG algorithm [33, 34, 35, 36, 37]. The TCC method, however, suffers from a systematic error due to the CAS wavefunction being computed once and for all, and not being computed self-consistently with the external correction. The present work could therefore pave the way for updating schemes for the matrix-product state that reduce or altogether eliminate such systematic errors.
We will derive three fully variational formulations of CASCC, which we call Lagrangians even if some additional variables introduced are not strictly speaking Lagrange multipliers. After a brief discussion of the CASCC method in Section II, the first Lagrangian is introduced in Section III. The Lagrangian has a similar interpretation as the standard SRCC Lagrangian as a constrained optimization of the CC energy. This Lagrangian is applicable to standard CASCC calculations where the full CASCI model vector is used, and where the linked form of the amplitude equations is used. (These are equivalent to the unlinked amplitude equations in just about any practical truncation scheme.) Just as the SRCC Lagrangian features a Lagrange multiplier operator that supplements the exponential ansatz, the CASCC Lagrangian features, in a very natural manner, a “Lagrange multiplier” CASCI model vector, i.e., we have two model vectors in addition to two cluster operators. The working equations for the model vectors are a right-eigenvalue problem for the CAS projection of the connected Hamiltonian, a non-Hermitian operator, yielding the CASCI model vector as solution, and an inhomogenous left-eigenvalue problem, giving the dual CASCI model vector as solution. The asymmetry of the rôles of these vectors can be traced to the “ket-centric” method of development of CASCC as a projection technique.
In Section IV we introduce a differenetiable submanifold of the CAS in which to approximate the CASCI model vectors. (This manifold must not be confused with what some authors call “the projection manifold” in CC theory.) Thus, both model vectors are assumed to lie on the manifold, and the corresponding variational equations are derived. However, the asymmetry of the equations for the two model vectors implies that available codes for the variational approximation of non-Hermitian eigenvalue problems cannot be reused without modification. This motivates the second Lagrangian, introduced in Section V, where the symmetry is restored by force. The price to pay is additional Lagrange multipliers, which now live in the tangent spaces of the variational manifold at the approximate CASCI vectors. These are not necessary to solve for when the ground-state energy is the sole point of interest. In the full CASCI limit, the first and second Lagrangians are equivalent, but when restricted to the variational manifold, the Lagrangians effecively describe different approximations.
The third Lagrangian, introduced in Section VI, is equivalent to the second Lagrangian, but has the added benefit that it is on the form of an expectation value. It is therefore suitable for formulating an action integral, while it also gives additional insight into the second Lagrangan’s structure. The action-integral formulations of the first and third Lagrangians are discussed in Section VII, including their equations of motion.
Section VIII discusses some possible applications of the variational formulations, including iterative improvement of DMRG-tailored CC. Section IX sums up and presents some future perspectives.
## II The CASCC method
We discuss CASCC in sufficiently general terms as to be applicable to molecular systems, nuclear structure calculations, and even number-conserving boson systems with only small modifications. As usual in CC theory papers, we need to specify our notation and index conventions, which we attempt to keep at an economical level. We assume a standard setup with a computational Hilbert space $H$ generated by a fixed set of single-particle functions (SPFs), i.e., $H$ is the span of all possible linearly independent determinants (or permanents in the boson case) constructed with these functions. The computational Hilbert space is split into a model space $H_0$ and its orthogonal complement $H_ext$ , the external space, i.e., $H=H_0⊕H_ext$ . The model space is assumed to be a standard complete-active space (CAS), and for simplicity, we assume that no SPFs are inactive. The determinantal basis for $H$ is denoted $\{\ket{φ_ρ}\}$ , while the determinantal bases for $H_0$ and $H_ext$ are denoted $\{\ket{φ_α}\}$ and $\{\ket{φ_μ}\}$ , respectively. Thus, we will use the convention that $α,β,…$ denote model space detereminant indices, $μ,ν,…$ determine external-space indices, and $ρ,σ,…$ are of either type.
In the fermion case, a set of $N$ CAS functions constitute the formal reference, while in the boson case one single $φ_0⊂\{φ_p\}$ with occupation number $N$ determines the formal reference, in both cases denoted $\ket{φ_0}$ . We denote excitation operators relative to $\ket{φ_0}$ by $X_α$ (excitatons within the model space), and $X_μ$ (external excitations, i.e., excitations with at least one non-CAS label). In general, we allow $α=0$ in a model-space cluster operator, for which $X_α=I$ , the identity operator. A general excitation of either type is denoted $X_ρ$ . The algebra of general cluster operators $T=∑_ρτ_ρX_ρ$ is denoted $T$ , and the subalgebra of model-space excitation operators $T_0=∑_ατ_αX_α$ is denoted $T_0$ . The set of external cluster operators $T=∑_μτ_μX_μ$ is denoted $T_ext$ , and as linear spaces $T=T_0⊕T_ext$ .
We have the resolution of identity $I=P+Q$ with
$$
P=∑_α\ket{φ_α}\bra{φ_α}, Q=∑_μ\ket{
φ_μ}\bra{φ_μ}.
$$
Any $\ket{ψ}∈H$ such that $\braket{φ_0}{ψ}≠ 0$ can now be written uniquely as $\ket{ψ}=\ket{ψ_0}+\ket{ψ_ext}=(S_0+S)\ket{φ_0}$ , with $S_0∈T_0$ , and $S∈T_ext$ . Since any particle-number conserving cluster operator is nilpotent, there is a unique $T_0+T∈T$ such that $\exp(T_0+T)=\exp(T)\exp(T_0)=S_0+S$ . Let $C=\exp(T_0)∈T_0$ (since $T_0$ is a subalgebra). It follows that
$$
\ket{ψ}=e^T\ket{ψ_0}=e^TC\ket{φ_0}, \tag{1}
$$
a unique decomposition and the basic ansatz for CASCC, and at the current stage it is equivalent to untruncated standard single-reference CC theory.
Consider now the (partially) similarity transformed Schrödinger equation,
$$
e^-THe^T\ket{ψ_0}=E\ket{ψ_0}, \tag{2}
$$
where $H$ is the system Hamiltonian. Defining $\bar{H}_CAS=P\bar{H}P=Pe^-THe^TP=PHe^TP$ , and projecting onto CAS and external determinants we obtain a non-Hermitian CASCI eigenvalue problem and external CC amplitude equations, respectively,
$$
\displaystyle\bar{H}_CAS\ket{ψ_0} \displaystyle=E\ket{ψ_0}, \displaystyle\braket{φ_μ}{\bar{H}}{ψ_0} \displaystyle=0.
$$
Equations (3a) and (3b) are the linked CASCC working equations and fulfilled if and only if the exact Schrödinger equation $H\ket{ψ}=E\ket{ψ}$ is satisfied, as at this point, no cluster operator truncations have been introduced. We remark, that some implementations of semilinear CASCC employs the equivalent unlinked formulation where the Schrödinger equation is not similarity transformed, and where Eq. (3b) instead becomes $\braket{φ_μ}{(H-E)\exp(T)}{ψ_0}=0$ .
The standard CASCC approach to approximate treatments is to truncate the external amplitude spaces (for both $T$ and $Λ$ ) to a discrete space $T_d⊂T_ext$ that includes, say, singles, doubles, and the first-order interaction space [5]. Furthermore, $T_0$ (and $Λ_0$ ) are truncated such that it includes a selected set of higher-order excitations that together allow spanning a sufficiently flexible model vector. However, if the semilinear ansatz is considered, the CAS cluster operators must be untruncated. If not, there will be nonlinear constraints for the linear CAS amplitudes. (It is then typical to do a full similarity transform of $H$ , replacing Eqs. (II) with the fully linked SRCC equations $\braket{φ_ρ}{e^-T-T_0He^T+T_0}{φ_0}=Eδ_ρ,0$ .) Thus, we assume in the following that the CASCI vectors are kept untruncated. As in standard SRCC, for each amplitude $t_μ$ included in the approximate cluster operators $T∈T_d$ , the corresponding projective equation is included, constituting exactly one equation per unknown amplitude.
## III The first Lagrangian
We observe that for any $\bra{\tilde{ψ}_0}∈H_0^†$ , (equivalently, for any $\tilde{C}∈T_0^†$ , with $\bra{\tilde{ψ}_0}=\bra{φ_0}\tilde{C}$ ), such that $\braket{\tilde{ψ}_0}{ψ_0}≠ 0$ , we have
$$
E=E_\bar{H_CAS}(\tilde{ψ}_0,ψ_0)≡\frac{
\braket{\tilde{ψ}_0}{\bar{H}_CAS}{ψ_0}}{\braket{\tilde{ψ
}_0}{ψ_0}} \tag{4}
$$
if Eq. (3a) is fulfilled. Indeed, it is readily verified that the bivariate Rayleigh quotient $E_\bar{H_CAS}$ is stationary with respect to variations in $\bra{\tilde{ψ}_0}$ if and only if Eq. (3a) is satisfied, with the additional constraint that $\braket{\tilde{ψ}_0}{ψ_0}≠ 0$ . The amplitude equations (3b) are then additional constraints, leading to the introduction of a Lagrange multiplier cluster operator $Λ=∑_μλ_μX_μ^†∈T_ext^†$ ( $∈T_d^†$ in the truncated case) and a Lagrangian
$$
L≡\braket{\tilde{ψ}_0}{ψ_0}^-1(\bra{\tilde{ψ}_0}+\bra{
φ_0}Λ)\bar{H}\ket{ψ_0}. \tag{5}
$$
This is the first CASCC Lagrangian, providing a fully variational formulation whenever the model functions are full CASCI vectors. Indeed, it is readily verified, that optimization of $L$ with respect to $λ_μ$ and $\bra{\tilde{ψ}_0}$ reproduces the CASCC working equations (3a) and (3b), and that at the solution $L$ evaluates to the energy, regardless of the values of $Λ$ and $\bra{\tilde{ψ}_0}$ . Differentiation of $L$ with respect to the variables $\ket{ψ_0}$ and $τ_μ$ gives an inhomogenous CASCI eigenvalue problem and a linear system for $Λ$ , respectively,
$$
\displaystyle\bra{\tilde{ψ}_0}\bar{H}_CAS+\bra{φ_0}Λ
\bar{H}P \displaystyle=E\bra{\tilde{ψ}_0}, \displaystyle(\bra{\tilde{ψ}_0}+\bra{φ_0}Λ)[\bar{H},X_μ]
\ket{ψ_0} \displaystyle=0.
$$
It is interesting to note that Eq. (5) is on the form $L(ψ,\tilde{ψ})=\braket{\tilde{ψ}}{H}{ψ}/\braket{\tilde{ψ}}{ψ}$ , with
$$
\bra{\tilde{ψ}}=(\bra{\tilde{ψ}_0}+\bra{φ_0}Λ)e^-T. \tag{7}
$$
Since $e^-T$ is invertible, $\bra{\tilde{ψ}}∈H^†$ is completely general. This shows that the Lagrangian $L$ can be interpreted in terms of Arponen’s bivariational principle [15], and that the bra and ket will be exact left and right eigenvectors of $H$ in the full, untruncated limit. Moreover, the quantum mechanical state is assigned the density operator
$$
ρ=\braket{\tilde{ψ}}{ψ}^-1\ket{ψ}\bra{\tilde{ψ}}. \tag{8}
$$
When requiring the Hellmann–Feynman theorem to be valid, this implies that expectation values of observables are given by $\braket{Ω}=(ρΩ)=\braket{\tilde{ψ}}{Ω}{ ψ}/\braket{\tilde{ψ}}{ψ}$ whenever all variables $(T,Λ,ψ_0,\tilde{ψ}_0)$ are variationally determined as a critical point of $L$ .
The equations (6a) and (6b) are coupled, which in a practical calculation can be resolved as follows: Let $\bra{ψ^\prime_0}∈H_0^†$ be such that
$$
\bra{ψ^\prime_0}\bar{H}_CAS=E\bra{ψ^\prime_0},
$$
and we assume that the eigenvalue $E$ is isolated, and that $\braket{ψ_0^\prime}{ψ_0}=1$ . Introduce the projector $P_0=\ket{ψ_0}\bra{ψ^\prime_0}$ , and observe that $\bar{H}_CAS=EP_0+(1-P_0)\bar{H}_CAS(1-P_0)$ . We assume now that $E$ is the leftmost eigenvalue of $\bar{H}_CAS$ . Since $E$ is isolated, $\bar{H}_⊥≡(1-P_0)(\bar{H}_CAS-E)(1-P_0)$ has a right inverse as operator on the subspace $H_⊥^†=H_0^†(1-P_0)$ . We now decompose $\bra{\tilde{ψ}_0}=\bra{ψ^\prime_0}+\bra{χ}$ with $\bra{χ}∈H_⊥^†$ , so that $\bra{\tilde{ψ}_0}$ is normalized as $\braket{\tilde{ψ}_0}{ψ_0}=\braket{ψ^\prime_0}{ψ_0}=1$ . Inserting the decomposition into Eq. (6a) gives
$$
\bra{χ}=-\bra{φ_0}Λ\bar{H}(1-P_0)\bar{H}_⊥^-1. \tag{9}
$$
Equation (6b) becomes
$$
(\bra{ψ^\prime_0}-\bra{φ_0}Λ\bar{H}(1-P_0)\bar{H}_⊥^
-1+\bra{φ_0}Λ)[\bar{H},X_μ]\ket{ψ_0}=0, \tag{10}
$$
which is now a linear system for $Λ$ in terms of known quantities.
## IV Variational approximations of the first Lagrangian
Equation (5) describes a map $L:T_d×T^†_d×H_0× H_0^†→ℂ$ , whose critical points are equivalent to solutions of the left and right Schrödinger equations when $T_d=T_ext$ is not truncated. The CASCI vectors are kept untruncated, which incurs an exponential scaling of computational cost with respect to the system size, also for truncated cluster operators. We therefore introduce a differentiable manifold $M⊂H_0$ in which we wish to approximate $\ket{ψ_0}∈M$ and $\bra{\tilde{ψ}_0}∈M^†$ (where $\bra{ω}∈M^†$ if and only if $\ket{ω}∈M$ ), leading to a projected Lagrangian $L:T_d×T_d^†×M× M^†→ℂ$ . Relevant examples to keep in mind are the Hartree–Fock manifold, i.e., the set of Slater determinants defined by $N$ linearly independent SPFs, considered here for its simplicity, and the manifold of matrix-product states (MPSs) of fixed bond dimension [38, 39]. We will henceforth let $M$ be understood from context, and in particular the original exact case $M=H_0$ is a special case.
Near any $\ket{ψ_0}∈M$ we can choose local coordinates $z∈ V≡ℂ^n$ , where $n$ is the dimension of $M$ . That is, there is a smooth map $χ:V→M$ with a smooth inverse, with $\ket{ψ_0}=\ket{χ(z_0)}$ for some $z_0∈ V$ . (Strictly speaking, the map $χ$ may depend on the neighborhood of $ψ_0$ considered, as there may not be a global coordinatization of $M$ . However, this has no bearing on the results.) The tangent space $T_ψ_0M$ is the $n$ -dimensional subspace of $H_0$ spanned by the tangent vectors $\ket{t_k}≡\ket{∂_kχ(z_0)}$ , with $∂_k≡∂/∂_z_{k}$ . Without loss of generality, we may assume that these tangent vectors form an orthonormal system. Thus, the orthogonal projector onto $T_ψ_0M$ is $P_ψ_0≡∑_k\ket{t_k}\bra{t_k}$ . It follows from orthonormality that $\ket{∂_lt_k}=\ket{∂_l∂_kχ(z_0)}=(I-P_ψ_{ 0})\ket{∂_lt_k}$ . Similarly, we define local coordinates $\tilde{z}∈ V$ such that $\bra{\tilde{ψ}_0}=\bra{χ(\tilde{z}_0)}∈M^†$ , with corresponding tangent vectors $\bra{\tilde{t}_k}≡\bra{∂_kχ(\tilde{z}_0)}$ and projector $P_\tilde{ψ_0}=∑_k\ket{\tilde{t}_k}\bra{\tilde{t}_k}$ . We will assume without much loss of generality that $M$ contains rays, i.e., $\ket{ψ_0}∈M$ implies $c\ket{ψ_0}∈M$ for all scalars $c$ . This again implies that $\ket{ψ_0}∈ T_ψ_0M$ . For tangent vectors we will use the notation $\ket{v}=∑_kv_k\ket{t_k}$ . However, since $\ket{v}$ depends on its attachment point, a more proper notation would be $\ket{v;ψ_0}$ or similar. We will let it be clear from context whether a ket is in $M$ or in $T_ψ_0M$ .
We will shortly consider the stationary conditions of $L_M$ , for which we need to introduce arbitrary variations of of the model functions. We first take a detour, however, and consider the bivariational approximation of the left and right eigenvectors of an arbitrary non-Hermitian operator $K$ by vectors in $\ket{ψ_0}∈M$ and $\bra{ψ_0^\prime}∈M^†$ . The different notation for the bra compared to the above is deliberate. We introduce the bivariate Rayleigh quotient
$$
E_K(ψ_0,ψ_0^\prime)≡\frac{\braket{ψ_0^
\prime}{K}{ψ_0}}{\braket{{ψ}_0^\prime}{ψ_0}}. \tag{11}
$$
It is straightforward to compute its infinitesimal variation as
$$
\braket{{ψ}_0^\prime}{ψ_0}δE_K(ψ_0,{ψ}_0
^\prime)=\bra{δ{ψ}_0^\prime}(K-E)\ket{ψ_0}\\
+\bra{ψ_0^\prime}(K-E)\ket{δψ_0}, \tag{12}
$$
where $E=E_K(ψ_0,{ψ}_0^\prime)$ . In the exact case, where $M$ is the full linear space, $E_K$ is clearly stationary if and only if $K\ket{ψ_0}=E\ket{ψ_0}$ , $\bra{ψ_0^\prime}K=E\bra{ψ_0^\prime}$ , and $\braket{ψ_0^\prime}{ψ_0}≠ 0$ , which is an expression of the bivariational principle [15, 40], i.e., a variational formulation of the simultaneous solution of the left and right eigenvalue problem for $K$ . For the approximation on the manifold, $δE_K=0$ for all possible variations in the model functions if and only if $\braket{ψ_0^\prime}{ψ_0}≠ 0$ and
$$
\displaystyle P_{ψ_0^\prime}(K-E)\ket{ψ_0} \displaystyle=0, \displaystyle\bra{{ψ}_0^\prime}(K-E)P_ψ_0 \displaystyle=0.
$$
This is a coupled set of nonlinear equations for the bra and ket model functions. In the case of the Hartree–Fock manifold, Eq. (IV) becomes the non-Hermitian Hartree-Fock method of Froelich and Löwdin [41]. In the case of an MPS, one algorithm to solve Eq. (IV) is the alternating linear scheme (ALS) [39, 42] suitably generalized to independent bra and kets, i.e., to non-Hermitian problems. A similar generalization can be done for the modified ALS (MALS), which, for the non-Hermitian eigenvalue problem, collapses to the non-Hermitian density-matrix renormalization group method (DMRG) [12]. The DMRG algorithm, however, modifies the bond ranks of the MPS during its iterations, strictly speaking leaving the differentiable manifold world. This is a technical issue we will presently ignore. We further remark, that in typical situations, the tangent space projectors are not used on explicit form. Many manifolds $M$ of interest are also parameterized locally with coordinates having redundancies, a statement true for both the Hartree–Fock and MPS examples [43, 39]. This is resolved by introducing the mathematical structure of a principal bundle, and usually leads to straight-forward formulations of tangent-space equations [44, 28].
The stationary conditions $δ L=0$ are now straightforwardly computed as
$$
\displaystyle P_\tilde{ψ_0}(\bar{H}_CAS-E)\ket{ψ_0} \displaystyle=0 \displaystyle\bra{\tilde{ψ}_0}(\bar{H}_CAS-E)P_ψ_0 \displaystyle=-\bra{φ}Λ\bar{H}P_{ψ_0}. \displaystyle\braket{φ_μ}{\bar{H}}{ψ_0} \displaystyle=0 \displaystyle(\bra{\tilde{ψ}_0}+\bra{φ_0}Λ)[\bar{H},X_μ]
\ket{ψ_0} \displaystyle=0.
$$
Again, we note that the nonlinearities lead to couplings, such that in contrast to the standard CASCC equation (3a), the ket model vector cannot be solved independently of the bra model vector. Moreover any algorithm developed for the non-Hermitian eigenvalue problem, i.e., numerical solutions of Eq. (IV) such as the non-Hermitian DMRG algorithm, are not reusable for Eq. (IV) due to the inhomogeneity in Eq. (14b). This motivates the second Lagrangian to be presented in the next section.
## V The second Lagrangian
The second Lagrangian is based on the observation that, in the full CASCI version of CASCC, the bra model vector $\bra{\tilde{ψ}_0}$ can be chosen arbitrarily and still produce the right energy by projection, see Eq. (4). We therefore devise a scheme that forces it to be the variational left-eigenvector of $\bar{H}_CAS$ , i.e., we force both Eqs. (IV) to hold with $K=\bar{H}_CAS$ . Thus, we obtain additional constraints with Lagrange multipliers $\ket{u}∈ T_ψ_0M$ and $\bra{\tilde{u}}∈ T_\tilde{ψ_0}M^†$ , producing the Lagrangian
$$
\tilde{L}≡ L+∂_\tilde{ψ_0}E(ψ_0,\tilde{ψ}
_0;\tilde{u})+∂_ψ_0E(ψ_0,\tilde{ψ}_0;u), \tag{15}
$$
with $E=E_\bar{H_CAS}$ . The last two terms are directional derivatives given by
$$
∂_\tilde{ψ_0}E(ψ_0,\tilde{ψ}_0;\tilde{u})=
\braket{\tilde{ψ}_0}{ψ_0}^-1\bra{\tilde{u}}(\bar{H}_CAS-
E(ψ_0,\tilde{ψ}_0))\ket{ψ_0}, \tag{16}
$$
and
$$
∂_ψ_0E(ψ_0,\tilde{ψ}_0;u)=\braket{\tilde{
ψ}_0}{ψ_0}^-1\bra{\tilde{ψ}_0}(\bar{H}_CAS-
E(ψ_0,\tilde{ψ}_0))\ket{u}, \tag{17}
$$
respectively. In terms of the local coordinates, $\ket{u}=∑_k\ket{t_k}\braket{t_k}{u}=∑_k\ket{t_k}u_k$ and $\bra{\tilde{u}}=∑_k\tilde{u}_k\bra{\tilde{t}_k}$ , with $u,\tilde{u}∈ V$ being the $2n$ complex multipliers and the additional degrees of freedom introduced.
Differentiation with respect to the new Lagrange multipiers $\tilde{u}$ reproduces Eq. (14a), and differentiation with respect to $u$ produces the equation
$$
\bra{\tilde{ψ}_0}(\bar{H}_CAS-E)P_ψ_0=0, \tag{18}
$$
replacing Eq. (14b). The amplitude equations (14c) are still obtained by differentiation with respect to $λ_μ$ . Thus, a non-Hermitian left/right eigenvalue calculation is coupled to the CASCC external amplitude equations. Once solved, the ground-state energy is obtained straightforwardly.
The multipliers $Λ$ , $u$ , and $\tilde{u}$ are required if properties other than the ground-state energy are to be computed. The equations for the multipliers are obtained by by differentiation with respect to $T$ , $\bra{\tilde{ψ}_0}$ , and $\ket{ψ_0}$ , resulting in coupled linear systems, detailed in coordinate-free form in the Appendix.
## VI The third Lagrangian
The second Lagrangian $\tilde{L}$ solved the problem of the inhomogeneity of the left eigenvalue problem. However, the functional is not of the form of an expectation value. Some squinting at Eq. (15) allows us to come up with the following functional, depending on the CASCC variables in addition to multipliers $\ket{v}∈ T_ψ_0M$ and $\bra{\tilde{v}}∈ T_\tilde{ψ_0}M^†$ ,
$$
\hat{L}≡\frac{1}{2}\braket{φ_0}{Λ\bar{H}}{ψ_0}+\frac{1}{2
}\frac{\braket{\tilde{v}}{\bar{H}_CAS}{ψ_0}}{\braket{\tilde{v}}{
ψ_0}}+\frac{1}{2}\frac{\braket{\tilde{ψ}_0}{\bar{H}_CAS}{v}
}{\braket{\tilde{ψ}_0}{v}}. \tag{19}
$$
It is straightforward to show, that differentiation with respect to the tangent vectors reproduce Eqs. (14a) and (18) by means of the bivariational principle. Differentiation with respect to $λ_μ$ also reproduces Eq. (14c). The factors $1/2$ now ensure that the Lagrangian evaluates to the energy at the critical point, and that we have $\hat{L}=(Hρ)$ with the rank-2 non-Hermitian density operator $ρ=(ρ_1+ρ_2)/2$ , where
$$
ρ_1=\braket{\tilde{ψ_0}}{v}^-1e^T\ket{v}\bra{\tilde{ψ}_0}e^
{-T} \tag{20}
$$
and
$$
ρ_2=\braket{\tilde{v}}{ψ_0}^-1e^T\ket{ψ_0}(\bra{\tilde{v}}+
\braket{\tilde{v}}{ψ_0}\bra{φ_0}Λ)e^-T. \tag{21}
$$
It is readily verified, that $(ρ_i)=1$ and that $ρ_i^2=ρ_i$ .
Thus, the Lagragians $\tilde{L}$ and $\hat{L}$ are equivalent, in the sense that the working equations are identical for $\ket{ψ_0}$ , $\bra{\tilde{ψ}_0}$ , and $T$ , the values of which are sufficient to compute the energy using either Lagrangian. This also implies that the Hellmann–Feynman approach to expectation values and other properties yield identical results, even if the Lagrange multipliers have different values for the two Lagrangians.
In the Appendix, the Lagrange multiplier equations are derived, and they are closely related to those of the second Lagrangian $\tilde{L}$ , with basically the same computational footprint and details.
## VII Time-dependent formalism
Real-time propagation of coupled-cluster theory is gaining interest in the research community. It offers a computational alternative to response theory and excited states, and there is a demand for accurate first-principles simulations of molecular systems subject to intense and ultrashort laser pulses, see Refs. [45, 27] for reviews. In this section, we provide an overview of the explicitly time-dependent generalizations of the variational formulations. A detailed exposition including response theory based on this formalism is relegated to future work.
### VII.1 Action functional for the first Lagrangian
The time-dependent Schrödinger equation and its dual are the critical point conditions of the bivariate action functional [46, 15]
$$
S≡∫_t_0^t_1\braket{\tilde{ψ}(t)}{i∂_t-H}{ψ(t)}
dt, \tag{22}
$$
depending on the history of the system state for time $t∈[t_0,t_1]$ . As usual, any infinitesimal variation is assumed vanishing at the endpoints of the interval. Thus, $δ S=0$ with respect to variations in $\bra{\tilde{ψ}}$ and $\ket{ψ}$ if and only if, respectively,
$$
H\ket{ψ(t)}=i∂_t\ket{ψ(t)}, and \bra{\tilde{ψ}(t)}H
=-i∂_t\bra{\tilde{ψ}(t)}, \tag{23}
$$
in addition to boundary conditions. The functional $S$ generalizes the expectation-value functional to the time domain by treating the energy as generator for time evolution. Interestingly, the equations are canonical, i.e., on the form of classical Hamiltonian mechanics,
$$
i∂_t\ket{ψ}=∂_{\tilde{ψ}}K, -i∂_
t\bra{\tilde{ψ}}=∂_{ψ}K, \tag{24}
$$
with the Hamiltonian function $K=\braket{\tilde{ψ}}{H}{ψ}$ . We remark, that the overlap $\braket{\tilde{ψ}}{ψ}$ is preserved during time evolution.
The functional (22) is appropriate also for approximate states that contain rays, which is true for the CASCC bra and ket: For arbitrary complex scalars $c$ and $\tilde{c}$ , $\ket{ψ}→ c\ket{ψ}$ and $\bra{\tilde{ψ}}→\tilde{c}\bra{\tilde{ψ}}$ under the transformations $\ket{ψ_0}→ c\ket{ψ_0}$ and $(Λ,\bra{\tilde{ψ}_0})→(\tilde{c}Λ,\tilde{c}\bra{\tilde{ψ }_0})$ . Thus, $ρ$ in Eq. (8) is invariant. We can therefore assume that $ρ(t_0)=1$ , which will be preserved during evolution. For approximate methods that do not contain rays, the functional must be made explicitly invariant with respect to time-local phase and normalization by dividing the integrand by $\braket{\tilde{ψ}(t)}{ψ(t)}$ , as detailed for traditional variational methods in Ref. [47].
Inserting the CASCC bra and ket vectors from Eqs. (1) and (7), we obtain the functional
$$
S=∫_0^Ti\braket{φ_0}{Λ\dot{T}}{ψ_0}+i\braket{\tilde{
ψ}_0}{\dot{ψ}_0}-K(τ,λ,\tilde{ψ}_0,ψ_0)
dt, \tag{25}
$$
where the dot denotes a time derivative, and where
$$
K(T,Λ,ψ_0,\tilde{ψ}_0)=(\bra{φ_0}Λ+\bra{
\tilde{ψ}_0})\bar{H}\ket{ψ_0}, \tag{26}
$$
is the energy expressed in the CASCC variables, i.e., the Lagrangian (5) without the denominator. We are assuming that at all times $t$ , the model vectors are $\ket{ψ_0(t)}∈M$ , and $\bra{\tilde{ψ}_0(t)}∈M^†$ . Requiring $S$ to be stationary with respect to all independent variables now reveals equations of motion of the form
$$
\displaystyle iP_\tilde{ψ_0}∂_t\ket{ψ_0} \displaystyle=P_\tilde{ψ_0}\bar{H}\ket{ψ_0} \displaystyle iX(ψ_0)\dot{τ} \displaystyle=F(τ,ψ_0), \displaystyle-i(∂_t\bra{\tilde{ψ}_0})P_ψ_0 \displaystyle=-i\bra{φ_0}Λ\dot{T}P_ψ_0+(\bra{\tilde{ψ}_0
}+\bra{φ_0}Λ)\bar{H}P_ψ_0. \displaystyle-iX(ψ_0)^T\dot{λ} \displaystyle=G_0(τ,λ,ψ_0,\tilde{ψ}_0)+G(τ,λ,
ψ_0,\tilde{ψ}_0),
$$
the derivation of which is very similar to the stationary conditions for the Lagrangian $L$ . The vector-valued functions $F$ and $G$ are the Lagrangian derivatives with respect to $λ$ and $τ$ , respectively, under the condition that $\braket{\tilde{ψ}}{ψ}=1$ (which is preserved during evolution), i.e.,
$$
\displaystyle F(τ,ψ_0)_μ \displaystyle=\braket{φ_μ}{\bar{H}}{ψ_0}, \displaystyle G(τ,λ,ψ_0,\tilde{ψ}_0)_μ \displaystyle=(\bra{\tilde{ψ}_0}+\bra{φ_0}Λ)[\bar{H},X_μ]
\ket{ψ_0}. \tag{28}
$$
The vector-valued function $G_0$ is given by
$$
G_0(τ,λ,ψ_0,\tilde{ψ}_0)_μ=\braket{φ_0}{Λ X
_μJ\bar{H}}{ψ_0}. \tag{30}
$$
Here, we have introduced an oblique projection operator $J=J(ψ_0,\tilde{ψ}_0)$ discussed in the Appendix. The occurence of $\dot{T}$ in Eq. (27c) can be eliminated by means of Eq. (27b). In the case where $M=H_0$ , i.e., the full CASCI model spaces, $P_ψ_0=P_\tilde{ψ_0}=P_CAS$ . The equations of motion should be compared with the CASCC working equations. In particular, we note the non-Hermitian right time-dependent Schrödinger equation (27a) and an inhomogenous left equation (27c). We may “invert” the projector $P_\tilde{ψ_0}$ by means of $J$ (see the Appendix), so that Eq. (27a) becomes
$$
i∂_t\ket{ψ_0}=J\bar{H}\ket{ψ_0}, \tag{31}
$$
and Eq. (27c) becomes
$$
-i∂_t\bra{\tilde{ψ}_0}=-i\bra{φ_0}Λ\dot{T}J+(\bra{
\tilde{ψ}_0}+\bra{φ_0}Λ)\bar{H}J. \tag{32}
$$
We also note, that with a critical point of $L$ as initial condition, e.g., the ground state, the solutions to Eqs. (VII.1) are easily integrated to a stationary state given by $\ket{ψ(t)}=\exp(-iEt)\ket{ψ(0)}$ and $\bra{\tilde{ψ}(t)}=\exp(iEt)\bra{\tilde{ψ}(0)}$ . Linear perturbation theory will yield response-theory type results [48].
### VII.2 Action functional for the third Lagrangian
Turning to the third Lagrangian, the pure state is replaced by the density operator $ρ=(ρ_1+ρ_2)/2$ , with $ρ_i=\ket{ψ_i}\bra{\tilde{ψ}_i}$ being pure-state operators, see Eqs. (20) and (21). For a general convex combination $ρ=∑_kp_kρ_k$ of this form, where $0≤ p_k≤ 1$ and $∑_kp_k=1$ , we define an action functional
$$
\hat{S}=∑_k∫_0^Tp_k\bra{\tilde{ψ}_k}i∂_t-H\ket{
ψ_k} dt. \tag{33}
$$
Equations of motion for $ρ_i$ are now obtained by forcing $\hat{S}$ to be critical. It is readily seen, that if $E=(ρ H)$ is critical, then $ρ(t)=ρ$ gives a critical point of $\hat{S}$ , as expected. Moreover, if all the pure states are independent and allowed to vary freely in the full computational Hilbert space $H$ , then $\hat{S}$ is critical if and only if the left and right time-dependent Schrödinger equations for each pair $\ket{ψ_i}$ and $\bra{\tilde{ψ}_i}$ are satisfied. Finally, for approximate bra and ket wavefunctions whose normalization is not fixed, i.e., the manifold of $(ψ_i,\tilde{ψ}_i)$ contains rays, one can show that $\braket{\tilde{ψ}_i}{ψ_i}$ is preserved during evolution, and in particular that $(ρ)$ is a constant of motion. In conclusion, the action functional $\hat{S}$ is a reasonable starting point for defining approximate time evolution.
For the ansatz $ρ=(ρ_1+ρ_2)/2$ of the third Lagrangian, we may now assume $\braket{\tilde{ψ}_0}{v}=\braket{\tilde{v}}{ψ_0}=1$ , which yields
$$
\hat{S}=\frac{1}{2}∫_0^Ti\braket{φ_0}{Λ\dot{T}}{ψ_0}+i
\braket{\tilde{v}}{\dot{ψ}_0}+i\braket{\tilde{ψ}_0}{\dot{v}}+i
\braket{\tilde{ψ}_0}{\dot{P}_ψ_0}{v}-2\hat{K}, \tag{34}
$$
with
$$
\hat{K}=\frac{1}{2}≤ft(\braket{φ_0}{Λ\bar{H}}{ψ_0}
+\braket{\tilde{v}}{\bar{H}}{ψ_0}+\braket{\tilde{ψ}_0}{\bar{H}}{v}\right) \tag{35}
$$
being obtained from the Lagrangian $\hat{L}$ with the assumed normalizations. Using integration by parts, we obtain the equivalent action
$$
\hat{S}=\frac{1}{2}∫_0^Ti\braket{φ_0}{Λ\dot{T}}{ψ_0}-i
\braket{\dot{\tilde{v}}}{{ψ}_0}-i\braket{\dot{\tilde{ψ}}_0}{v}-i
\braket{\tilde{v}}{\dot{P}_\tilde{ψ_0}}{ψ_0}-2\hat{K}, \tag{36}
$$
indicating that the equations of motion of the bras and kets will exhibit a certain symmetry.
Using similar techniques as outlined in the Appendix, we can now derive equations of motion. Variation of the Lagrange multipliers $Λ$ and $\tilde{v}$ reproduces Eqs. (27a) and (27b). Variation of the Lagrange multiplier $v$ gives the equation
$$
\displaystyle-i\bra{\dot{\tilde{ψ}}_0}P_ψ_0 \displaystyle=\bra{\tilde{ψ}_0}\bar{H}P_ψ_0, T \displaystyle-iA(c)^T\dot{λ} \displaystyle=\frac{1}{2}G_0(τ,λ,ψ_0,\tilde{ψ}_0)+\hat{G}
(τ,λ,ψ_0,\tilde{ψ}_0,v,\tilde{v}),
$$
which replaces Eq. (27d). Here,
$$
\hat{G}_μ∂_τ_{μ}\hat{K}=\frac{1}{2}\braket{φ_0
}{\bar{H}_μ}{ψ_0}+\frac{1}{2}\braket{\tilde{v}}{\bar{H}_μ}{ψ_
{0}}+\frac{1}{2}\braket{\tilde{ψ}_0}{\bar{H}_μ}{v}, \tag{38}
$$
with $\bar{H}_μ=[\bar{H},X_μ]$ , see the Appendix.
It remains to determine equations of motion for the Lagrange multipliers $v$ and $\tilde{v}$ , obtained by varying $\bra{\tilde{ψ}_0}$ and $\ket{ψ_0}$ . Unlike the time-independent case, the multipliers are necessary to compute any expectation value, i.e., any physical prediction. The equations of motion, whose derivation is outlined in the Appendix, read
$$
\begin{split}-i\bra{\dot{\tilde{v}}}&=\bra{φ_0}Λ(\bar{H}-i\dot{T})J
+\bra{\tilde{v}}\bar{H}J\\
& -\bra{\tilde{ψ}_0}\bar{H}(JP^\prime_\tilde{ψ_0}(\tilde{v})
+[J,P^\prime_ψ_0(v)])J,\end{split} \begin{split}i\ket{\dot{v}}&=J\bar{H}\ket{v}-J(P^\prime_ψ_0(v)J+[P^
\prime_\tilde{ψ_0}(\tilde{v}),J])\bar{H}\ket{ψ_0}.\end{split}
$$
The occurence of $\dot{T}$ in Eq. (39a) can be eliminated using the amplitude equation of motion.
While these equations seem fairly complicated, they are expressed on coordinate-free form, and are expected to simplify considerably when concrete manifold parameterizations are employed. Moreover, if a manifold $M$ with low dimension compared to the full model CAS, the overall cost is small compared to evolution of the amplitude equations.
## VIII Possible applications
While the Lagrangians derived in this article should be useful for standard CASCC calculations, especially when properties and explicit time evolution are considered, the generalization to CASCC manifold approximation can have several additional uses. We here outline only a few ideas to be explored in future research.
### VIII.1 Correction of tailored CC
The tailored CC method [32] is an approach for external correction of a fixed variational model wavefunction. A variational model wavefunction $\ket{ψ_0}∈M⊂H_0$ is computed using, say CASCI as in Ref. [32], or the quantum chemical DMRG method, QC-DMRG [9, 10, 11, 49, 13], and inserted into the CC amplitude equations for an external cluster operator $T∈T_d$ , producing the DMRG-TCC method [33, 34, 35, 36]. An important caveat is that only CAS singles and doubles are used in practice due to the cost of extracting triples from the DMRG wavefunction. Thus, one first solves
$$
E_0=\braket{ψ_0,*}{H_CAS}{ψ_0,*}=\min≤ft\{\braket{ψ}
{H_CAS}{ψ}\midψ∈M\right\}, \tag{40}
$$
or equivalently,
$$
P_ψ_{0,*}(H_CAS-E_0)\ket{ψ_0,*}=0, \tag{41}
$$
extracts an approximation $\ket{ψ_0}≈\ket{ψ_0,*}$ due to practical considerations, and subsequently solves Eq. (3b) for the external $T∈T_d$ . Here, $H_CAS=PHP$ . The total energy is then given by
$$
E_TCC=E_0+E_corr=E_\bar{H_CAS}(ψ_
0,ψ_0),
$$
see Eq. (4). In many cases, it is the fully exponential parameterization which is actually implemented, i.e., $\ket{ψ_0}=\exp[T_0]\ket{φ_0}≈\ket{ψ_0,*}$ , with $T_0$ containing, say, only singles and doubles extracted from $\ket{ψ_0,*}$ .
Clearly, the TCC method can be viewed as a first approximation to the self-consistent solution of Eq. (IV), the stationary conditions of the first Lagrangian. Assume that $\ket{ψ_0}=\ket{ψ_0,*}$ , i.e., no approximations to the TCC CAS reference is done when inserting into the external amplitude equations. In order to improve upon the TCC calculation, then, it is easy to see that the current estimate $(\ket{ψ_0},T)$ must be fed back into the critical point condition and iterated until convergence. Note, however, that subsequent iterations involve both a non-hermitian right eigenvalue problem and an inhomogenous left eigenvalue problem.
Alternatively, one may resort to the second Lagrangian, where the inhomogenous left eigenvalue problem is replaced with a homogenous eigenvalue problem at the expense of additional Lagrange multipliers. However, if only the ground-state energy is required, these multipliers are not needed. A two-sided non-Hermitian eigenvalue problem can then be computed variationally, selfconsistently with the $T$ -amplitude equations, and the energy computed via
$$
E_TCC=E_\bar{H_CAS}(\tilde{ψ}_0,ψ_0),
$$
generalizing the previous formula. This approach allows existing implementations of non-Hermitian DMRG algorithms to be reused for corrected DMRG-TCC calculations. It must be stressed, however, that if additional properties need to be calculated, the Lagrange multipliers $Λ$ , $u$ and $\tilde{u}$ need to be computed.
The ideal approach, however, would be to use the first Lagrangian and adapt the DMRG algorithm to optimization of $L$ . The microiterations in the DMRG sweeps would be smaller-dimension versions of Eqs. (14a) and (14b), and would need novel computer implementations. We relegate a study of this algorithm for future study.
If some approximation to $\ket{ψ_0,*}$ is done, such as for current DMRG-TCC implementations, further considerations must be made. The semilinear CASCC Lagrangians still provide the exact result, and it is likely that correction schemes still can be found based on iterations or, e.g., perturbation theory.
### VIII.2 Attosecond dymanics
With the advent of technology producing ultra-short laser pulses, it has been possible to experimentally probe the laws of quantum mechanics at a time-scale resolving the motion of electrons in complex molecules [50]. In particular, even for an initial condition of single-reference nature, time evolution may produce strong multireference character. Moreover, description of unbound systems mandate the use of global resolution of single-particle basis functions in the form of grids, discrete-variable representations, or similar. Currently, the state-of-the-art methods for such descriptions with time-dependent single-particle bases functions are time-dependent multiconfigurational Hartree–Fock-type methods (MCTDHF) [51, 52, 53, 54, 55], or single-reference coupled-cluster-based approaches [21, 22, 23, 24]. The CASCC Lagrangians may be excellent starting points for development of low-scaling methods applicable to larger molecules and with higher-quality static correlation description than before.
Let us sketch the procedure for the derivation of the orbital-adaptive time-dependent version of the first Lagrangian, for simplicity restricting our attention to the case $M=H_0$ , i.e., we use the full, linear CASCI expansion for model vectors. The outline closely follows that of Ref. [21] in the general principles.
In an orbital-adaptive bivariational method such as we are going to describe for CASCC, the bra and ket are constructed from biorthogonal but otherwise independent sets $\{\tilde{φ}_x\}$ and $\{φ_x\}$ of SPFs, respectively, i.e., $\braket{\tilde{φ}_x}{φ_y}≡δ_xy$ , and correspondingly the creation and annihilation operators associated with these functions satisfy $\{\tilde{c}_x,c^†_y\}=δ_xy$ . For bosons, the anticommutator is replaced by a commutator. We stress, that the individual SPFs are typically resolved in a much larger computational basis, e.g., some kind of grid or discrete-variable representation. This is mandatory feature of real-time propagation methods whenever an unbound system of electrons are described, as a given finite and rather small basis set $\{φ_x\}$ and $\tilde{φ}_x$ is always localized,. Thus, the operator $P=∑_x\ket{φ_x}\bra{\tilde{φ}_x}$ projects onto the current single-particle basis span, and $Q=1-P$ projects onto the remaining degrees of freedom. Typically, the range of $Q$ has much higher dimension than the relatively compact single-particle space, and is not localized.
The SPF spaces are divided into occupied and unoccupied CAS functions $\{φ_i\}$ and $\{φ_a\}$ , respectively, and external SPFs $\{φ_α\}$ , i.e., $\{φ_x\}=\{φ_i\}∪\{φ_a\}∪\{φ_x\}$ , and similarly $\{\tilde{φ}_x\}=\{\tilde{φ}_i\}∪\{\tilde{φ}_a\}∪ \{\tilde{φ}_x\}$ . An infinitesimal change in the SPFs that respect the biorthogonality constraint now reads [21]
$$
\displaystyleδ\ket{φ_y} \displaystyle=∑_xU_yx\ket{φ_x}+\ket{χ_y}, \displaystyleδ\bra{\tilde{φ}_x} \displaystyle=-∑_y\bra{φ_y}{U}_yx+\bra{\tilde{χ}_x}, \tag{42}
$$
where $U$ is a matrix, and where $\ket{χ_y}$ and $\bra{\tilde{χ}_x}$ are arbitrary SPFs such that $P\ket{χ_y}=0$ and $\bra{\tilde{χ}_x}P=0$ . Thus, the infinitesimal variations can be decomposed into a part along the basis itself, and an “orthogonal” part. This orthogonal part is responsible for describing ionization dynamics, for example, as it allows the SPFs to travel away from the local region in space described by the current value for the of SPF basis.
The CASCC state $ρ=\ket{ψ}\bra{\tilde{ψ}}$ is invariant under arbitrary transformations of the individual collections of occupied $\{φ_i\}$ CAS SPFs, the unoccupied $\{φ_a\}$ CAS SPFs, and the external $\{φ_α\}$ SPFs. This means that the corresponding diagonal blocks in $U$ are redundant and can be set to arbitrary values, fixing gauge choices for the SPF variations. The remaining blocks and the functions $\bra{\tilde{χ}}$ and $\ket{χ}$ are independent variables of the variations, and when introduced into the action they will lead to equations of motion for $P\ket{φ_x}$ , $Q\ket{φ_x}$ , $\bra{\tilde{φ}_x}P$ , and $\bra{\tilde{φ}_x}Q$ .
The equations of motion are obtained as follows. First, we write $\ket{ψ}=e^TC\ket{φ_0}$ where $C∈T_0$ is a model-space cluster operator $C=∑_αc_αX_α$ , and observe that the time derivative can be separated in terms of amplitude derivatives and SPF derivatives as
$$
\frac{d}{dt}\ket{ψ}=(D_amp+D_SPF)\ket{ψ},
$$
with $D_amp=∑_α\dot{c}_α\frac{∂}{∂ c_α }+∑_μ\dot{τ}_μ\frac{∂}{∂τ_μ}$ and $D_SPF=∑_x\dot{c}^†_x\tilde{c}_x$ . When inserted into the action functional (22), we get
$$
S=∫_t_0^t_1i\braket{\tilde{ψ}}{\dot{ψ}_amp}-\braket{
\tilde{ψ}}{H-iD_SPF}{ψ} dt, \tag{44}
$$
where $\ket{\dot{ψ}_amp}=D_amp\ket{ψ}$ . We further mote note that $iD_SPF$ , often called the Coriolis operator, acts as a one-body operator shift of the Hamiltonian, simply due to the time-dependence of the basis. Derivation of equations of motion for the amplitudes of $T$ , $Λ$ , and $\ket{ψ_0}$ and $\bra{\tilde{ψ}_0}$ now follows as before. For the derivation of the equations of motion for the SPFs, we rewrite the action as
$$
S=∫_t_0^t_1i({ρ^(1)η})-∑_n=1,2
(ρ^(n)h^(n))+i\braket{\tilde{ψ}}{\dot{ψ}_
amp} dt, \tag{1}
$$
where $ρ^(n)$ are reduced density matrices, being functions of the the amplitudes only due to Wick’s theorem, and where $η=[η_xy]=[\braket{\tilde{φ}_x}{\dot{φ}_y}]$ , $h^(1)=[\braket{\tilde{φ}_x}{H^(1)}{φ_y}]$ and $h^(2)=[\braket{\tilde{φ}_x\tilde{φ}_x^\prime}{H^(2)}{ φ_yφ_y^\prime}]$ , the integrals of the one-plus-two-body Hamiltonian $H=H^(1)+H^(2)$ , and all these matrices are functions of the SPFs only. Finally, $\braket{\tilde{ψ}}{\dot{ψ}_amp}$ is a function of the amplitudes only – its variation is identically zero and it can be omitted from the action functional in this context. Variations of $\bra{\tilde{φ}_x}$ and $\ket{φ_y}$ on the form discussed above now straightforwardly produces (partially coupled) equations for $\ket{\dot{φ}_x}$ and $\bra{\dot{φ}_y}$ , mean-field type equations formally similar to MCTDHF equations of motion.
The linear CAS model vector expansions may of course be replaced by, say, time-dependent matrix-product states. In that case, one has to carefully analyze the SPF invariance properties of the resulting state before deriving equations of motion.
## IX Conclusion
We have introduced fully variational formulations of CASCC, in three different formulations: The standard CASCC method in its semilinear linked formulation, including a formulation where the CAS model functions are approximated on a smooth submanifold, and two equivalent yet different formulations, where the inhomogenous non-Hermitian left eigenvalue problem for the CAS problem is replaced by a homogenous eigenvalue problem. This has the benefit that it partially decouples nonlinear working equations and allows reuse of already existing codes for the non-Hermitian two-sided eigenvalue problem, but the drawback that it introduces additional variables, Lagrange multipliers living in tangent space of the CAS manifold.
The equations derived are fairly abstract, yet very general, and in actual applications, the CAS manifold, such as a matrix-product state, will have a quite complex definition and with a non-trivial description of the tangent space due to gauge degrees of freedom [28]. Thus, the equations are expected to take on a different form in concrete realizations.
As for future research, several interesting projects are immediate. From the theory side, a fully variational framework can be used to derive response theory and excited states. Another interesting line of research is to investigate the extension of the present formalism from CASCC to the somewhat more general approach of subsystem-embedding subalgebras by Kowalski [8].
On the practical side, the iterative correction of DMRG-TCC is an easy target, and can be approached by either the first or second/third Lagrangian, depending on whether one wants to reuse existing non-Hermitian DMRG implementations. One interesting approach is to generalize the DMRG algorithm to direct optimization of the first CASCC Lagrangian, by optimizing $L$ in sweeps, freezing all but uo to two blocks of the bra and ket MPSs in each microiteration. This approach would be a generalization of the alternating linear scheme (ALS) and the modified ALS algorithm [42, 39] to bivariational problems.
In this article, no attempt at evaluating the computational cost of methods like CASCC with matrix-product states/DMRG has been attempted. Such considerations are complex and depend on implementation details of both CC amplitude equations and the DMRG algorithm for non-Hermitian operators. It is possible that implementations will be challenging, and that pragmatic considerations will force through approximations not considered here. Moreover, even iterative corrections to DMRG-TCC using the Lagrangians may turn out to be prohibitively expensive. Even so, it is likely that some corrections can be extracted that are computationally affordable.
Acknowledgements. This work has received funding from the Research Council of Norway (RCN) under CoE Grant No 262695 (Hylleraas Centre for Quantum Molecular Sciences) and from ERC-STG-2014 under grant agreement No 639508 (BIVAQUM).
## Appendix A Working equations for the second Lagrangian
A general variation in $\ket{ψ_0}=\ket{χ(z_0)}∈M$ can be written
$$
\ket{δψ_0}=∑_k\ket{t_k}δ z_k, \tag{46}
$$
and similarly
$$
\bra{δ\tilde{ψ}_0}=∑_kδ\tilde{z}_k\bra{\tilde{t}_k}. \tag{47}
$$
General tangent vectors are $\ket{u}=∑_k\ket{t_k}u_k$ and $\bra{\tilde{u}}=∑_k\tilde{u}_k\bra{\tilde{t}_k}$ .
We begin by computing the variation of the second Lagrangian $\tilde{L}$ in Eq. (15) with respect to $\bra{\tilde{ψ}_0}$ . The variation can be split into the variation of $L$ and of the directional derivative terms in Eqs. (16) and (17). The variation of $L$ with respect to $\bra{\tilde{ψ}_0}$ is
$$
δ L=\braket{\tilde{ψ}_0}{ψ_0}^-1≤ft(\bra{δ\tilde{ψ}
_0}\bar{H}\ket{ψ_0}-\braket{δ\tilde{ψ}_0}{ψ_0}E
(\tilde{ψ}_0,ψ_0)\right)=0, \tag{48}
$$
due to Eq. (16) being zero at the solution. Next, we vary the constraint terms, where we note that any tangent vector $\bra{\tilde{u}}$ depends on $\bra{\tilde{ψ}_0}$ , and similarly $\ket{u}$ depends on $\ket{ψ_0}$ . This dependence can be taken into account by writing $\bra{\tilde{u}}=\bra{\tilde{u}}P_\tilde{ψ_0}$ and $\ket{u}=P_ψ_0\ket{u}$ , respectively, and hence we need to vary only the projection operators.
We are ready to compute the variations of Eq. (16) and (17) to obtain
$$
\braket{\tilde{ψ}_0}{ψ_0}δ\tilde{L}=\bra{\tilde{u}}δ P_
\tilde{ψ_0}(\bar{H}-E(\tilde{ψ}_0,ψ_0))\ket{ψ_0
}\\
+\braket{δ\tilde{ψ}_0}{(\bar{H}-E(\tilde{ψ}_0,ψ_0
))}{u}, \tag{49}
$$
and
$$
\braket{\tilde{ψ}_0}{ψ_0}δ(\tilde{L}-L)=\bra{\tilde{ψ}_0}(
\bar{H}-E(\tilde{ψ}_0,ψ_0))δ P_ψ_0\ket{u}\\
+\braket{\tilde{u}}{(\bar{H}-E(\tilde{ψ}_0,ψ_0))}{δ{
ψ}_0}. \tag{50}
$$
We have omitted additional terms that vanish at the critical point. We also compute
$$
δ L=\braket{φ_0}{Λ\bar{H}}{δψ_0},
$$
which does not vanish in general.
In order to resolve the variations in the projection operators, a small lemma is now useful: Let $P^\prime_ψ_0(u)$ be the directional derivative of $P_ψ_0$ in the direction $\ket{u}∈ T_ψ_0M$ . Then, for any $\ket{v}∈ T_ψ_0M$ ,
$$
P_ψ_0^\prime(u)\ket{v}=P_ψ_0^\prime(v)\ket{u}, \bra{\tilde{v}}P_\tilde{ψ_0}^\prime(\tilde{u})=\bra{\tilde{u}}P_
\tilde{ψ_0}^\prime(\tilde{v}).
$$
The proof is as follows: The directional derivative of the projector is
$$
P^\prime_{ψ_0}(u)=∑_klu_l\ket{∂_lt_k}\bra{t_k}+u_
{l}^*\ket{t_k}\bra{∂_lt_k}.
$$
It is straightforward to show, that any differential of a projection operator $P$ satisfies $δ P=Pδ PP+(1-P)δ P(1-P)$ , which implies that application of the directional derivative to $\ket{v}∈ T_ψ_0M$ eliminates the complex conjugates and gives
$$
\begin{split}P^\prime_ψ_0(u)\ket{v}&=∑_klu_lv_k\ket{∂
_l∂_kχ(z_0)}\\
&=∑_klu_lv_k\ket{∂_k∂_lχ(z_0)}=P^\prime_
ψ_0(v)\ket{u},\end{split}
$$
by symmetry of mixed partial derivatives. We now obtain the useful formula
$$
δ P_ψ_0\ket{v}=P^\prime_ψ_0(v)\ket{δψ_0},
$$
and a similar formula for the bra projector variation.
We now assume that the model functions are normalized as $\braket{\tilde{ψ}_0}{ψ_0}=1$ , and we set $E=E(\tilde{ψ}_0,ψ_0)$ , the eigenvalue of Eq. (3a) and Eq. (6a). We define $\bar{H}_E=\bar{H}-E$ . Equating Eqs. (49) and (50) to zero gives
$$
\displaystyle P_\tilde{ψ_0}[P^\prime_\tilde{ψ_0}(\tilde{u})
\bar{H}_E\ket{ψ_0}+\bar{H}_E\ket{u}] \displaystyle=0, \displaystyle[\bra{\tilde{ψ}_0}\bar{H}_EP^\prime_ψ_0(u)+\bra{
\tilde{u}}\bar{H}_E+\bra{φ_0}Λ\bar{H}]P_ψ_0 \displaystyle=0,
$$
We finally compute the variation in $\tilde{L}$ with respect to $T$ , which will give an equation linear in $Λ$ . We note that $∂_t_{μ}\bar{H}=\bar{H}_μ≡[\bar{H},X_μ]$ . Let $F_μ(x,y)≡E_\bar{H_μ}(x,y)=\braket{x}{\bar{H }_μ}{y}/\braket{x}{y}$ . It is straightforward to differentiate and get
$$
∂_t_{μ}\tilde{L}=\braket{φ}{Λ\bar{H}_μ}{ψ_0}+
F_μ(\tilde{ψ}_0,ψ_0)\\
+\braket{\tilde{u}}{ψ_0}(F_μ(\tilde{u},ψ_0)-F
_μ(\tilde{ψ}_0,ψ_0))\\
\braket{\tilde{ψ}_0}{u}(F_μ(\tilde{ψ}_0,u)-F
_μ(\tilde{ψ}_0,ψ_0)), \tag{53}
$$
which when equated to zero becomes a linear system for $Λ$ , to be compared with the $Λ$ equations for the standard CASCC equations, Eq. (6b).
Together, Eqs. (A) and (53) form a coordinate-independent linear system for $(\ket{u},\bra{\tilde{u}},Λ)$ . For the sake of concreteness, we expand Eq. (A) on component form to obtain
$$
\displaystyle Au+B\tilde{u} \displaystyle=0, \displaystyle A^T\tilde{u}+Cu+Dλ \displaystyle=0, \tag{54}
$$
where $A_kl=\braket{\tilde{t}_k}{\bar{H}_E}{t_l}$ and $B_kl=\braket{∂_k\tilde{t}_l}{\bar{H}_E}{ψ_0}$ , a symmetric matrix. Furthermore, $C_kl=\braket{\tilde{ψ}_0}{\bar{H}_E}{∂_kt_l}$ is symmetric, and $D_kμ=\braket{φ_μ}{\bar{H}}{t_k}$ , and we have collected the $Λ$ amplitudes in a vector as $λ=(λ_μ)$ .
## Appendix B Working equations for the third Lagrangian
Variation of $\hat{L}$ with respect to $\bra{\tilde{ψ}_0}$ yields
$$
2δ\hat{L}=\braket{\tilde{v}}{ψ_0}^-1\bra{\tilde{v}}δ P_
\tilde{ψ_0}(\bar{H}-E(\tilde{v},ψ_0))\ket{ψ_0}\\
+\braket{\tilde{ψ}_0}{v}^-1\bra{δ\tilde{ψ}_0}(\bar{H}-
E(\tilde{v},ψ_0))\ket{v}, \tag{56}
$$
while variation with respect to $\ket{ψ_0}$ yields
$$
2δ\hat{L}=\braket{φ_0}{Λ\bar{H}}{δψ_0}+\braket{
\tilde{v}}{ψ_0}^-1\bra{\tilde{v}}(\bar{H}-E(\tilde{v},ψ_0
))\ket{δψ_0}\\
+\braket{\tilde{ψ}_0}{v}^-1\bra{\tilde{ψ}_0}(\bar{H}-E(
\tilde{v},ψ_0))δ P_ψ_0\ket{v}. \tag{57}
$$
Since $\hat{L}$ is scale invariant with respect to $v$ and $\tilde{v}$ , we may choose $\braket{\tilde{v}}{ψ_0}=\braket{\tilde{ψ}_0}{v}=1$ , and it is readily seen that we obtain a linear system for $v$ and $\tilde{v}$ which is identical to that obtained for the second Lagrangian $\tilde{L}$ above, i.e., Eq. (A). However, we have not yet determined the equation for $Λ$ . Differentiation of $\hat{L}$ with respect to $t_μ$ gives
$$
2∂_t_{μ}\hat{L}=\braket{φ_0}{\bar{H}_μ}{ψ_0}+\mathcal
{F}_μ(\tilde{v},ψ_0)+F_μ(\tilde{ψ}_0,v), \tag{58}
$$
which is not the same as Eq. (53). Thus, Eq. (54), (55) (with $u=v$ and $\tilde{u}=\tilde{v}$ ), and Eq. (58) equated to zero form a linear system of equations for the Lagrange multiplier set $(Λ,v,\tilde{v})$ .
## Appendix C The operator $J$
Consider a tangent vector $\ket{x}∈ T_ψ_0M$ satisfying an equation
$$
P_\tilde{ψ_0}\ket{x}=P_\tilde{ψ_0}\ket{y}, \tag{59}
$$
where $\ket{y}∈H$ is arbitrary. We desire to solve for $\ket{x}$ . Similarly, we also consider a bra tangent vector $\bra{\tilde{x}}∈ T_\tilde{ψ_0}M^†$ satisfying the
$$
\bra{\tilde{x}}P_ψ_0=\bra{\tilde{y}}P_ψ_0. \tag{60}
$$
Starting with the first equation (59), projection onto $\bra{\tilde{t}_k}$ gives the equivalent condition $\braket{\tilde{t}_k}{x}=\braket{\tilde{t}_k}{y}$ . Introduce the overlap matrix $S_kl=\braket{\tilde{t}_k}{t_l}$ , and it is clear that for $\ket{x}$ to be unique, $S$ must be invertible. We expand $\ket{x}$ in the basis for $T_ψ_0M$ , i.e., $\ket{x}=∑_l\ket{t_l}\braket{t_l}{x}$ . Thus,
$$
\braket{\tilde{t}_k}{x}=∑_lS_kl\braket{t_l}{x}.
$$
Application of $S^-1$ from the left gives
$$
\braket{t_l}{x}=∑_kS^-1_lk\braket{\tilde{t}_k}{x}.
$$
Multiplication with $\ket{t_l}$ from the left and summing gives
$$
\ket{x}=J\ket{y},
$$
where
$$
J=∑_kl\ket{t_k}S^-1_kl\bra{\tilde{t}_l},
$$
which therefore acts as a solution operator to Eq. (59). A similar argument shows that the operator $J$ also acts as a solution operator for the bra equation (60), i.e., $\bra{\tilde{x}}=\bra{\tilde{y}}J$ .
It is straightforward to see that $J^2=J$ , but in general we have $J^†≠ J$ , so that $J$ is an oblique projector. The projector happens to be orthogonal whenever $\ket{ψ_0}=\ket{\tilde{ψ}_0}$ . We note that for any $\ket{y}∈H$ , $J\ket{y}∈ T_ψ_0M$ , and for any $\bra{\tilde{y}}∈H^†$ , we have $\bra{\tilde{y}}J∈ T_\tilde{ψ_0}M^†$ .
## Appendix D Equations of motion for the third Lagrangian
Starting from Eq. (36), variation of $ψ_0(t)∈M$ gives
$$
\begin{split}2δ\hat{S}&=∫_t_0^t_1i\braket{φ_0}{Λ
\dot{T}}{δψ_0}-i\braket{\dot{\tilde{v}}}{δ{ψ_0}}-i\braket{
\dot{\tilde{ψ}}_0}{δ P_ψ_0}{v}\\
& -i\braket{\tilde{v}}{\dot{P}_\tilde{ψ_0}}{δψ_0}-2
δ\hat{K} dt\\
&=∫_t_0^t_1i\braket{φ_0}{Λ\dot{T}}{δψ_0}-i
\braket{\dot{\tilde{v}}}{δ{ψ_0}}-i\braket{\dot{\tilde{ψ}}_0}{P^
{\prime}_ψ_0(v)}{δψ_0}\\
& -i\braket{\dot{\tilde{ψ}}_0}{P^\prime_\tilde{ψ_0}(\tilde{
v})}{δψ_0}-2δ\hat{K} dt\end{split} \tag{61}
$$
Further simplification reads
$$
\begin{split}2δ\hat{S}&=∫_t_0^t_1i\braket{φ_0}{Λ
\dot{T}}{δψ_0}-i\braket{\dot{\tilde{v}}}{δ{ψ_0}}\\
& +\braket{\tilde{ψ}_0}{\bar{H}J(P^\prime_\tilde{ψ_0}(
\tilde{v})+P^\prime_ψ_0(v))}{δψ_0}-2δ\hat{K}
dt\end{split} \tag{62}
$$
Here, we used Eqs. (A) and the equation of motion (37a). The variation in $\hat{K}$ is
$$
2δ\hat{K}=\braket{φ_0}{Λ\bar{H}}{δψ_0}+
\bra{\tilde{v}}\bar{H}\ket{δψ_0}+\bra{\tilde{ψ}_0}\bar{H}P^
\prime_ψ_0(v)\ket{δψ_0}. \tag{63}
$$
Since $δ\hat{S}=0$ is required for arbitrary $δψ(t)∈ T_ψ_0M$ , we obtain the Euler–Lagrange equations
$$
\begin{split}-i\bra{\dot{\tilde{v}}}&=\bra{φ_0}Λ(\bar{H}-i\dot{T})J
+\bra{\tilde{v}}\bar{H}J\\
& -\bra{\tilde{ψ}_0}\bar{H}(JP^\prime_\tilde{ψ_0}(\tilde{v})
+[J,P^\prime_ψ_0(v)])J.\end{split} \begin{split}i\ket{\dot{v}}&=J\bar{H}\ket{v}-J(P^\prime_ψ_0(v)J+[P^
\prime_\tilde{ψ_0}(\tilde{v}),J])\bar{H}\ket{ψ_0}.\end{split} \tag{34}
$$
## References
- [1] N. Oliphant and L. Adamowicz. Multireference coupled-cluster method using a single-reference formalism. J. Chem. Phys., 94:1229–1235, 1991.
- [2] N. Oliphant and L. Adamowicz. The implementation of the multireference coupled-cluster method based on the single-reference formalism. J. Chem. Phys., 96:3739–3744, 1992.
- [3] P. Piecuch, N. Oliphant, and L. Adamowicz. A state-selective multireference coupled- cluster theory employing the single-reference formalism. J. Chem. Phys., 99:1875–1900, 1993.
- [4] L.Z. Stolarczyk. Complete active space coupled-cluster method. extension of single-reference coupled-cluster method using the casscf wavefunction. Chem. Phys. Lett., 217(1):1–6, 1994.
- [5] D.I. Lyakh, M. Musiał, V.F. Lotrich, and R. J. Bartlett. Multireference Nature of Chemistry: The Coupled-Cluster View. Chem. Rev., 112:182–243, 2012.
- [6] L. Adamowicz and J.-P. Malrieu. Multireference self-consistent size-extensive state-selective configuration interaction. J. Chem. Phys., 105:9240–9247, 1996.
- [7] L. Adamowicz, J.-P. Malrieu, and V.V. Ivanov. New approach to the state-specific multireference coupled-cluster formalism. J. Chem. Phys., 112:10075–10084, 2000.
- [8] K. Kowalski. Properties of coupled-cluster equations originating in excitation sub-algebras. J. Chem. Phys., 148:094104, 2018.
- [9] S.R. White and R.L. Martin. Ab initio quantum chemistry using the density matrix renormalization group. J. Chem. Phys., 110:4127–4130, 1999.
- [10] Ö. Legeza, J. Röder, and B.A. Hess. QC-DMRG study of the ionic-neutral curve crossing of LiF. Mol. Phys., 101:2019–2028, 2003.
- [11] G.K.-L. Chan, M. Kállay, and J. Gauss. State-of-the-art density matrix renormalization group and coupled cluster theory studies of the nitrogen binding curve. J. Chem. Phys., 121(13):6110–6116, 2004.
- [12] G.K.-L. Chan and T. Van Voorhis. Density-matrix renormalization-group algorithms with nonorthogonal orbitals and non- Hermitian operators, and applications to polyenes. J. Chem. Phys., 122:204101, 2005.
- [13] S. Wouters and D. Van Neck. The density matrix renormalization group for ab initio quantum chemistry. Eur. Phys. J. D, 68:272, 2014.
- [14] Hendrik J. Monkhorst. Calculation of properties with the coupled-cluster method. Int. J. Quant. Chem., 12(S11):421–432, 1977.
- [15] J.S. Arponen. Variational principles and linked-cluster exp s expansions for static and dynamic many-body problems. Ann. Phys., 151(2):311–382, 1983.
- [16] L. Adamowicz, W. D. Laidig, and R. J. Bartlett. Analytical gradients for the coupled-cluster method. Int. J. Quant. Chem., 26(S18):245–254, 1984.
- [17] George Fitzgerald, Robert J. Harrison, and Rodney J. Bartlett. Analytic energy gradients for general coupled‐cluster methods and fourth‐order many‐body perturbation theory. J. Chem. Phys., 85(9):5143–5150, November 1986.
- [18] T. Helgaker and P. Jørgensen. Analytical Calculation of Geometrical Derivatives in Molecular Electronic Structure Theory. Adv. Quant. Chem., 19:183–245, 1988.
- [19] E. A. Salter, Gary W. Trucks, and Rodney J. Bartlett. Analytic energy derivatives in many‐body methods. I. First derivatives. J. Chem. Phys., 90(3):1752–1766, February 1989.
- [20] T. Helgaker and P. Jørgensen. Configuration-interaction energy derivatives in a fully variational formulation. Theor. Chim. Acta, 75:111–127, 1989.
- [21] S. Kvaal. Ab initio quantum dynamics using coupled-cluster. J. Chem. Phys., 136(19):194109, 2012.
- [22] T. Sato, H. Pathak, Y. Orimo, and K.L. Ishikawa. Communication: Time-dependent optimized coupled-cluster method for multielectron dynamics. J. Chem. Phys., 148:051101, 2018.
- [23] H. Pathak, T. Sato, and K.L. Ishikawa. Time-dependent optimized coupled-cluster method for multielectron dynamics. III. a second-order many-body perturbation approximation. J. Chem. Phys., 153:034110, 2020.
- [24] H. Pathak, T. Sato, and K.L. Ishikawa. Time-dependent optimized coupled-cluster method for multielectron dynamics. II. a coupled electron-pair approximation. J. Chem. Phys., 152:124115, 2020.
- [25] M.B. Hansen, N.K. Madsen, A. Zoccante, and O. Christiansen. Time-dependent vibrational coupled cluster theory: Theory and implementation at the two-mode coupling level. J. Chem. Phys., 151:154116, 2019.
- [26] T.B. Pedersen and S. Kvaal. Symplectic integration and physical interpretation of time-dependent coupled-cluster theory. J. Chem. Phys., 150:144106, 2019.
- [27] Benedicte Sverdrup Ofstad, Einar Aurbakken, Øyvind Sigmundson Schøyen, Håkon Emil Kristiansen, Simen Kvaal, and Thomas Bondo Pedersen. Time-dependent coupled-cluster theory. WIREs Comp. Mol. Sci., n/a(n/a):e1666, May 2023. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/wcms.1666.
- [28] C. Lubich, I.V. Oseledets, and B. Vandereycken. Time integration of tensor trains. SIAM J. Numer. Anal., 53:917–941, 2015.
- [29] L.-H. Frahm and D. Pfannkuche. Ultrafast ab initio quantum chemistry using matrix product states. J. Chem. Theory Comp., 15:2154–2165, 2019.
- [30] A. Baiardi and M. Reiher. Large-scale quantum dynamics with matrix product states. J. Chem. Theory Comp., 15:3481–3498, 2019.
- [31] Tilmann Bodenstein and Simen Kvaal. A state-specific multireference coupled-cluster method based on the bivariational principle. J. Chem. Phys., 153(2):024106, July 2020. Publisher: AIP Publishing.
- [32] T. Kinoshita, O. Hino, and R.J. Bartlett. Coupled-cluster method tailored by configuration interaction. J. Chem. Phys., 123(7):074106, 2005.
- [33] L. Veis, A. Antalík, J. Brabec, F. Neese, Ö Legeza, and J. Pittner. Coupled cluster method with single and double excitations tailored by matrix product state wave functions. J. Phys. Chem. Lett., 7:4072–4078, 2016.
- [34] F.M. Faulstich, A. Laestadius, Ö. Legeza, R. Schneider, and S. Kvaal. Analysis of the tailored coupled-cluster method in quantum chemistry. SIAM J. Numer. Anal., 57:2579–2607, 2019.
- [35] F.M. Faulstich, M. Máté, M.A. Laestadius, A.and Csirik, L. Veis, A. Antalik, J. Brabec, R. Schneider, J. Pittner, S. Kvaal, and Ö. Legeza. Numerical and theoretical aspects of the dmrg-tcc method exemplified by the nitrogen dimer. J. Chem. Theory Comp., 15:2206–2220, 2019.
- [36] J. Lang, A. Antalík, L. Veis, J. Brandejs, J. Brabec, Ö Legeza, and J. Pittner. Near-linear scaling in DMRG-based tailored coupled clusters: An implementation of DLPNO-TCCSD and DLPNO-TCCSD(t). J. Chem. Theory Comp., 16:3028–3040, 2020.
- [37] Maximilian Mörchen, Leon Freitag, and Markus Reiher. Tailored coupled cluster theory in varying correlation regimes. The Journal of Chemical Physics, 153(24):244113, December 2020.
- [38] S. Holtz, T. Rohwedder, and R. Schneider. On manifolds of tensors of fixed TT-rank. Numer. Math., 120:701–731, 2011.
- [39] M. Bachmayr, R. Schneider, and A. Uschmajew. Tensor networks and hierarchical tensors for the solution of high-dimensional partial differential equations. Foundations of Computational Mathematics, 16:1423–1472, 2016.
- [40] P.-O. Löwdin. On the stability problem of a pair of adjoint operators. J. Math. Phys., 24:70–87, 1983.
- [41] P. Froelich and P.-O. Löwdin. On the hartree–fock scheme for a pair of adjoint operators. J. Math. Phys., 24:88, 1983.
- [42] S. Holtz, T. Rohwedder, and R. Schneider. The alternating linear scheme for tensor optimization in the tensor train format. SIAM J. Sci. Comp., 34:A683–A713, 2012.
- [43] C. Lubich, T. Rohwedder, R. Schneider, and B. Vandereycken. Dynamical approximation by hierarchical Tucker ad tensor-train tensors. SIAM J. Matrix Anal. Appl., 34:476–494, 2015.
- [44] C. Lubich. From Quantum to Classical Molecular Dynamics: Reduced Models and Numerical Analysis. European Mathematical Society, 2008.
- [45] Xiaosong Li, Niranjan Govind, Christine Isborn, A. Eugene DePrince, and Kenneth Lopata. Real-Time Time-Dependent Electronic Structure Theory. Chem. Rev., 120(18):9951–9993, August 2020. Publisher: American Chemical Society (ACS).
- [46] P.R. Chernoff and J.E. Marsden. Properties of Infinite Dimensional Hamiltonian Systems. Springer, 1974.
- [47] P. Kramer and M. Saraceno. Geometry of the Time-Dependent Variational Principle in Quantum Mechanics, volume 140 of Lecture Notes In Physics. Springer, 1981.
- [48] R. McWeeny. Methods of molecular quantum mechanics. Academic Press, San Diego, San Francisco, New York, Boston, London, Sydney, Tokyo, 2nd edition, 1992.
- [49] Konrad Heinrich Marti and Markus Reiher. The Density Matrix Renormalization Group Algorithm in Quantum Chemistry. Zeitschrift für Physikalische Chemie, 224(3-4):583–599, April 2010. Publisher: De Gruyter (O).
- [50] K. Keisuke, Y. Ninota, and T. Sekikawa. Time-resolved high-harmonic spectroscopy of ultrafast photoisomerization dynamics. Opt. Expr., 26:31039, 2018.
- [51] J. Zanghellini, M. Kitzler, Ch. Fabian, T. Brabec, and A. Scrinzi. A MCTDHF approach to multi-electron dynamics in laser fields. Laser Physics, 13:1064–1068, 2003.
- [52] O.E. Alon and L.S. Streltsov, A.I. andCederbaum. Unified view on multiconfigurational time propagation for systems consisting of identical particles. J. Chem. Phys., 127:154103, 2007.
- [53] H. Miyagi and L.B. Madsen. Time-dependent restricted-active-space self-consistent-field theory for laser-driven many-electron dynamics. Phys. Rev. A, 87:062511, 2013.
- [54] T. Sato and K.L. Ishikawa. Time-dependent multiconfiguration self-consistent-field method based on the occupation-restricted multiple-active-space model for multielectron dynamics in intense laser fields. Phys. Rev. A, 91, 2015.
- [55] D. Hochstuhl, C.M. Hinz, and M. Bonitz. Time-dependent multiconfiguration methods for the numerical simulation of photoionization processes of many-electron atoms. The European Physical Journal Special Topics, 223(2):177–336, 2014.