{'fraction_non_alphanumeric': 0.04597274510843629, 'fraction_numerical': 0.010631528713081288, 'mean_word_length': 4.366960728514513, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Three Lagrangians for the complete-active space coupled-cluster method

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 couped-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. 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 densitymatrix 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 in the spirit of Helgaker and Jørgensen's Lagrangian formulation of traditional singlereference coupled-cluster (SRCC) theory [14,15] have been fruitless. The Lagrangian formulation of SRCC is today standard and allows easy access to response theory and molecular properties. Thus, as CASCC produces excellent absolute energies even for difficult multireference problems [16], the method is restricted to ground-state energies without a Lagrangian, since conversion of the * simen.kvaal@kjemi.uio.no semilinear form to standard SRCC exponential form can be impractical. 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 [17][18][19][20][21][22]. As the MPS ansatz is also being employed for realtime propagation of quantum chemical systems [23][24][25], one can speculate that a time-dependent CASCC method with a MPS-based model space description can be far superior to existing approaches.
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 lefteigenvalue 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 groundstate 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 {|φ ρ }, while the determinantal bases for H 0 and H ext are denoted {|φ α } and {|φ µ }, 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 |φ 0 . We denote excitation operators relative to |φ 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 Any |ψ ∈ H such that φ 0 |ψ = 0 can now be written uniquely as |ψ = |ψ 0 + |ψ ext = (S 0 + S) |φ 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 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, where H is the system Hamiltonian. DefiningH CAS = PHP = P e −T He T P = P e T P , and projecting onto CAS and external determinants we obtain a non-Hermitian CASCI eigenvalue problem and external CC amplitude equations, respectively, 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]. 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 ψ 0 | ∈ H † 0 , (equivalently, for anyC ∈ T † 0 , with ψ 0 | = φ 0 |C), such that ψ 0 |ψ 0 = 0, we have if Eq. (3a) is fulfilled. Indeed, it is readily verified that the bivariate Rayleigh quotient EH CAS is stationary with respect to variations in ψ 0 | if and only if Eq. (3a) is satisfied, with the additional constraint that ψ 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 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 ψ 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 ψ 0 |. Differentiation of L with respect to the variables |ψ 0 and τ µ gives an inhomogenous CASCI eigenvalue problem and a linear system for Λ, respectively, It is interesting to note that Eq. (5) is on the form L(ψ,ψ) = ψ |H|ψ / ψ |ψ , with Since e −T is invertible, ψ | ∈ H † is completely general. This shows that the Lagrangian L can be interpreted in terms of Arponen's bivariational principle [26], 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 When requiring the Hellmann-Feynman theorem to be valid, this implies that expectation values of observables are given by Ω = Tr(ρΩ) = ψ |Ω|ψ / ψ |ψ whenever all variables (T, Λ, ψ 0 ,ψ 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 and we assume that the eigenvalue E is isolated, and that ψ ′ 0 |ψ 0 = 1. Introduce the projector P 0 = |ψ 0 ψ ′ 0 |, and observe thatH CAS which is now a linear system for Λ in terms of known quantities.

IV. VARIATIONAL APPROXIMATIONS OF THE FIRST LAGRANGIAN
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 |ψ 0 ∈ M and ψ 0 | ∈ M † (where ω| ∈ M † if and only if |ω ∈ M), leading to a projected Lagrangian 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 [27,28]. 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 |ψ 0 ∈ M we can choose local coordinates z ∈ V ≡ C n , where n is the dimension of M. That is, there is a smooth map χ : V → M with a smooth inverse, with |ψ 0 = |χ(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 ψ0 M is the n-dimensional subspace of H 0 spanned by the tangent vectors |t k ≡ |∂ 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 ψ0 M is P ψ0 ≡ k |t k t k |. It follows from orthonormality that |∂ l t k = |∂ l ∂ k χ(z 0 ) = (I − P ψ0 ) |∂ l t k . Similarly, we define local coordinatesz ∈ V such that ψ 0 | = χ(z 0 )| ∈ M † , with corresponding tangent vectors t k | ≡ ∂ k χ(z 0 )| and projector Pψ 0 = k |t k t k |. We will assume without much loss of generality that M contains rays, i.e., |ψ 0 ∈ M implies c |ψ 0 ∈ M for all scalars c. This again implies that |ψ 0 ∈ T ψ0 M. For tangent vectors we will use the notation |v = k v k |t k . However, since |v depends on its attachment point, a more proper notation would be |v; ψ 0 or similar. We will let it be clear from context whether a ket is in M or in T ψ0 M.
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 |ψ 0 ∈ M and ψ ′ 0 | ∈ M † . The different notation for the bra compared to the above is deliberate. We introduce the bivariate Rayleigh quotient It is straightforward to compute its infinitesimal variation as This is a coupled set of nonlinear equations for the bra and ket model functions. In the case of the Hartree-Fock manifold, Eq. (13) becomes the non-Hermitian Hartree-Fock method of Froelich and Löwdin [30]. In the case of an MPS, one algorithm to solve Eq. (13) is the alternating linear scheme (ALS) [28,31] 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 densitymatrix 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 [28,32]. This is resolved by introducing the mathematical structure of a principal bundle, and usually leads to straight-forward formulations of tangentspace equations [23,33]. The stationary conditions δL = 0 are now straightforwardly computed as 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. (13) such as the non-Hermitian DMRG algorithm, are not reusable for Eq. (14) 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 ψ 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 lefteigenvector ofH CAS , i.e., we force both Eqs. (13) to hold with K =H CAS . Thus, we obtain additional constraints with Lagrange multipliers |u ∈ T ψ0 M and ũ| ∈ Tψ 0 M † , producing the Lagrangiañ with E = EH CAS . The last two terms are directional derivatives given by and respectively. In terms of the local coordinates, |u = k |t k t k |u = k |t k u k and ũ| = kũ k t k |, with u,ũ ∈ V being the 2n complex multipliers and the additional degrees of freedom introduced.
Differentiation with respect to the new Lagrange multipiersũ reproduces Eq. (14a), and differentiation with respect to u produces the equation 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ũ 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 , ψ 0 |, and |ψ 0 , resulting in coupled linear systems, detailed in coordinate-free form in the Appendix.

VI. THE THIRD LAGRANGIAN
The second LagrangianL 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 |v ∈ T ψ0 M and ṽ| ∈ Tψ 0 M † , It is straightforward to show, that differentiation with respect to the tangent vectors reproduce Eqs. (14a) and (14b') 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 havê L = Tr(Hρ) with the rank-2 non-Hermitian density operator ρ = (ρ 1 + ρ 2 )/2, where and It is readily verified, that Tr(ρ i ) = 1 and that ρ 2 i = ρ i . Thus, the LagragiansL andL are equivalent, in the sense that the working equations are identical for |ψ 0 , ψ 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 LagrangianL, with basically the same computational footprint and details.

VII. TIME-DEPENDENT FORMALISM
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.

A. 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 [26,34] 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 ψ | and |ψ if and only if, respectively, 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, with the Hamiltonian function K = ψ |H|ψ . We remark, that the overlap ψ |ψ is preserved during time evolution.
The functional (21) is appropriate also for approximate states that contain rays, which is true for the CASCC bra and ket: For arbitrary complex scalars c andc, |ψ → c |ψ and ψ | →c ψ | under the transformations |ψ 0 → c |ψ 0 and (Λ, ψ 0 |) → (cΛ,c ψ 0 |). Thus, ρ in Eq. (8) is invariant. We can therefore assume that Tr ρ(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 ψ (t)|ψ(t) , as detailed for traditional variational methods in Ref. [35].
Inserting the CASCC bra and ket vectors from Eqs. (1) and (7), we obtain the functional where the dot denotes a time derivative, and where 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 |ψ 0 (t) ∈ M, and ψ 0 (t)| ∈ M † . Requiring S to be stationary with respect to all independent variables now reveals equations of motion of the form 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 ψ |ψ = 1 (which is preserved during evolution), i.e., The vector-valued function G 0 is given by Here, we have introduced an oblique projection operator J = J(ψ 0 ,ψ 0 ) discussed in the Appendix. The occurence ofṪ in Eq. (26c) can be eliminated by means of Eq. (26b). In the case where M = H 0 , i.e., the full CASCI model spaces, P ψ0 = Pψ 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 (26a) and an inhomogenous left equation (26c). We may "invert" the projector Pψ 0 by means of J (see the Appendix), so that Eq. (26a) becomes i∂ t |ψ 0 = JH |ψ 0 , (26a') and Eq. (26c) becomes We also note, that with a critical point of L as initial condition, e.g., the ground state, the solutions to Eqs. (26) are easily integrated to a stationary state given by |ψ(t) = exp(−iEt) |ψ(0) and ψ (t)| = exp(iEt) ψ (0)|. Linear perturbation theory will yield response-theory type results [36].

B. 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 = |ψ i ψ i | being pure-state operators, see Eqs. (19) and (20). For a general convex combination ρ = k p k ρ k of this form, where 0 ≤ p k ≤ 1 and k p k = 1, we define an action functional Equations of motion for ρ i are now obtained by forcinĝ S to be critical. It is readily seen, that if E = Tr(ρH) is critical, then ρ(t) = ρ gives a critical point ofŜ, as expected. Moreover, if all the pure states are independent and allowed to vary freely in the full computational Hilbert space H, thenŜ is critical if and only if the left and right time-dependent Schrödinger equations for each pair |ψ i and ψ i | are satisfied. Finally, for approximate bra and ket wavefunctions whose normalization is not fixed, i.e., the manifold of (ψ i ,ψ i ) contains rays, one can show that ψ i |ψ i is preserved during evolution, and in particular that Tr(ρ) is a constant of motion. In conclusion, the action functionalŜ is a reasonable starting point for defining approximate time evolution. For the ansatz ρ = (ρ 1 + ρ 2 )/2 of the third Lagrangian, we may now assume ψ 0 |v = ṽ|ψ 0 = 1, which yieldŝ being obtained from the LagrangianL with the assumed normalizations. Using integration by parts, we obtain the equivalent action ) 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ṽ reproduces Eqs. (26a) and (26b). Variation of the Lagrange multiplier v gives the equation which replaces Eq. (26c), in a similar manner as for the time-independent case. Continuing, variation of T gives the equation which replaces Eq. (26d). Here, It remains to determine equations of motion for the Lagrange multipliers v andṽ, obtained by varying ψ 0 | and |ψ 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 and The occurence ofṪ in Eq. (35a) 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.

A. Correction of tailored CC
The tailored CC method [37] is an approach for external correction of a fixed variational model wavefunction. A variational model wavefunction |ψ 0 ∈ M ⊂ H 0 is computed using, say the quantum chemical DMRG method, QC-DMRG [9][10][11]13], and inserted into the CC amplitude equations for an external cluster operator T ∈ T d , producing the DMRG-TCC method [38][39][40][41]. Thus, one first solves or equivalently, and subsequently solves Eq. (3b) for the external T ∈ T d .
Here, H CAS = P HP . The total energy is then given by see Eq. (4). Clearly, the TCC method can be viewed as a first approximation to the self-consistent solution of Eq. (14), the stationary conditions of the first Lagrangian. In order to improve upon the TCC calculation, however, it is easy to see that the current estimate (|ψ 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 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ũ 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.

B. 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 [42]. In particular, even for an initial condition of singlereference 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-Focktype methods (MCTDHF) [43][44][45][46][47], or single-reference coupled-cluster-based approaches [17][18][19][20]. 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. [17] 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 {φ x } and {ϕ x } of SPFs, respectively, i.e., φ x |ϕ y ≡ δ xy , and correspondingly the creation and annihilation operators associated with these functions satisfy {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φ x is always localized,. Thus, the operator P = x |ϕ x φ x | projects onto the current singleparticle 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., An infinitesimal change in the SPFs that respect the biorthogonality constraint now reads [17] δ |ϕ y = where U is a matrix, and where |χ y and χ x | are arbitrary SPFs such that P |χ y = 0 and χ 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 ρ = |ψ ψ | 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 χ| and |χ are independent variables of the variations, and when introduced into the action they will lead to equations of motion for P |ϕ x , Q |ϕ x , φ x | P , and φ x | Q.
The equations of motion are obtained as follows. First, we write |ψ = e T C |φ 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 where |ψ amp = D amp |ψ . We further mote note that iD SPF , often called the Coriolis operator, acts as a onebody 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 |ψ 0 and ψ 0 | now follows as before. For the derivation of the equations of motion for the SPFs, we rewrite the action as where ρ (n) are reduced density matrices, being functions of the the amplitudes only due to Wick's theorem, and (2) |ϕ y ϕ y ′ ], the integrals of the oneplus-two-body Hamiltonian H = H (1) + H (2) , and all these matrices are functions of the SPFs only. Finally, ψ |ψ 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 φ x | and |ϕ y on the form discussed above now straightforwardly produces (partially coupled) equations for |φ x and φ 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 [23]. 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 [28,31]  A general variation in |ψ 0 = |χ(z 0 ) ∈ M can be written and similarly General tangent vectors are |u = k |t k u k and ũ| = kũ k t k |. We begin by computing the variation of the second La-grangianL in Eq. (15) with respect to ψ 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 ψ 0 | is δL = ψ 0 |ψ 0 −1 δψ 0 |H |ψ 0 − δψ 0 |ψ 0 E(ψ 0 , ψ 0 ) = 0, (A3) due to Eq. (16) being zero at the solution. Next, we vary the constraint terms, where we note that any tangent vector ũ| depends on ψ 0 |, and similarly |u depends on |ψ 0 . This dependence can be taken into account by writing ũ| = ũ| Pψ 0 and |u = P ψ0 |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 and We have omitted additional terms that vanish at the critical point. We also compute δL = φ 0 |Λ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 ′ ψ0 (u) be the directional derivative of P ψ0 in the direction |u ∈ T ψ0 M. Then, for any |v ∈ T ψ0 M, and by symmetry, The proof is as follows: The directional derivative of the projector is It is straightforward to show, that any differential of a projection operator P satisfies δP = P δP P + (1 − P )δP (1 − P ), which implies that application of the directional derivative to |v ∈ T ψ0 M eliminates the complex conjugates and gives by symmetry of mixed partial derivatives. We now obtain the useful formula δP ψ0 |v = P ′ ψ0 (v) |δψ 0 , and a similar formula for the bra projector variation.
Appendix C: The operator J Consider a tangent vector |x ∈ T ψ0 M satisfying an equation where |y ∈ H is arbitrary. We desire to solve for |x . Similarly, we also consider a bra tangent vector x| ∈ Tψ 0 M † satisfying the x| P ψ0 = ỹ| P ψ0 .
Starting with the first equation (C1), projection onto t k | gives the equivalent condition t k |x = t k |y . Introduce the overlap matrix S kl = t k |t l , and it is clear that for |x to be unique, S must be invertible. We expand |x in the basis for T ψ0 M, i.e., |x = l |t l t l |x . Thus, t k |x = l S kl t l |x . A similar argument shows that the operator J also acts as a solution operator for the bra equation (C2), i.e., x| = ỹ| J.
Since δŜ = 0 is required for arbitrary δψ(t) ∈ T ψ0 M, we obtain the Euler-Lagrange equations