# Second response theory: A theoretical formalism for the propagation of
quantum superpositions

Martín A. Mosquera  Department of Chemistry and Biochemistry, Montana State
University, Bozeman, MT 59717, USA<EMAIL_ADDRESS>
###### Abstract

The propagation of general electronic quantum states provides information of
the interaction of molecular systems with external driving fields. These can
also offer understandings regarding non-adiabatic quantum phenomena. Well
established methods focus mainly on propagating a quantum system that is
initially described exclusively by the ground state wavefunction. In this
work, we expand a previously developed size-extensive formalism within coupled
cluster theory, called second response theory, so it propagates quantum
systems that are initially described by a general linear combination of
different states, which can include the ground state, and show how with a
special set of time-dependent cluster operators such propagations are
performed. Our theory shows strong consistency with numerically exact results
for the determination of quantum mechanical observables, probabilities, and
coherences. We discuss unperturbed non-stationary states within second
response theory and their ability to predict matrix elements that agree with
those found in linear and quadratic response theories. This work also
discusses an approximate regularized methodology to treat systems with
potential instabilities in their ground-state cluster amplitudes, and compares
such approximations with respect to reference results from standard unitary
theory.

##  1 Introduction

The quantitative prediction of the behavior of quantum systems is a critical
task in electronic structure theory and a main motivator to find improved
algorithms and approaches for advanced computing architectures [1, 2, 3], such
as quantum computers [4, 5, 6, 7, 8]. The field of electronic structure has
advanced significantly in recent years where it is now possible to perform
routinary molecular drug design simulations [9, 10], materials high-throughput
screening [11, 12], simulations of biochemical systems [13, 14], among others.
In the field of materials research, for instance, localized states in
semiconductors that arise from impurities can behave as pseudo-atomic states
[15] that are excited by thermal fluctuations or related phenomena. Because of
systems like these, a constant target of theory development focuses on
computation of energies and forces. However, other properties are also crucial
for the engineering of quantum systems [16]. This includes variables such as
electrical and magnetic multipolar moments, and their related derivatives. A
more advanced type of computation is the direct propagation of a quantum
system that involves excited states in its initial state. This is an important
goal since non-linear propagations reveal behaviors that are inaccessible by
the standard analysis of energetic properties. Real time propagations may also
be more computationally efficient than techniques based on
eigenvalues/eigenvectors only [17, 18], and offer the advantage of being able
to capture the response of the system to strong and ultra-strong driving
fields, as well as their non-adiabatic couplings to nuclear motion [19].

There is a wide set of electronic structure theories available for different
types of applications [20, 21, 22, 23, 24, 25]. Among these, coupled cluster
(CC) theory is a very reliable formalism to study quantum systems [26, 27,
28]. CC methods enjoy the ability to predict accurate quantities that feature
size-extensivity [29, 30, 31, 32]; meaning they, in principle, scale correctly
to the macroscopic scale. They are systematically improvable [33], so they are
guaranteed to reach exact numerical values if computational resources allow
it. These points have been strong motivators for the scientific community to
explore the potential application of CC methods via quantum computers. Unitary
CC theory [34, 35, 36, 37, 38], for example, has gained significant attention
for predicting ground-state energies through quantum algorithms, with
promising results [39, 40, 41, 42]. Therefore, it can be expected that, as
quantum and standard architectures continue to advance, the applicability of
CC-based theory and algorithms will expand as well.

Ground-state methods in CC theory have reached a highly productive stage for
the design and understanding of molecular systems. Algorithms to study
periodic systems have also been deployed recently [43, 44]. Excited-state
methods, on the other hand, are the subject of on-going theory development.
Furthermore, time-dependent (TD) CC formalisms are being designed and
implemented to understand quantum systems beyond equilibrium [45, 46, 47, 48,
49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62]. This field started
with early theoretical work that explored the potential use of time-dependent
CC operators for the calculation of observable quantities in electronic
structure and nuclear dynamics [63, 64, 65, 66, 67], and it has provided
support to such formalisms that continue making progress [46]. These have made
significant advances for propagations that start from the ground state, and in
combination with non-Hermitian (standard) CC, variational CC, or unitary CC.
There are opportunities, however, to develop formalisms that predict the
evolution of general linear combinations of quantum states (which may or may
not include the ground state of the system). As common goals with standard CC
methods, extended TD techniques must be based on linked diagrammatic
expressions and systematic improvability. Quantum superpositions of ground and
excited states are fundamental in the study of coherent phenomena. These take
place in molecular and materials excitonics, quantum computers, ultracold
atoms, quantum spectroscopy, nuclear magnetic and electron paramagnetic
resonance experiments, spin lattices, among other physically relevant cases.
For these reasons, theoretical techniques capable of predicting time-dependent
and frequency-space observable signatures of coherent phenomena are critical.
Although we focus in this work on electronic degrees of freedom, the
techniques discussed here can be extended to the study of coherent quantum
vibrations, or vibronic states.

Motivated by the points above, a TD CC approach was proposed in Ref. [68].
This formalism, called second response (SR) theory, can be used to predict the
time-dependent behavior of systems that are initially in a linear combination
of excited states, while maintaining the desired feature of diagrammatic
connectedness (ground states are included through an extra set of equations).
We showed the formalism to be consistent with linear and quadratic CC response
theories, and validated it successfully with model systems [68]. SR theory is
an extension of previous work focused on the linear response of excited states
[69, 70, 71]. In this work, we show further developments of SR theory in
which, by using three TD excitation vectors, the expectation value of a
quantum observable as a function of time can be predicted, and this includes
probabilities and coherences. Our theory propagates general quantum
superpositions that can now include the ground state of the system. We then
show the solution of the SR equations that give cluster operators that
describe unperturbed non-stationary states. Finally, we discuss an
approximated method to treat systems that require regularization to achieve
stability. These developments are then verified through a three-
electron/three-level system that can be solved exactly using standard unitary
quantum propagations.

##  2 Theory

###  2.1 Motivation

This work focuses on non-relativistic Coulombic electronic systems. The free
(or unperturbed) Hamiltonian, denoted as $\hat{H}_{0}$, is the sum of the
kinetic, repulsion, electron-nuclear attraction (assuming fixed nuclei)
components. We are interested in electronic systems that are subject to an
external perturbation, so the TD Hamiltonian of interest reads
$\hat{H}(t)=\hat{H}_{0}+\hat{V}(t)$, where $\hat{V}(t)$ is the TD perturbing
potential operator. It is important to distinguish two viewpoints: First,
standard TD quantum mechanics provides expressions and rules to determine a TD
observable of interest. We denote the initial state of the system in the
standard quantum mechanical picture as $\Psi(0)$. Second, our goal is to
develop coupled-cluster equations that propagate the CC analogue of $\Psi(0)$.
So throughout this work we will establish connections between standard quantum
mechanical observables and their CC equivalents. For convenience we will often
use Heisenberg operators: If $\hat{A}$ denotes the operator corresponding to
some observable $A$ of interest, then
$\hat{A}^{\rm{H}}(t)=\hat{U}^{\dagger}(t)\hat{A}\hat{U}(t)$, where
$\hat{U}(t)=\mathcal{T}\exp[-\rmi\int_{0}^{t}\hat{H}(s)\rmd s]$ (with
$\mathcal{T}$ denoting the time-ordering superoperator).

We define the auxiliary wavefunction:

$|\Psi_{\rm{R}}(t=0;g_{\rm{R}})\rangle=|\Psi_{0}\rangle+g_{\rm{R}}|\Psi(0)\rangle~{},$
(1)

where:

$|\Psi(0)\rangle=S|\Psi_{0}\rangle+\sum_{N\neq 0}C_{N}|\Psi_{N}\rangle~{},$
(2)

$\Psi_{0}$ denotes the ground state of the system, $g_{\rm{R}}$ a coefficient,
and $\Psi_{N}$ an excited state ($N\neq 0$), the sum ($\sum_{N}$) runs over
excited states. We also define the wavefunction: $\Psi_{\rm{E}}=\sum_{N\neq
0}C_{N}\Psi_{N}$. The wavefunctions $\\{\Psi_{J}\\}$ satisfy
$\hat{H}_{0}|\Psi_{J}\rangle=E_{J}|\Psi_{J}\rangle$. $\Psi(0)$ is the initial
state of interest to be propagated. The exact TD wavefunction of the system
reads $|\Psi(t)\rangle=\hat{U}(t)|\Psi(0)\rangle$. If known, $\Psi(t)$ would
determine all the measurable quantities of the system. Previous work [68] on
SR CC theory showed how to propagate general arbitrary states by computing
separately the CC representations of
$\langle\Psi_{\rm{E}}|\hat{A}^{\rm{H}}(t)|\Psi_{0}\rangle$, and
$\langle\Psi_{\rm{E}}|\hat{A}^{\rm{H}}(t)|\Psi_{\rm{E}}\rangle$. In this
current work we show how, through an integrated approach, an arbitrary initial
state that includes a contribution from the ground state is propagated in SR
CC theory. Both approaches are equivalent and give the same numerical results,
but the present formulation may offer advantages for potential
implementations.

SR theory is motivated by the following property of $\Psi_{\rm{R}}$ in
standard TD quantum mechanics:

$\lim_{g_{\rm{R}}\rightarrow 0}\frac{\partial}{\partial
g_{\rm{R}}}|\Psi_{\rm{R}}(0;g_{\rm{R}})\rangle=|\Psi(0)\rangle~{}.$ (3)

The derivative with respect to the coefficient $g_{\rm{R}}$ is equivalent to
the wavefunction $\Psi(0)$, which is a linear combination of quantum states.
This combination includes, in principle, any desired number of eigen-states
(with $\\{C_{N}\\}$ coefficients), which may include the ground-state
wavefunction of the system (with $S$ coefficient).

Similarly, we define the conjugate auxiliary wavefunction that depends on a
second parameter, $g_{\rm{L}}$:

$\langle\Psi_{\rm{L}}(0;g_{\rm{L}})|=\langle\Psi_{0}|+g_{\rm{L}}\langle\Psi(0)|~{},$
(4)

which follows $\lim_{g_{\rm{L}}\rightarrow 0}\partial/\partial
g_{\rm{L}}\langle\Psi_{\rm{L}}(0;g_{\rm{L}})|=\langle\Psi(0)|$. For CC
developments, we denote the operators $\hat{X}^{J}$ and $\hat{T}$ in terms of
the particle-hole excitation operators $\\{\hat{\tau}_{\mu}\\}$ (for an
excitation labeled $\mu$, when applied to $|0\rangle$, $\hat{\tau}_{\mu}$
promotes electrons from occupied levels into virtual ones and returns the
corresponding wavefunction determinant, [72]), so
$\hat{T}=\sum_{\mu}t_{\mu}\hat{\tau}_{\mu}$ and
$\hat{X}^{J}=\sum_{\mu}X_{\mu}^{J}\hat{\tau}_{\mu}$. $\hat{T}$ is the standard
ground-state cluster operator; $X^{J}_{\mu}$ is defined later below.

Even though the limit procedure is redundant for standard linear combinations,
for a CC wavefunction of the form:

$|\Phi(t=0;g_{\rm{R}})\rangle=\exp\Big{[}\hat{T}+g_{\rm{R}}\Big{(}S\hat{X}^{0}+\sum_{N}C_{N}\hat{X}^{N}\Big{)}\Big{]}|0\rangle~{},$
(5)

where $|0\rangle$ is the reference state, it leads to the CC analogue of
equation (3), which resembles $\Psi(0)$:

$\lim_{g_{\rm{R}}\rightarrow 0}\frac{\partial}{\partial
g_{\rm{R}}}|\Phi(t=0;g_{\rm{R}})\rangle=\Big{(}S\hat{X}^{0}+\sum_{N}C_{N}\hat{X}^{N}\Big{)}e^{\hat{T}}|0\rangle~{},$
(6)

where $\hat{X}^{N}$ is the CC excitation operator from equation-of-motion CC
(EOM-CC); $\hat{X}^{0}$ is simply the unity in this case [73]. The derivative
and limit then yield a linear combination of CC wavefunctions, in analogy with
the standard quantum mechanical approach. For instance, the analogue of the
eigen-function $|\Psi_{N}\rangle$ is the CC wavefunction
$\hat{X}^{N}e^{\hat{T}}|0\rangle$. Although these two wavefunctions can in
principle lead to the same excited state energies, they are not exactly the
same objects, due to the asymmetric (non-Hermitian) nature of standard CC
theory [32, 74]. For example, the CC analogue of $\langle\Psi_{N}|$ is
$\langle 0|\hat{\Lambda}^{N}\exp(-\hat{T})$, where $\hat{\Lambda}^{N}$ is the
conjugate excitation operator from EOM-CC, and it is expanded as
$\hat{\Lambda}^{N}=\sum_{\mu}\hat{\tau}^{\dagger}_{\mu}\Lambda^{N}_{\mu}$.
Therefore, the conjugate version is different (not simply a Hermitian
conjugate) because the “lambda” vectors must compensate for the term
$\exp[-\hat{T}-g_{\rm{R}}(S\hat{X}^{0}+\sum_{N}C_{N}\hat{X}^{N})]$. Hence we
introduce the wavefunction (the CC analogue of $\langle\Psi_{\rm{L}}|$):

$\langle\Upsilon(t=0;g_{\rm{L}},g_{\rm{R}})|=\langle
0|\hat{\lambda}_{\rm{m}}(0;g_{\rm{L}},g_{\rm{R}})\exp[-\hat{x}_{\rm{m}}(0;g_{\rm{R}})]~{},$
(7)

where
$\hat{x}_{\rm{m}}(0)=\hat{T}+g_{\rm{R}}(S\hat{X}^{0}+\sum_{N}C_{N}\hat{X}^{N})$,
and
$\hat{\lambda}_{\rm{m}}=\sum_{\mu}\hat{\tau}^{\dagger}_{\mu}\lambda_{\rm{m},\mu}$.
Our cluster operator expansions here include an identity term, labeled
$\mu=0$, where $\hat{\tau}_{0}=1$. Due to the non-Hermiticity, the left CC
wavefunction has information about the right-handed problem.

Returning to the standard quantum mechanical picture, we consider the
expectation value of the auxiliary TD observable $\tilde{A}(t)$:

$\tilde{A}(t)=\frac{\langle\Psi_{\rm{L}}(0)|\hat{A}^{\rm{H}}(t)|\Psi_{\rm{R}}(0)\rangle}{\mathcal{N}^{2}}~{},$
(8)

where
$\mathcal{N}^{2}=\langle\Psi_{\rm{L}}(0)|\Psi_{\rm{R}}(0)\rangle=1+g_{L}S^{*}+g_{R}S+g_{L}g_{R}$.
This function satisfies:

$\displaystyle\partial_{\rm{LR}}^{2}\tilde{A}(t)$
$\displaystyle=\lim_{g_{\rm{L}},g_{\rm{R}}\rightarrow
0}\frac{\partial^{2}}{\partial g_{\rm{L}}\partial g_{\rm{R}}}\tilde{A}(t)$
$\displaystyle=\langle\Psi_{\rm{E}}|\hat{A}^{\rm{H}}(t)|\Psi_{\rm{E}}\rangle-(1-|S|^{2})\langle\Psi_{0}|\hat{A}^{\rm{H}}(t)|\Psi_{0}\rangle~{},$

and

$\displaystyle\partial_{\rm{L}}\tilde{A}(t)$
$\displaystyle=\lim_{g_{\rm{L}},g_{\rm{R}}\rightarrow
0}\frac{\partial}{\partial
g_{\rm{L}}}\tilde{A}(t)=\langle\Psi_{0}|\hat{A}^{\rm{H}}(t)|\Psi_{\rm{E}}\rangle~{},$
(10) $\displaystyle\partial_{\rm{R}}\tilde{A}(t)$
$\displaystyle=\lim_{g_{\rm{L}},g_{\rm{R}}\rightarrow
0}\frac{\partial}{\partial
g_{\rm{R}}}\tilde{A}(t)=\langle\Psi_{\rm{E}}|\hat{A}^{\rm{H}}(t)|\Psi_{0}\rangle~{}.$

The last two equations define the operators
$\partial_{\rm{LR}}^{2}=\lim_{g_{\rm{L}},g_{\rm{R}}\rightarrow
0}\partial^{2}/\partial g_{\rm{L}}\partial g_{\rm{R}}$,
$\partial_{\rm{L}}=\lim_{g_{\rm{L}},g_{\rm{R}}\rightarrow 0}\partial/\partial
g_{\rm{L}}$, and $\partial_{\rm{R}}=\lim_{g_{\rm{L}},g_{\rm{R}}\rightarrow
0}\partial/\partial g_{\rm{R}}$. Using these definitions we obtain that the
left CC wavefunction must obey:

$\langle
0|\Big{(}S^{*}\hat{L}^{0}+\sum_{N}C_{N}^{*}\hat{\Lambda}^{N}\Big{)}e^{-\hat{T}}=\langle
0|\partial_{\rm{L}}\hat{\lambda}_{\rm{m}}(0)e^{-\hat{T}}~{},$ (11)

giving a term similar to $\langle\Psi(0)|$. We express
$\hat{L}_{0}=1+\hat{\Lambda}$, where $\hat{\Lambda}$ is the “lambda” operator,
which is the ground-state left analogue of $\hat{T}$.

From these two expressions we observe that a TD observable of a system that
initiates at a general quantum state, $\Psi(0)$, is obtained (which is a
desired result):

$\langle\Psi(t)|\hat{A}|\Psi(t)\rangle=\partial_{\rm{LR}}^{2}\tilde{A}+S\partial_{\rm{L}}\tilde{A}+S^{*}\partial_{\rm{R}}\tilde{A}+\langle\Psi_{0}|\hat{A}^{\rm{H}}(t)|\Psi_{0}\rangle~{}.$
(12)

In this work, motivated by this relation, we analyze the auxiliary observable
$\tilde{A}$ using cluster operators and use the results from such analysis to
represent $\langle\Psi(t)|\hat{A}|\Psi(t)\rangle$. In section 6 we apply such
expression to study a model 3-level/3-electron system.

###  2.2 Second response theory

The standard TD equations from CC theory are not limited to propagations from
the ground state. For example, using these standard CC equations, one can
imagine the situation where a system is driven from the ground-state to
another different quantum state. If the simulation were to be restarted at
such new quantum state, the system is indeed now being propagated from a state
that is not the ground state. However, the evolved “right” excitation TD
vector is not connected in a straightforward fashion to the EOM-CC eigen-
vectors, nor the conjugate TD excitation vector. This motivates our goal,
which is to describe the auxiliary observable shown in equation (8) through CC
quantities, then use of the differential step followed by the limiting
procedure presented in section 2.1 to compute
$\langle\Psi(t)|\hat{A}|\Psi(t)\rangle$ (similarly as in Ref. [68]), without
calculating the wavefunction $\Psi(t)$ of course. It allows for a connection
between EOM-CC eigen-vectors and general TD propagations, where one sets the
values of the initial-state coefficients $S$, and $\\{C_{N}\\}$. For this
reason, we examine differentiation with respect to the superposition
parameters $g_{\rm{L}}$ and $g_{\rm{R}}$ of the standard CC motion equations:

$\rmi\partial_{t}x_{\rm{m},\mu}(t)=\langle\hat{\tau}^{\dagger}_{\mu}e^{-\hat{x}_{\rm{m}}(t)}\hat{H}(t)e^{+\hat{x}_{\rm{m}}(t)}\rangle_{0}~{},$
(13)

and

$-\rmi\partial_{t}\lambda_{\rm{m},\mu}(t)=\langle\hat{\lambda}_{\rm{m}}(t)e^{-\hat{x}_{\rm{m}}(t)}[\hat{H}(t),\hat{\tau}_{\mu}]e^{+\hat{x}_{\rm{m}}(t)}\rangle_{0}~{},$
(14)

where the subscript “$\rm{m}$” indicates the vectors are modified objects as
they account for the new initial conditions (which depend on $g_{\rm{L}}$ and
$g_{\rm{R}}$), $\hat{x}_{\rm{m}}=\sum_{\mu}x_{\rm{m},\mu}\hat{\tau}_{\mu}$,
and $\langle\cdot\rangle_{0}$ denotes average with respect to the reference
($|0\rangle$): $\langle\cdot\rangle_{0}=\langle 0|\cdot|0\rangle$. The terms
$x_{\rm{m},0}$ and $\lambda_{\rm{m},0}$ are time-independent (this is a
consequence of our TD variational method, [68]). The phase follows the
relation:
$\partial_{t}\phi(t)=\langle\hat{\lambda}_{\rm{m}}(t)e^{-\hat{x}_{\rm{m}}(t)}\hat{H}(t)e^{+\hat{x}_{\rm{m}}}(t)\rangle_{0}$.
This number, $\phi(t)$, does not influence the calculation of observable
expectation values, so it will not be considered in detail in this work (its
properties are presented in Ref. [68]). These modified operators determine the
TD behavior of the following ket and its left conjugate:

$\displaystyle|\Phi(t;g_{\rm{R}})\rangle=\exp[\hat{x}_{\rm{m}}(t)-\rmi\phi(t)]|0\rangle~{},$
(15) $\displaystyle\langle\Upsilon(t;g_{\rm{R}},g_{\rm{L}})|=\langle
0|\hat{\lambda}_{\rm{m}}(t)\exp[-\hat{x}_{\rm{m}}(t)+\rmi\phi(t)]~{}.$ (16)

The inner product of these two CC wavefunctions is
$\langle\Upsilon(t)|\Phi(t)\rangle=\langle\hat{\lambda}_{\rm{m}}(t)\rangle_{0}$.

As shown in Ref. [68], there are four (reduced to three in this work) TD
excitation vectors that are needed to compute TD observables. These are:
$\hat{x}_{\rm{r}}(t)=\partial_{\rm{R}}\hat{x}_{\rm{m}}(t)$,
$\hat{\lambda}_{\rm{l}}(t)=\partial_{\rm{L}}\hat{\lambda}_{\rm{m}}(t)$,
$\hat{\lambda}_{\rm{r}}(t)=\partial_{\rm{R}}\hat{\lambda}_{\rm{m}}(t)$, and
$\hat{\lambda}_{\rm{l,r}}(t)=\partial_{\rm{LR}}^{2}\hat{\lambda}_{\rm{m}}(t)$.
In addition to such vectors, the excitation vectors for a propagation that
starts exclusively from the ground state are required as well in our
theoretical analysis (for this reason our formalism is called “second response
theory”). They are denoted as $\hat{x}(t)$, and $\hat{\lambda}(t)$, so
$\hat{x}(0)=\hat{T}$ and $\hat{\lambda}(0)=\hat{L}_{0}$, and they satisfy
equations (13) and (14), respectively. In the limit $g_{\rm{L}}$ and
$g_{\rm{R}}$ tend to zero, $\hat{x}_{\rm{m}}(t)$ and
$\hat{\lambda}_{\rm{m}}(t)$ tend to $\hat{x}(t)$ and $\hat{\lambda}(t)$,
respectively. Given the importance of $\lim_{g_{\rm{R}}\rightarrow
0}\hat{x}_{\rm{m}}(t)=\hat{x}(t)$, we define the quantity:

$\bar{O}_{x}(t)=e^{-\hat{x}(t)}\hat{O}(t)e^{+\hat{x}(t)}~{},$ (17)

where $\hat{O}(t)$ is some observable such as $\hat{H}(t)$ or $\hat{A}$. The
term $\bar{O}_{x,\tau,\mu}(t)$ is used to denote:

$\bar{O}_{x,\tau,\mu}(t)=[\bar{O}_{x}(t),\hat{\tau}_{\mu}]~{},$ (18)

and $\bar{O}_{T}=e^{-\hat{T}}\hat{O}e^{+\hat{T}}$. The $T$-transformed free
Hamiltonian is denoted as
$\bar{H}^{0}_{T}=e^{-\hat{T}}\hat{H}_{0}e^{+\hat{T}}$. The Jacobian matrix of
EOM-CC is defined as
$\mathcal{A}_{\mu\nu}=\langle\hat{\tau}^{\dagger}_{\mu}[\bar{H}_{T}^{0},\hat{\tau}_{\nu}]\rangle_{0}$.
In this way, the vectors $\mathbf{X}^{J}$ and $\bm{\Lambda}^{J}$ satisfy the
equation $\mathcal{A}\mathbf{X}^{J}=\Omega_{J}\mathbf{X}^{J}$ and
$(\bm{\Lambda}^{J})^{\rm{T}}\mathcal{A}=(\bm{\Lambda}^{J})^{\rm{T}}\Omega_{J}$,
where $(\mathbf{X}^{J})_{\mu}=X^{J}_{\mu}$,
$(\bm{\Lambda}^{J})_{\mu}=\Lambda^{J}_{\mu}$, and $\Omega_{J}$ is the
excitation energy (from ground to excited state, $\Omega_{J}=E_{J}-E_{0}$).

Using the operators defined above and based on section 2.1, the initial
conditions for the modified vectors are:

$\hat{x}_{\rm{m}}(0;g_{\rm{R}})=\hat{x}(0)+g_{\rm{R}}\hat{x}_{\rm{r}}(0)~{},$
(19)

and

$\hat{\lambda}_{\rm{m}}(0;g_{\rm{L}},g_{\rm{R}})=\hat{\lambda}(0)+g_{\rm{L}}\hat{\lambda}_{\rm{l}}(0)+g_{\rm{R}}\hat{\lambda}_{\rm{r}}(0)+g_{\rm{L}}g_{\rm{R}}\hat{\lambda}_{\rm{l,r}}(0)~{},$
(20)

where we will be interested in the limit where $g_{\rm{R}}$ and $g_{\rm{L}}$
tend to zero, as discussed later on. These limits are necessary to eliminate
higher-order terms in powers of $g_{\rm{R}}$ and $g_{\rm{L}}$ that the
modified vectors develop at $t>0$. As mentioned before, the asymmetric
dependencies on $g_{\rm{L}}$ and $g_{\rm{R}}$ by the excitation operators
$\hat{x}_{\rm{m}}(t)$ and $\hat{\lambda}_{\rm{m}}(t)$ are due to the non-
Hermitian attributes of the formalism (these would not take place in unitary
CC theory). To elaborate further on this, the vector $\hat{\lambda}_{\rm{m}}$
has different features than $\hat{x}_{\rm{m}}$ for two reasons: i), the left
expression has information about the right-handed problem; hence there must be
a vector $\hat{\lambda}_{\rm{r}}$. ii), The auxiliary quantum mechanical
observable $\tilde{A}(t)$ also depends on the product $g_{\rm{L}}g_{\rm{R}}$,
but the right CC wavefunction is independent of $g_{\rm{L}}$ and could not
introduce such dependency if $\hat{A}$ were the identity operator. Therefore
this dependency on $g_{\rm{L}}g_{\rm{R}}$ by the CC expression of
$\tilde{A}(t)$ requires the term $\hat{\lambda}_{\rm{l,r}}$, equation (20).

The operators $\hat{x}_{\rm{r}}$, $\hat{\lambda}_{\rm{l}}$,
$\hat{\lambda}_{\rm{r}}$, and $\hat{\lambda}_{\rm{l,r}}$ follow the relations
below [obtained from differentiating the motion equations of
$\hat{x}_{\rm{m}}$ and $\hat{\lambda}_{\rm{m}}$, equations (13) and (14)]
[68]:

$\displaystyle\rmi\partial_{t}x_{\rm{r},\mu}(t)=\langle\hat{\tau}^{\dagger}_{\mu}[\bar{H}_{x}(t),\hat{x}_{\rm{r}}(t)]\rangle_{0}~{},$
(21) $\displaystyle-$
$\displaystyle\rmi\partial_{t}\lambda_{\rm{l},\mu}(t)=\langle\hat{\lambda}_{\rm{l}}(t)[\bar{H}_{x}(t),\hat{\tau}_{\mu}]\rangle_{0}~{},$
(22) $\displaystyle-$
$\displaystyle\rmi\partial_{t}\lambda_{\rm{r},\mu}(t)=\langle\hat{\lambda}(t)[\bar{H}_{x,\tau,\mu}(t),\hat{x}_{\rm{r}}(t)]+\hat{\lambda}_{\rm{r}}(t)\bar{H}_{x,\tau,\mu}(t)\rangle_{0}~{},$
(23) $\displaystyle-$
$\displaystyle\rmi\partial_{t}\lambda_{\rm{l,r},\mu}(t)=\langle\hat{\lambda}_{\rm{l,r}}(t)\bar{H}_{x,\tau,\mu}(t)+\hat{\lambda}_{\rm{l}}(t)[\bar{H}_{x,\tau,\mu}(t),\hat{x}_{\rm{r}}(t)]\rangle_{0}~{}.$
(24)

With the equations developed above one can propagate the CC analogue of the
state $\Psi(0)=S\Psi_{0}+\sum_{N}C_{N}\Psi_{N}$ (or
$\Psi(0)=S\Psi_{0}+\Psi_{\rm{E}}$). This requires that the derivative vectors
be expressed in terms of other vectors that ensure solution to equations
(21-24) and that initially describe a term like
$\langle\Psi(0)|\hat{A}|\Psi(0)\rangle$. First, $\hat{x}_{\rm{r}}(t)$ follows
the initial condition:

$\hat{x}_{\rm{r}}(0)=S\hat{X}^{0}+\sum_{N}C_{N}\hat{X}^{N}~{},$ (25)

and is propagated as expressed in equation (21). Now, to further simplify the
solution of these equations, and generalize Ref. [68] to include a
contribution from the ground state, we write

$\hat{\lambda}_{\rm{l}}(t)=\hat{\lambda}^{\rm{E}}_{\rm{l}}(t)+S^{*}\hat{\lambda}(t)~{},$
(26)

where $\hat{\lambda}_{\rm{l}}^{\rm{E}}(0)=\sum_{N}C_{N}^{*}\hat{\Lambda}^{N}$,
$\hat{\lambda}(0)=\hat{L}_{0}$, and $t\geq 0$. For the right derivative,

$\hat{\lambda}_{\rm{r}}(t)=\hat{\lambda}_{\rm{r}}^{\rm{E}}(t)+S\hat{\lambda}(t)~{},$
(27)

where
$\hat{\lambda}_{\rm{r}}^{\rm{E}}(0)=-\sum_{N,I}C_{N}F^{NI}\hat{\Lambda}^{I}/(\Omega_{N}+\Omega_{I})$,
$t\geq 0$, and

$F^{NI}=\sum_{\mu\nu}\langle\hat{L}_{0}\big{[}[\bar{H}_{T}^{0},\hat{\tau}_{\mu}],\hat{\tau}_{\nu}\big{]}\rangle_{0}X^{N}_{\mu}X^{I}_{\nu}~{}.$
(28)

For the term $\hat{\lambda}_{\rm{l,r}}$ we have:

$\hat{\lambda}_{\rm{l,r}}(t)=\hat{\lambda}_{\rm{l,r}}^{\rm{E}}(t)+S\hat{\lambda}_{\rm{l}}^{\rm{E}}(t)+\hat{\lambda}(t)+S^{*}\hat{\lambda}_{\rm{r}}^{\rm{E}}(t)~{},$
(29)

where $\hat{\lambda}_{\rm{l,r}}^{\rm{E}}(t=0)=\sum_{J}Y^{J}\hat{\Lambda}^{J}$,
$t\geq 0$, and

$Y_{J}=\sum_{N,I}C_{N}^{*}C_{I}\frac{\langle\hat{\Lambda}^{N}\big{[}[\bar{H}^{0}_{T},\hat{X}^{J}],\hat{X}^{I}\rangle_{0}}{\Omega_{N}-\Omega_{J}-\Omega_{I}}~{}.$
(30)

This expression and that for $\hat{\lambda}_{\rm{l}}^{\rm{E}}(0)$ are needed
for the estimation of the excited-state component of the observable
$\tilde{A}(t)$ (as shown in section 3). By substituting equations (26-29) into
equations (22-24), we observe that
$\hat{\lambda}_{\rm{l}}^{\rm{E}},~{}\hat{\lambda}_{\rm{r}}^{\rm{E}}$ and
$\hat{\lambda}_{\rm{l,r}}^{\rm{E}}$ follow the same equations as
$\hat{\lambda}_{\rm{l}},~{}\hat{\lambda}_{\rm{r}}$ and
$\hat{\lambda}_{\rm{l,r}}$ (supplemental material). Furthermore, these
operators simplify the final expression for a TD observable.

We now show that the normalized auxiliary observable below provides excited-
state information:

$\tilde{A}_{\rm{N}}(t)=\frac{\langle\Upsilon(t)|\hat{A}|\Phi(t)\rangle}{\mathcal{N}_{2}^{2}}~{},$
(31)

where
$\mathcal{N}^{2}_{2}=\langle\Upsilon(t)|\Phi(t)\rangle=1+g_{\rm{L}}S^{*}+g_{\rm{R}}S+g_{L}g_{R}$.
The term $\langle\Upsilon(t)|\hat{A}|\Phi(t)\rangle$ is equivalent to
$\langle\hat{\lambda}_{\rm{m}}(t)\exp[-\hat{x}_{\rm{m}}(t)]\hat{A}\exp[+\hat{x}_{\rm{m}}(t)]\rangle_{0}$.
This norm is preserved because the relevant scalar amplitudes added to our CC
operators do not evolve as a function of time. This auxiliary quantity,
$\tilde{A}_{\rm{N}}$, is the CC analogue of the previously defined
$\tilde{A}(t)$ TD average.

We observe the following:

$\displaystyle\partial_{\rm{LR}}^{2}\tilde{A}_{\rm{N}}(t)$
$\displaystyle=\langle\hat{\lambda}_{\rm{l}}(t)[\bar{A}_{x}(t),\hat{x}_{\rm{r}}(t)]\rangle_{0}$
$\displaystyle~{}~{}~{}~{}-S^{*}\langle\\{\hat{\lambda}(t)[\bar{A}_{x}(t),\hat{x}_{\rm{r}}(t)]+\hat{\lambda}_{\rm{r}}(t)\bar{A}_{x}(t)\\}\rangle_{0}$
$\displaystyle~{}~{}~{}~{}+\langle\\{\hat{\lambda}_{\rm{l,r}}(t)-S\hat{\lambda}_{\rm{l}}(t)\\}\bar{A}_{x}(t)\rangle_{0}+(2|S|^{2}-1)\langle\hat{\lambda}(t)\bar{A}_{x}(t)\rangle_{0}~{},$
$\displaystyle=\langle\hat{\lambda}_{\rm{l}}^{\rm{E}}(t)[\bar{A}_{x}(t),\hat{x}_{\rm{r}}(t)]+\hat{\lambda}^{\rm{E}}_{\rm{l,r}}(t)\bar{A}_{x}(t)\rangle_{0}~{}.$

We then have that:

$\displaystyle\langle\Psi_{\rm{E}}|A^{\rm{H}}(t)|\Psi_{\rm{E}}\rangle=\langle\hat{\lambda}_{\rm{l}}^{\rm{E}}$
$\displaystyle(t)[\bar{A}_{x}(t),\hat{x}_{\rm{r}}(t)]+\hat{\lambda}^{\rm{E}}_{\rm{l,r}}(t)\bar{A}_{x}(t)\rangle_{0}+$
$\displaystyle(1-|S|^{2})\langle\hat{\lambda}(t)\bar{A}_{x}(t)\rangle_{0}~{}.$

The first-order derivatives follow:

$\displaystyle\partial_{\rm{L}}\tilde{A}_{\rm{N}}(t)$
$\displaystyle=\langle\hat{\lambda}_{\rm{l}}^{\rm{E}}(t)\bar{A}_{x}(t)\rangle_{0}~{},$
(34) $\displaystyle\partial_{\rm{R}}\tilde{A}_{\rm{N}}(t)$
$\displaystyle=\langle\hat{\lambda}_{\rm{r}}^{\rm{E}}(t)\bar{A}_{x}(t)+\hat{\lambda}(t)[\bar{A}_{x}(t),\hat{x}_{\rm{r}}(t)]\rangle_{0}~{},$
(35)

so
$\langle\Psi_{0}|\hat{A}^{\rm{H}}(t)|\Psi_{\rm{E}}\rangle=\partial_{\rm{R}}\tilde{A}_{\rm{N}}(t)$
and
$\langle\Psi_{\rm{E}}|\hat{A}^{\rm{H}}(t)|\Psi_{0}\rangle=\partial_{\rm{L}}\tilde{A}_{\rm{N}}(t)$.
Combining expressions and using equation (12) we obtain the formula:

$\langle\Psi(t)|\hat{A}|\Psi(t)\rangle=\langle\hat{\lambda}_{\rm{l}}(t)[\bar{A}_{x}(t),\hat{x}_{\rm{r}}(t)]+\hat{\lambda}_{\rm{l,r}}(t)\bar{A}_{x}(t)\rangle_{0}~{}.$
(36)

We also denote the right hand side of the above equation as
$\langle\hat{A}\rangle_{\rm{sr}}(t)$, with the “sr” subscript implying the
average is computed using the SR cluster operators. From this expression we
note that in order to predict the evolution of the observable one needs to
solve for three TD cluster operators. This equation above is also more compact
than that used in Ref. [68], which would require computing matrix elements
separately in the case where $S\neq 0$.

The relation shown above, equation (36), is not an exact equality, but a
rigorous assignment due to the non-Hermitian nature of the present CC
formalism. This aspect is also present in linear and quadratic CC response
theories [74]. So a deviation between standard quantum mechanics and the CC
formalism is expected to take place, but in our simulations this deviation is
small, as discussed later. This is an issue that can be solved by restoring
Hermiticity and unitarity, where our theoretical foundation remains
applicable.

###  2.3 Determining wavefunction probabilities and coherences

As discussed previously in Ref. [68], the CC wavefunctions do not provide
straightforwardly the probability of observing a given eigenstate at a
particular point time. In standard wavefunction theory this is accomplished by
projecting the eigenstate on the total wavefunction and taking square modulus
of the result, i.e., $p_{I}(t)=|\langle\Psi_{I}|\Psi(t)\rangle|^{2}$. As we
show in section 6, an alternative to this is to compute such amplitude as a
quantum mechanical observable, for example, by defining:

$\hat{\mathcal{P}}_{IJ}=e^{+\hat{T}}\hat{X}^{I}|0\rangle\langle
0|\hat{\Lambda}^{J}e^{-\hat{T}}~{},$ (37)

which is the CC analogue of $|\Psi_{I}\rangle\langle\Psi_{J}|$, setting
$\hat{A}=\hat{\mathcal{P}}_{II}$, and using equation (36) to compute
$p_{I}(t)$. More generally, we define
$\tilde{p}_{IJ}(t)=\langle\hat{\mathcal{P}}_{IJ}\rangle_{\rm{sr}}(t)$ [that
is, taking $\hat{A}$ as $\hat{\mathcal{P}}_{IJ}$ in equation (36) and
computing accordingly]. The wavefunction-based coherence matrix element
$p_{\rm{wf},IJ}(t)=C_{I}^{*}(t)C_{J}(t)$ has the evident property:
$p_{\rm{wf},IJ}^{*}(t)=p_{\rm{wf},JI}(t)$. To satisfy this we compute the
coherence as follows:

$\displaystyle{\rm Re}~{}p_{IJ}(t)$
$\displaystyle=\frac{1}{2}\Big{[}\tilde{p}_{IJ}(t)+\tilde{p}_{JI}(t)\Big{]}~{},$
(38) $\displaystyle{\rm Im}~{}p_{IJ}(t)$
$\displaystyle=\frac{1}{2}\Big{[}\tilde{p}_{IJ}(t)-\tilde{p}_{JI}(t)\Big{]}~{}.$

In section 6 we compare the complex number $p_{IJ}(t)$ to the object
$p_{\rm{wf},IJ}(t)$ defined above.

##  3 Unperturbed non-stationary states

In this section we show that well-known matrix elements from linear and
quadratic response CC theory can be obtained on the basis of unperturbed SR
non-stationary states. We present this development as an example of the
consistency of SR theory with respect to LR and QR theory results.

The initial conditions of the excitation vectors can be derived under two
considerations: i), The SR theory averages must predict very accurately the
exact unitary theory results, and, ii), if the system was originally described
by an excited state wavefunction, it should remain in that state when no
external perturbation is applied. This is known as a stationary state.
Therefore, we need to examine solutions to the set shown in equations (21-24)
is absence of perturbations, so $\bar{H}_{x}(t)=\bar{H}^{0}_{T}$. For
$\hat{x}_{\rm{r}}$ and $\hat{\lambda}_{\rm{l}}$, we have that:
$\hat{x}_{\rm{r}}(t)=\hat{X}^{N}\exp(-\rmi\Omega_{N}t)$ and
$\hat{\lambda}_{\rm{l}}(t)=\hat{\Lambda}^{N}\exp(+\rmi\Omega_{N}t)$. These
ensure solution to the Jacobian eigenvalue problem:
$\rmi\partial_{t}\mathbf{x}_{\rm{r}}=\Omega_{N}\mathbf{x}_{\rm{r}}=\mathcal{A}\mathbf{x}_{\rm{r}}$,
and
$-\rmi\partial_{t}\bm{\lambda}_{\rm{l}}^{\rm{T}}=\Omega_{N}\bm{\lambda}_{\rm{l}}^{\rm{T}}=\bm{\lambda}^{\rm{T}}_{\rm{l}}\mathcal{A}$
(where “$\rm{T}$” denotes transpose operation).

Now suppose the system is described by a linear superposition of states,
leading to an unperturbed non-stationary state, so

$\displaystyle\hat{x}_{\rm{r}}(t)$
$\displaystyle=S\hat{X}^{0}+\sum_{N}C_{N}\hat{X}^{N}\exp(-\rmi\Omega_{N}t)~{},$
(39) $\displaystyle\hat{\lambda}_{\rm{l}}(t)$
$\displaystyle=S^{*}\hat{L}_{0}+\sum_{N}C_{N}^{*}\hat{\Lambda}^{N}\exp(+\rmi\Omega_{N}t)~{}.$
(40)

With the two solutions above, we can find $\hat{\lambda}_{\rm{r}}$ and
$\hat{\lambda}_{\rm{l,r}}$. First, we express
$\hat{\lambda}_{\rm{r}}(t)=\sum_{J}d_{J}(t)\hat{\Lambda}^{J}+S\hat{L}_{0}$,
where $\\{d_{J}\\}$ are TD coefficients (from section 2.2 we note that
$\hat{\lambda}(t)=\hat{L}_{0}$ and $\hat{x}(t)=\hat{T}$ in absence of external
perturbations). By using the bi-orthogonality property of the excitation
vectors $\\{\hat{X}^{J},\hat{\Lambda}^{J}\\}$ we observe that

$-\rmi\partial_{t}d_{J}=\sum_{N}C_{N}F^{JN}e^{-\rmi\Omega_{N}t}+\Omega_{J}d_{J}(t)~{},$
(41)

so
$d_{J}(t)=-\sum_{N}C_{N}F^{NJ}\exp(-\rmi\Omega_{N}t)/(\Omega_{N}+\Omega_{J})$.
This gives:

$\hat{\lambda}_{\rm{r}}^{\rm{E}}(t)=-\sum_{N,J}\frac{F^{NJ}}{\Omega_{N}+\Omega_{J}}\hat{\Lambda}^{J}e^{-\rmi\Omega_{N}t}~{},$
(42)

which in turn leads to the TD dependency of $\hat{\lambda}_{\rm{r}}$, that is,
$\hat{\lambda}_{\rm{r}}(t)=\hat{\lambda}^{\rm{E}}_{\rm{r}}(t)+S\hat{L}_{0}$.
In a similar fashion, if we write
$\hat{\lambda}_{\rm{l,r}}(t)=\sum_{J}y_{J}(t)\hat{\Lambda}^{J}+S\hat{\lambda}_{\rm{l}}(t)+\hat{\lambda}(t)+S^{*}\hat{\lambda}^{\rm{E}}_{\rm{r}}(t)$,
we have that

$\displaystyle-\rmi\partial_{t}$ $\displaystyle y_{J}=$
$\displaystyle\Omega_{J}y_{J}(t)+\sum_{M,N}C^{*}_{M}C_{N}\exp[-\rmi(\Omega_{N}-\Omega_{M})t]\langle\hat{\Lambda}^{M}\big{[}[\bar{H}^{0}_{T},\hat{X}^{J}],\hat{X}^{N}\big{]}\rangle_{0}~{},$

where the sums run over excited states. The solution to the above equation is:

$y_{J}(t)=\sum_{M,N}C^{*}_{M}C_{N}\exp[-\rmi(\Omega_{N}-\Omega_{M})t]\frac{\langle\hat{\Lambda}^{M}\big{[}[\bar{H}^{0}_{T},\hat{X}^{J}],\hat{X}^{N}\big{]}\rangle_{0}}{\Omega_{M}-\Omega_{J}-\Omega_{N}}~{}.$
(44)

Using the last relation we have that $Y_{J}=y_{J}(t=0)$.

The expressions above for perturbation-free evolutions lead to results
consistent with linear and quadratic CC response theories. Now, suppose the TD
wavefunction is given by $\Psi(t)=\exp(-\rmi
E_{0}t)[S\Psi_{0}+\sum_{N}C_{N}\Psi_{N}\exp(-\rmi\Omega_{N}t)]$. We note that:

$\displaystyle\langle\Psi(t)|\hat{A}|\Psi(t)\rangle=|S|^{2}$
$\displaystyle\langle\Psi_{0}|\hat{A}|\Psi_{0}\rangle+\Big{[}S^{*}\sum_{N}e^{-\rmi\Omega_{N}t}\langle\Psi_{0}|\hat{A}|\Psi_{N}\rangle+\rm{c.c.}\Big{]}$
(45)
$\displaystyle+\sum_{M,N}C_{M}^{*}C_{N}e^{-\rmi(\Omega_{N}-\Omega_{M})t}\langle\Psi_{M}|\hat{A}|\Psi_{N}\rangle~{}.$
(46)

Examining the CC analogue of the above equation, based on common phasor terms,
we identify the required matrix elements. These are, the term

$\langle\Psi_{N}|\hat{A}|\Psi_{0}\rangle=\langle\hat{\Lambda}^{N}\bar{A}_{T}\rangle_{0}~{},$
(47)

its conjugate:

$\langle\Psi_{0}|\hat{A}|\Psi_{N}\rangle=\langle\hat{L}_{0}[\bar{A}_{T},\hat{X}^{N}]\rangle_{0}-\sum_{J}\frac{F^{NJ}}{\Omega_{J}+\Omega_{N}}\langle\hat{\Lambda}^{J}\bar{A}_{T}\rangle_{0}~{},$
(48)

and the excited-state/excited-state element:

$\displaystyle\langle\Psi_{M}|\hat{A}|\Psi_{N}\rangle=\delta_{MN}$
$\displaystyle\langle\hat{L}_{0}\bar{A}_{T}\rangle_{0}+\langle\hat{\Lambda}^{M}[\bar{A}_{T},\hat{X}^{N}]\rangle_{0}$
(49)
$\displaystyle+\sum_{J}\frac{\langle\hat{\Lambda}^{M}\big{[}[\bar{H}^{0}_{T},\hat{X}^{J}],\hat{X}^{N}\big{]}\rangle_{0}}{\Omega_{M}-\Omega_{J}-\Omega_{N}}\langle\hat{\Lambda}^{J}\bar{A}_{T}\rangle_{0}~{}.$
(50)

Equations (47) and (48) are identical to those from CC linear response theory,
and equation (50) was derived originally within the context of CC quadratic
response theory [32]. These equations clearly illustrate the non-Hermiticity
of the present formalism, where interchanging labels ($N\leftrightarrow M$,
for example) gives different formulas [74].

##  4 Approximated scheme based on regularization

Even though standard CC theory is designed to address dynamically correlated
systems, or more challenging cases through spin symmetry breaking, there can
be scenarios where solving the CC equations is difficult due to low-energy
gaps or related instabilities in the numerical solution methods used to
determine $\\{t_{\mu}\\}$ (these can occur in perturbation theory applications
[75, 76]). The CC variational method is not necessarily restricted in use to
non-relativistic Coulombic Hamiltonians, but it can also be applied to double-
hybrid DFT Hamiltonians, and model systems such as those described by Hubbard
Hamiltonians (or similar), where this issue may occur. To circumvent numerical
instabilities (i.e., divergences in $\\{t_{\mu}\\}$) when they take place and
facilitate the determination of a wavefunction of the form
$\exp(\hat{T})|0\rangle$, regularization can be an appealing approach. It may
not provide a completely physically meaningful wavefunction in some cases, but
such wavefunction could be used as a starting point for alternative formalisms
that could need it. In the context of DFT, we showed for the dissociation of
fundamental diatomic molecules that regularization can eliminate instabilities
and give results with consistent physical meaning [77].

As we proposed it previously [78], one can regularize the equations by adding
an additional term, and find the stationary point of the function (inspired
from the standard variational method [79, 80]):

$\mathcal{E}(\alpha)=\binrel@{\rm{stat.}}\binrel@@{\mathop{\kern
0.0pt\rm{stat.}}\limits_{\bm{\Lambda},\mathbf{t}}}\Big{[}\langle(1+\hat{\Lambda})e^{-\hat{T}}\hat{H}_{0}e^{+\hat{T}}\rangle_{0}+\alpha\bm{\Lambda}\cdot\mathbf{t}\Big{]}~{},$
(51)

where $\alpha$ is a regularization number,
$\bm{\Lambda}\cdot\mathbf{t}=\sum_{\mu}\Lambda_{\mu}t_{\mu}$, and
“$\rm{stat.}$” stands for stationarization (the minimization principle does
not apply to this case). We denote the solution to the problem above as
$\\{\Lambda_{\mu}^{\prime}\\}$ and $\\{t_{\mu}^{\prime}\\}$ (the prime symbol
indicating these are regularized solutions). We also define
$\hat{L}^{\prime}=1+\hat{\Lambda}^{\prime}$. In this section we focus on the
eigen-value problem.

We wish to apply the SR method in this case for $t=0$ only. With the solution
regularized, the wavefunction $\Psi_{\rm{R}}$ takes a different form, where
the term that is not multiplied by $g_{\rm{R}}$ is a reference state that is
the standard quantum mechanical analogue of $\exp(\hat{T}^{\prime})|0\rangle$
(which is not the true ground state CC wavefunction). If we denote this
wavefunction as $\Psi_{0}^{\prime}$, then
$|\Psi_{R}(0,g_{R})\rangle=|\Psi_{0}^{\prime}\rangle+g_{\rm{R}}|\Psi(0)\rangle,$
where $\Psi(0)$ is a state we want to represent using a CC object, such as the
true ground state or an excited state of interest. The CC wavefunctions that
we seek here are of the form $\hat{X}^{N}\exp(\hat{T}^{\prime})|0\rangle$ and
$\langle 0|\exp(-\hat{T}^{\prime})\hat{\Lambda}^{N}$. So we take
$\hat{x}(0)=\hat{T}^{\prime}$ and $\lambda(0)=\hat{L}^{\prime}$.

For $\hat{A}=\hat{H}_{0}$, we then have, from equation (36), the
approximation:

$E=\langle\hat{\Lambda}[\bar{H}^{0}_{T^{\prime}},\hat{X}]\rangle_{0}+\langle\hat{L}^{\prime}\bar{H}^{0}_{T^{\prime}}\rangle_{0}~{},$
(52)

where $E$ represents $\langle\Psi(0)|\hat{H}|\Psi(0)\rangle$, and
$\bar{H}^{0}_{T^{\prime}}=e^{-\hat{T}^{\prime}}\hat{H}_{0}e^{+\hat{T}^{\prime}}$.
We treat $\hat{X}$ and $\hat{\Lambda}$ as variational quantities. We also
define
$\tilde{\mathcal{A}}_{\mu\nu}=\langle\hat{\tau}^{\dagger}_{\mu}[\bar{H}_{T^{\prime}},\hat{\tau}_{\nu}]\rangle_{0},$
which is now based on the regularized Hamiltonian. Note that the identity
operator $\hat{\tau}_{0}=1$ is not required to expand the excitation operators
($\hat{X}$ and $\hat{\Lambda}$) because it commutes with
$\bar{H}_{T^{\prime}}$ (but in a different formulation including
$\hat{\tau}_{0}$ may be necessary).

With these definitions we have:

$E=(\bm{\Lambda})^{\rm{T}}\tilde{\mathcal{A}}\mathbf{X}+\langle\hat{L}^{\prime}\bar{H}^{0}_{T^{\prime}}\rangle_{0}~{},$
(53)

and the eigenvalue problem takes the form:
$(\bm{\Lambda}^{N})^{\rm{T}}\tilde{\mathcal{A}}=(\bm{\Lambda}^{N})^{\rm{T}}\tilde{\Omega}_{N}$
and $\tilde{\mathcal{A}}\mathbf{X}^{N}=\tilde{\Omega}_{N}\mathbf{X}^{N}$,
where $\tilde{\Omega}_{I}=E_{I}-E^{\prime}_{0}$ and
$E^{\prime}_{0}=\langle\hat{L}^{\prime}\bar{H}^{0}_{T^{\prime}}\rangle_{0}$.
The estimated energies of interest are then
$E_{I}=\tilde{\Omega}_{I}+E^{\prime}_{0}$. As expected, by solving this
problem one also obtains approximated values for the excited-state energies
($\\{E_{I}\\},~{}I\neq 0$).

Assuming that $\hat{\lambda}(0)=\hat{L}^{\prime}$ renders the present method
as approximated in essence. It can be shown that in a more rigorous approach
$\hat{\lambda}(0)$ involves more terms, but this requires a careful response
theory analysis that includes the fact that $\exp(\hat{T}^{\prime})|0\rangle$
is a linear combination of the eigen-states of the Hamiltonian $\hat{H}_{0}$,
which is beyond the current scope.

##  5 Computational methodology

Figure 1: Electronic configurations used (left), and operators needed to
compute properties arising from propagations that start from a quantum linear
superposition (right). A quantum mechanical TD average is then computed using
the vectors shown in this figure, and equation (36). The reference is the
state $|0\rangle$.

To explore the SR CC theoretical propagation method developed in this work, we
examine a three-level/three-electron system. In the reference configuration
$|0\rangle$, the first energy level, labeled “$j$”, is doubly occupied, the
next level, “$i$”, is singly occupied, and the virtual level, “$a$”, is
unoccupied. We consider a non-trivial Hamiltonian that features two-electron
interactions, and a strong driving perturbation. The unperturbed Hamiltonian
has the form:

$\hat{H}_{0}=\hat{h}_{0}+\sum_{q}U_{q}\hat{n}_{q\uparrow}\hat{n}_{q\downarrow}+\Big{[}\sum_{\mu>\nu}W_{\mu\nu}\hat{\tau}_{\mu}\hat{\tau}_{\nu}+\rm{H.c.}\Big{]}~{},$
(54)

where the single particle Hamiltonian, $\hat{h}_{0}$, reads:
$\hat{h}_{0}=\sum_{p}\epsilon_{p}\hat{n}_{p}+\sum_{\mu}b_{\mu}[\hat{\tau}_{\mu}+\hat{\tau}^{\dagger}_{\mu}]$.
The number operator for a spin-level reads
$\hat{n}_{q\sigma}=\hat{q}^{\dagger}_{\sigma}\hat{q}_{\sigma}$, and
$\hat{n}_{p}=\hat{n}_{p\downarrow}+\hat{n}_{p\uparrow}$ is the number operator
of level $p$. We take $b_{\mu}=b$, $W_{\mu\nu}=w_{0}$, $U_{a}=U_{i}=U_{0}$,
and $U_{j}=0$. All the single-particle orbital levels are equally spaced by a
$\Delta\epsilon$ term, and we also set $\epsilon_{j}=0$ for convenience.

Having defined the free Hamiltonian of the system, the full TD energy operator
is defined as:

$\hat{H}(t)=\hat{H}_{0}-f(t)\hat{D}~{}.$ (55)

In this case we take the dipole operator as
$\hat{D}=d_{0}\sum_{\mu}[\hat{\tau}_{\mu}+\hat{\tau}^{\dagger}_{\mu}]$, and
the driving function corresponds to a rectangular pulse of the form
$f(t)=f_{0}$ if $0<t<t_{0}$ ($t_{0}=5~{}\rm{fs}$), $f(t)=0$, otherwise.

In section 6 we propagate standard and CC-based quantum superpositions and
compare the results. We then show that the regularized scheme leads to
consistent results with respect to both the reference and the SR CC results.
For these numerical simulations we assume that $\Delta\epsilon=1~{}\rm{eV}$,
$\epsilon_{j}=0$, $\epsilon_{i}=\Delta\epsilon$, and
$\epsilon_{a}=2\Delta\epsilon$. We take $b=0.1~{}\rm{eV}$,
$w_{0}=U_{0}=0.2~{}\rm{eV}$, $d_{0}=0.5$, and $f_{0}=0.04~{}\rm{au}$ (or
$2.1\times 10^{10}~{}{\rm V/m}$); other values are specified accordingly in
section 6. Figure 1 shows a basic schematics of the three-electron/three-level
system of interest in terms of its single-particle basis. There are nine
configurations that lead to nine full-body wavefunctions.

In physical terms, an $M_{\rm{s}}=+1/2$ system with equally “spaced” energy
levels, initially in a linear combination of states, is subject to an external
perturbation. We use our SR CC method to compute expectation values and
compare against numerically exact standard (unitary) quantum mechanics. In
figure 1 we show the steps required to perform a propagation based on SR
theory. First, one solves the ground- and excited-state CC problems, which
give the operators $\hat{T}$, $\hat{\Lambda}$, $\\{\hat{X}^{N}\\}$, and
$\\{\hat{\Lambda}^{N}\\}$. Then one determines the operator $\hat{x}(t)$,
through solution of the standard TD CC equations that propagate a system that
initiates at the ground-state. Following this step, one solves the SR TD
relations, equations (21, 22, and 24), for $\hat{x}_{\rm{r}}$,
$\hat{\lambda}_{\rm{l}}$, and $\hat{\lambda}_{\rm{l,r}}$, and uses these
vectors in equation (36) to compute the TD observable $A$. The initial
conditions are applied as well, as discussed in section 2.2, and which are
based on the $S$, and $\\{C_{N}\\}$ coefficients. To solve the TD SR
equations, we use the midpoint method (or second order Runge-Kutta technique)
in the range $[0-50~{}\rm{fs}]$; this is based on 50,000 steps.

##  6 Results and discussion

Figure 2: Standard propagation and driving pulse. a), Time-dependent dipole,
$D_{0}(t)=\langle\Psi_{0}|\hat{D}^{\rm{H}}(t)|\Psi_{0}\rangle$, for
propagation from the ground-state of the system: open circles denote the
standard TD CC propagation, purple line refers to the unitary reference result
(numerically exact). b), Shape of the pulse applied (shown between 0 and 10
fs), which is turned-off at $t=5~{}\rm{fs}$.

The standard eigenstates of the system are obtained by diagonalizing the
unperturbed Hamiltonian, such that
$\hat{H}_{0}|\Psi_{J}\rangle=E_{J}|\Psi_{J}\rangle$, $J=0,1,\ldots,8$ (details
about these states can be found in the supplemental material). In an analogous
fashion, we first solve the ground-state CC problem iteratively to find
$\hat{T}$ and $\hat{\Lambda}$, then the left and right eigenvalue problems
(with matrix $\mathcal{A}$) are solved to determine the set
$\\{\mathbf{\Lambda}^{I},~{}\mathbf{X}^{I}\\}$ ($I=1,\ldots,8$). Having these
eigenvalue problems solved, in order to generate reference quantum mechanical
data, we proceed to propagate the standard wavefunction of the system using
$\rmi\partial_{t}\Psi(t)=\hat{H}(t)\Psi(t)$ [as $\Psi(t+\delta
t)=\exp\\{-\rmi\delta t\hat{H}(t+\delta t/2)\\}\Psi(t)$]. Then, we solve the
second response equations for $\hat{\lambda}_{\rm{l}}$,
$\hat{\lambda}_{\rm{l,r}}$, and $\hat{x}_{\rm{r}}$, equations (21, 22, and
24), as mentioned in section 5. The TD observable is computed for both cases
using $\langle\Psi(0)|\hat{A}^{\rm{H}}(t)|\Psi(0)\rangle$ (which serves as
reference) and
$\langle\hat{\lambda}_{\rm{l}}(t)[\bar{A}_{x}(t),\hat{x}_{\rm{r}}(t)]+\hat{\lambda}_{\rm{l,r}}(t)\bar{A}_{x}(t)\rangle_{0}$.
We study the propagation of three different initial states, denoted “QS1”,
“QS2”, and “QS3”, respectively (“QS” standing for “quantum superposition”).
These initial states are defined as follows:

$\displaystyle\rm{QS1}$
$\displaystyle:~{}~{}\Psi(0)=\frac{1}{\sqrt{3}}\Psi_{0}+\frac{1}{\sqrt{3}}\Psi_{1}+\frac{1}{\sqrt{3}}\Psi_{2}~{},$
(56) $\displaystyle\rm{QS2}$
$\displaystyle:~{}~{}\Psi(0)=\frac{1}{\sqrt{2}}\Psi_{7}+\frac{\rmi}{\sqrt{2}}\Psi_{8}~{},$
(57) $\displaystyle\rm{QS3}$
$\displaystyle:~{}~{}\Psi(0)=\frac{1}{\sqrt{4}}\Psi_{0}+\frac{1}{\sqrt{4}}\Psi_{3}+\frac{\rm{1}}{\sqrt{2}}\Psi_{5}~{}.$
(58)

Note that QS2 features a complex-valued coefficient ($C_{8}=i/\sqrt{2}$).

Figure 3: Evolution of the system as initiating in the QS1 state in terms of
the dipole and probabilities, $D(t)=\langle{\rm QS}1|\hat{D}^{\rm{H}}(t)|{\rm
QS}1\rangle$, where $|\rm{QS1}\rangle$ denotes the first initial quantum
superposition shown in equation (56). a), TD dipole: purple line corresponds
to the numerically exact unitary result, open circles denote the SR CC theory
results. b), Probabilities of the unperturbed Hamiltonian eigenstates
$\Psi_{1}$ (purple line: exact result, cross symbol: SR theory), $\Psi_{2}$
(blue line: exact result, “+” symbol: SR theory), and $\Psi_{6}$ (red line:
unitary simulation, open circle: SR theory). Figure 4: a), TD population of
level $a$, purple line: exact result, cross symbols: SR theory. b), TD
population of level $i$, red line: unitary simulation, open circles: SR CC
theory.

The Hamiltonian $\hat{H}_{0}$ yields a ground state that is non-trivial, where
the contributions of the singles and doubles configurations are in the range
$[0.04-0.1]$. These are significant contributions to the reference of the
system. Figure 2.a shows the result of propagating the ground-state using the
standard unitary and the CC methods. As expected, the TD standard CC equations
reproduce quite well the unitary results. As mentioned before, from this
propagation one needs the vector $\hat{x}(t)$ to propagate quantum
superpositions; the operator $\hat{\lambda}(t)$ is not needed as
$\hat{\lambda}_{\rm{l,r}}(t)$ contains the information about this cluster
operator. To illustrate the applicability of the SR method, we begin
considering the initial state “QS1”,
$\Psi(0)=C_{0}\Psi_{0}+C_{1}\Psi_{1}+C_{2}\Psi_{2}$, where
$C_{0}=C_{1}=C_{2}=1/\sqrt{3}$ (also note that $S=C_{0}$). The state
$\Psi_{0}$ is dominated by the reference configuration ($|0\rangle$),
$\Psi_{1}$ by the first configuration ($|1\rangle$), and $\Psi_{2}$ by the
third configuration ($|3\rangle$); this also applies to the their
corresponding EOM-CC eigen-states. Figure 3.a displays the comparison between
the standard unitary and SR theory propagations for the dipole, and figure 3.b
the eigenstate probabilities. The agreement between these two propagations is
noticeable. Our propagations are entirely coherent due to the absence of
dissipative effects. For this reason, after the perturbation is turned off,
the dipole remains oscillating. The exact probabilities, as expected, reach a
steady value after $t=5~{}\rm{fs}$. The estimates obtained from the CC
operator $\hat{\mathcal{P}}_{II}$, except $p_{6}(t)$, however, show relatively
small deviations after $5~{}\rm{fs}$. The use of this operator, nonetheless,
yields a more desirable consistency than the other estimators discussed
previously in Ref. [68]. We attribute this again to the asymmetric
characteristics of the formalism and the truncation (up to the fourth order)
of $\exp[-\hat{x}(t)]\hat{\mathcal{P}}_{IJ}\exp[\hat{x}(t)]$. Figure 4
displays the TD evolution of the populations of the levels $i$ and $a$, where
a better overall agreement is observed in comparison to the eigenstate
probability case.

Figure 5: Propagation initiating at the QS2 state. a), TD dipole of the
system. Purple line denotes the unitary reference simulation, open circles
refer to the SR theory results. b), TD populations of eigenstates $\Psi_{7}$
(purple line: exact result, cross symbol: SR simulation) and $\Psi_{8}$ (blue
line: unitary reference result, “+” symbol: SR theory). c), Coherence of
states 7 and 8 (or $p_{78}(t)$) as a function of time; red line: exact result
for the real part of this coherence, open triangle: SR theory for the real
part; purple line: exact result for imaginary component, open circle: SR
theory for this imaginary part.

A similar consistent behavior between SR CC and unitary results holds for the
other two quantum superpositions, QS2 and QS3. For propagation starting at
QS2, in figure 5 we can note the excellent agreement for the TD dipole and
probabilities. QS2 is a 50 %-50 % combination of the last two excited states,
$\Psi_{7}$ and $\Psi_{8}$. $\Psi_{7}$ is dominated by the 5th and 7th
configurations, whereas $\Psi_{8}$ is dominated by the 8th configuration;
these are doubles determinants. In figure 5.a we observe that the oscillating
patterns of the dipole change after the external field is turned off. Even
though the field is on for a short period time, this range of time
($0-5~{}\rm{fs}$) is sufficient to drive the eigenstate TD probabilities in a
non-trivial fashion, as evidenced in figure 5, where a strong population
change occurs. Simulations where the external field is deactivated at longer
times were also explored in our present studies without noticing differences
in deviations with respect to those presented in this section. On the other
hand, the coherence term $p_{\rm{wf},78}=C_{7}^{*}(t)C_{8}(t)$ is also
described with a striking agreement by the CC operator
$\hat{\mathcal{P}}_{78}$ in conjunction with equation (38), figure 5.c. These
results suggest the SR formalism offers a promising theoretical method to
propagate quantum superpositions while preserving connectedness of quantum
mechanical expressions, a condition required to approach macroscopic/periodic
systems; where our developments, besides expanding Ref. [68] with the compact
inclusion of general quantum superpositions, also provide tools to explore
non-diagonal elements of the quantum mechanical density matrix, which depends
on an operator such as $\hat{\mathcal{P}}_{IJ}$.

In our simulations, the operator $\hat{\mathcal{P}}_{IJ}$ was expressed in
matrix form by explicitly computing the exponential matrices involved. The
resulting matrix was treated as a matrix associated to an observable, which
was then symmetrized. In more practical settings, however, computationally
efficient ways to use the operator $\hat{\mathcal{P}}_{IJ}$ would be needed.
This can be carried out by examining the expressions that result from
inserting this operator in equation 12, where a considerable aspect to address
in future work concerns truncation of the exponential operators. Furthermore,
the symmetrization step is not unique, as expected, as other procedures are
possible, and the trace is not strictly the unity if the exponential operator
is computed as a truncated power series; otherwise the approach is norm-
conserving. In a unitary version of SR CC, however, this issue may not arise.

Figure 6: Effect of regularization on dipole and population evolution,
starting from QS3. a), TD dipole where the purple and blue lines refer to the
numerically exact and the regularized SR theory results, respectively. b) and
c) show the comparison between the regularized (black line) and unitary
(purple line) reference for levels $a$ and $i$, correspondingly.

Our SR CC theory results provide support for equation (36), which describes
the evolution of a quantum mechanical observable. To further test equation
(36), we now examine the response of states that derive from the regularized
ground state and the diagonalization of $\tilde{\mathcal{A}}$. Our
regularization technique approximates energies and CC vectors, and thereby the
associated TD observables. However, it can provide useful results to analyze
and understand quantum mechanical quantities. We examine here the case where
$\alpha$ is relatively large (compared to $\Delta\epsilon$),
$\alpha=4.0~{}\rm{eV}$. For the iterative calculation of the $t$-amplitudes,
in the quasi-Newton step, which involves energy differences in a series of
denominators, the regularization number broadens such energy differences by
its value, for example, $\epsilon_{a}-\epsilon_{i}+\alpha$, and any other
energy difference between the configurations. So if the gap between virtual
and occupied levels is relatively small, the regularization removes potential
instabilities and yields convergent, approximate $t,~{}\Lambda$ amplitudes.
Tables S1 and S2 (supplemental material) show the regularized and exact energy
values for the system considered in figure 1. We note that with greater
regularization values, the errors increase slightly. The errors, have an upper
bound because the larger the $\alpha$ number is, the more negligible the
resulting $\hat{T}$ and $\hat{\Lambda}$ vectors become. Nonetheless, there are
elements in our theory to improve (not explored in this work) the accuracy of
a regularized method. Regularization, in our opinion, is convenient for cases
where instabilities arise in the determination of the ground-state cluster
amplitudes. To study the response of the system after the ground-state
regularized amplitudes are found, we follow the same steps as for the
unregularized propagations discussed before. Figure 6 shows the evolution of
the observable $\hat{D}$ for the regularized system with
$\alpha=4.0~{}\rm{eV}$. The system in this case begins from the state QS3. As
opposed to the unregularized calculations, small deviations are slightly more
noticeable in this case due to the effect of the regularization. This can be
observed for the dipole, figure 6.a, and the populations of the $a$ and $i$
levels, figures 6.b and .c, respectively. Such deviations are caused by
fluctuations in the excitation energies with respect to the exact ones, and in
the EOM-CC wavefunctions.

As mentioned in Ref. [68], the present SR CC approach is asymmetric and, in
its current form, is unable to reproduce exactly the standard wavefunction
theory results, despite these agree strikingly well. Nonetheless, the
robustness of the non-Hermitian SR CC model holds even if the strength of the
electron-electron interaction is increased, as shown in figure S1
(supplemental material). Repeating the steps outlined in section 5, but now
for $b=0.25~{}\rm{eV}$ and $w_{0}=U_{0}=0.5~{}\rm{eV}$, we observe the same
kind of overall agreement between the SR theoretical propagation and the
unitary reference evolution. The ground-state amplitudes are in this scenario
in the range $[0.08-0.2]$, reflecting a strong electron-electron interaction
(Tables S3 and S4, supplemental material). Similarly, the regularization
method yields the TD quantum dipole in close agreement with both the
unregularized SR and standard quantum mechanical results. Deviations with
respect to standard simulations occur in our theory because the SR CC
wavefunctions display a type of normalization that does not take place in the
parent quantum mechanical theory. Even though one can normalize wavefunctions
in a form that complies with the Hermitian theory, this would require
numerical operations that work against the practical purpose of asymmetric CC
theories. However, the unitary CC method is an alternative to the non-
Hermitian CC theory that is widely applied in quantum algorithms for
simulations in chemistry. The SR theory is applicable to expand TD unitary CC
theory, but it necessitates a separate careful analysis not covered in this
work. Nonetheless, the present SR theory captures accurately the TD behavior
of the observable represented by $\hat{A}$.

Returning to potential applications of the SR CC theory in realistic cases,
for strong external perturbations that are turned on for long enough such that
they induce a significant response by the quantum system, high-order
excitation vectors of fourth order or beyond could be necessary. In these
situations, an alternative approach to avoid such high-order excitation
vectors is the use of TD relaxed orbitals [81, 82, 49]. For example, a
possibility to expand the SR formalism in this direction is to employ
auxiliary reference TD orbitals that derive from a TD variational method. Such
orbitals could also be treated as functions of the parameters $g_{\rm{L}}$ and
$g_{\rm{R}}$. These would in turn provide a basis to find a different set of
TD equations for the excitation vectors that may be truncated at the third or
second order.

Our proposed SR CC theory derives from the standard (non-Hermitian yet size-
extensive) TD CC formalism [32], so it can propagate initial states described
by linear quantum superpositions. In practical TD propagations, SR CC methods
might inherit the numerical challenging technicalities of standard TD CC,
which arise from its non-linear dependencies on the standard TD excitation
vectors. TD EOM-CC theory, in contrast, offers a more direct approach to
propagations that does not rely on the non-linear aspects of standard TD CC
methods, but it is not size-extensive; TD EOM-CC operates similarly as TD
configuration-interaction theory [62, 58]. On the other hand, the non-
Hermiticity of standard TD CC [74] has the undesirable issue of breaking
strict time-reversibility, which is a fundamental physical constraint, but, as
mentioned before, it is a matter that can be resolved in unitary CC theories,
or lessened with an external symmetrization [83]. Finally, the starting
(initial-state) transition elements in our theory requires, in principle, the
application of conventional response theory, which can be computationally
demanding for large systems, but still remains size extensive. For simplicity
we expressed these elements in a sum-over-states fashion, which requires EOM-
CC eigenvectors and eigenvalues, but evidences the fact that these can be
computationally expensive. However, approximations or approaches that
circumvent the need for response-theory relaxation terms for initial
observable matrices, such as dipole matrices, could improve computational
efficiency and be the subject of future work.

##  7 Conclusion

This work examines in depth second response theory, a time-dependent coupled-
cluster formalism to propagate general quantum superpositions, and presents a
simple formula, not developed previously, that predicts accurately the
evolution of a quantum mechanical observable of interest, including
probabilities and coherences. Our theory gives connected expressions to
determine observables in the non-linear regime, and requires the solution to
the standard ground-state propagation problem and the time-dependency of three
response excitation vectors. With these elements provided, we studied its
application to a model three-electron/three-level system under different
initial quantum superpositions. This type of propagation is relevant to
further developing theoretical approaches and the study of systems that are
thermally and/or optically excited. Finally, we examined the use of the SR
energy expressions to approach regularized systems, which do not solve exactly
the ground-state problem, but offer stable approximations to its cluster
amplitudes, excited-state energies, and CC excitation vectors. These could
lead to further improvements to cases where determining the cluster amplitudes
is difficult because of low energy gaps, or related adverse effects.

M.A.M. thanks Montana State University, Bozeman, for startup funding. The
author also acknowledges support by the National Science Foundation through
the MonArk Quantum Foundry, DMR-1906383.

## References

## References

  * [1] Marzari N, Ferretti A and Wolverton C 2021 Nat. Mater. 20 736–749
  * [2] Morton S M, Silverstein D W and Jensen L 2011 Chem. Rev. 111 3962–3994
  * [3] Xie C, Yan D, Chen W, Zou Y, Chen R, Zang S, Wang Y, Yao X and Wang S 2019 Mater. Today 31 47–68
  * [4] Cao Y, Romero J, Olson J P, Degroote M, Johnson P D, Kieferová M, Kivlichan I D, Menke T, Peropadre B, Sawaya N P et al. 2019 Chem. Rev. 119 10856–10915
  * [5] McArdle S, Endo S, Aspuru-Guzik A, Benjamin S C and Yuan X 2020 Rev. Mod. Phys. 92 015003
  * [6] Bauer B, Bravyi S, Motta M and Chan G K L 2020 Chem. Rev. 120 12685–12717
  * [7] Kökcü E, Steckmann T, Wang Y, Freericks J, Dumitrescu E F and Kemper A F 2022 Phys. Rev. Lett. 129 070501
  * [8] Motta M and Rice J E 2022 WIRES Comput. Mol. Sci. 12 e1580
  * [9] Raha K, Peters M B, Wang B, Yu N, Wollacott A M, Westerhoff L M and Merz Jr K M 2007 Drug Discov. Today 12 725–731
  * [10] Cavalli A, Carloni P and Recanatini M 2006 Chem. Rev. 106 3497–3519
  * [11] Saal J E, Kirklin S, Aykol M, Meredig B and Wolverton C 2013 JOM 65 1501–1509
  * [12] Greeley J, Jaramillo T F, Bonde J, Chorkendorff I and Nørskov J K 2006 Nat. Mater. 5 909–913
  * [13] van der Kamp M W and Mulholland A J 2013 Biochem. 52 2708–2728
  * [14] Jones L O, Mosquera M A, Schatz G C and Ratner M A 2020 J. Am. Chem. Soc. 142 3281–3295
  * [15] Baskes M I 1992 Phys. Rev. B 46 2727
  * [16] Lachance-Quirion D, Tabuchi Y, Gloppe A, Usami K and Nakamura Y 2019 Appl. Phys. Express 12 070101
  * [17] Provorse M R and Isborn C M 2016 Int. J. Quantum Chem. 116 739–749
  * [18] Wang Z, Peyton B G and Crawford T D 2022 J. Chem. Theory Comput. 18 5479–5491
  * [19] Palacios A and Martín F 2020 WIRES Comput. Mol. Sci. 10 e1430
  * [20] Samanta P K, Mukherjee D, Hanauer M and Köhn A 2014 J. Chem. Phys. 140 134108
  * [21] Dral P O and Barbatti M 2021 Nat. Rev. Chem. 5 388–405
  * [22] Park J W, Al-Saadon R, MacLeod M K, Shiozaki T and Vlaisavljevich B 2020 Chem. Rev. 120 5878–5909
  * [23] Austin B M, Zubarev D Y and Lester Jr W A 2012 Chem. Rev. 112 263–288
  * [24] Vitillo J G, Cramer C J and Gagliardi L 2022 Isr. J. Chem. 62 e202100136
  * [25] von Lilienfeld O A, Müller K R and Tkatchenko A 2020 Nat. Rev. Chem. 4 347–358
  * [26] Sneskov K and Christiansen O 2012 WIRES Comput. Mol. Sci. 2 566–584
  * [27] Zhang I Y and Grüneis A 2019 Front. Mater. 6 123
  * [28] Liu J and Cheng L 2021 WIRES Comput. Mol. Sci. 11 e1536
  * [29] Bartlett R J and Purvis G D 1978 Int. J. Quantum Chem. 14 561–581
  * [30] Paldus J 1994 Algebraic approach to coupled cluster theory (Boston, MA: Springer US) pp 207–282
  * [31] Bartlett R J 2012 WIRES Comput. Mol. Sci. 2 126–138
  * [32] Koch H and Jørgensen P 1990 J. Chem. Phys. 93 3333–3344
  * [33] Bartlett R J and Musiał M 2007 Rev. Mod. Phys. 79 291–352
  * [34] Taube A G and Bartlett R J 2009 J. Chem. Phys. 130 144112
  * [35] Shen Y, Zhang X, Zhang S, Zhang J N, Yung M H and Kim K 2017 Phys. Rev. A 95 020501
  * [36] Harsha G, Shiozaki T and Scuseria G E 2018 J. Chem. Phys. 148 044107
  * [37] Lee J, Huggins W J, Head-Gordon M and Whaley K B 2018 J. Chem. Theory Comput. 15 311–324
  * [38] Evangelista F A, Chan G K and Scuseria G E 2019 J. Chem. Phys. 151 244112
  * [39] Romero J, Babbush R, McClean J R, Hempel C, Love P J and Aspuru-Guzik A 2018 Quantum Sci. Technol. 4 014008
  * [40] Xia R and Kais S 2020 Quantum Sci. Technol. 6 015001
  * [41] Ryabinkin I G, Yen T C, Genin S N and Izmaylov A F 2018 J. Chem. Theory Comput. 14 6317–6326
  * [42] Anand A, Schleich P, Alperin-Lea S, Jensen P W, Sim S, Díaz-Tinoco M, Kottmann J S, Degroote M, Izmaylov A F and Aspuru-Guzik A 2022 Chem. Soc. Rev. 51 1659–1684
  * [43] Neufeld V A, Ye H Z and Berkelbach T C 2022 J. Phys. Chem. Lett. 13 7497–7503
  * [44] Wang X and Berkelbach T C 2021 J. Chem. Theor. Comput. 17 6387–6394
  * [45] Folkestad S D, Sannes B S and Koch H 2023 The Journal of Chemical Physics 158
  * [46] Ofstad B S, Aurbakken E, Schøyen Ø S, Kristiansen H E, Kvaal S and Pedersen T B 2023 WIRES Comput. Mol. Sci. e1666
  * [47] Skeidsvoll A S, Moitra T, Balbi A, Paul A C, Coriani S and Koch H 2022 Physical Review A 105 023103
  * [48] Pathak H, Sato T and Ishikawa K L 2022 Front. Chem. 10 982120
  * [49] Pathak H, Sato T and Ishikawa K L 2021 J. Chem. Phys. 154 234104
  * [50] Skeidsvoll A S, Balbi A and Koch H 2020 Phys. Rev. A 102 023115
  * [51] Pedersen T B, Kristiansen H E, Bodenstein T, Kvaal S and Schøyen Ø S 2020 J. Chem. Theory Comput. 17 388–404
  * [52] Pedersen T B and Kvaal S 2019 J. Chem. Phys. 150 144106
  * [53] Koulias L N, Williams-Young D B, Nascimento D R, DePrince III A E and Li X 2019 J. Chem. Theory Comput. 15 6617–6624
  * [54] White A F and Chan G K L 2019 J. Chem. Theory Comput. 15 6137–6153
  * [55] White A F and Chan G K L 2018 J. Chem. Theory Comput. 14 5690–5700
  * [56] Park Y C, Perera A and Bartlett R J 2019 J. Chem. Phys. 151 164117
  * [57] Nascimento D R and DePrince III A E 2017 J. Phys. Chem. Lett. 8 2951–2957
  * [58] Nascimento D R and DePrince III A E 2016 J. Chem. Theor. Comput. 12 5834–5840
  * [59] Pigg D A, Hagen G, Nam H and Papenbrock T 2012 Phys. Rev. C 86 014308
  * [60] Wälz G, Kats D, Usvyat D, Korona T and Schütz M 2012 Phys. Rev. A 86 052519
  * [61] Kats D, Usvyat D and Schütz M 2011 Phys. Rev. A 83 062503
  * [62] Sonk J A, Caricato M and Schlegel H B 2011 The Journal of Physical Chemistry A 115 4678–4690
  * [63] Monkhorst H J 1977 Int. J. Quantum Chem. 12 421–432
  * [64] Schönhammer K and Gunnarsson O 1978 Phys. Rev. B 18 6606
  * [65] Hoodbhoy P and Negele J 1978 Phys. Rev. C 18 2380
  * [66] Hoodbhoy P and Negele J 1979 Phys. Rev. C 19 1971
  * [67] Arponen J 1983 Ann. Phys.–New York 151 311–382
  * [68] Mosquera M A 2022 Phys. Rev. A 106 052805
  * [69] Mosquera M A, Chen L X, Ratner M A and Schatz G C 2016 J. Chem. Phys. 144 204105
  * [70] Mosquera M A, Jackson N E, Fauvell T J, Kelley M S, Chen L X, Schatz G C and Ratner M A 2017 J. Am. Chem. Soc. 139 3728–3735
  * [71] Mosquera M A, Jones L O, Kang G, Ratner M A and Schatz G C 2021 J. Phys. Chem. A 125 1093–1102
  * [72] Helgaker T, Jørgensen P and Olsen J 2014 Molecular electronic-structure theory (John Wiley & Sons)
  * [73] Stanton J F and Bartlett R J 1993 J. Chem. Phys. 98 7029–7039
  * [74] Pedersen T B and Koch H 1997 J. Chem. Phys. 106 8059–8072
  * [75] Lawler K V, Parkhill J A and Head-Gordon M 2008 Mol. Phys. 106 2309–2324
  * [76] Stück D and Head-Gordon M 2013 J. Chem. Phys. 139 5203––5219
  * [77] Jacobson G, Marmolejo-Tejada J M and Mosquera M A 2023 ChemPhysChem 24 e202200592
  * [78] Mosquera M A 2021 J. Phys. Chem. A 125 8751–8763
  * [79] Helgaker T and Jørgensen P 1988 Adv. Quantum Chem. 19 183–245
  * [80] Helgaker T, Jørgensen P and Handy N C 1989 Theor. Chim. Acta 76 227–245
  * [81] Kvaal S 2012 J. Chem. Phys. 136 194109
  * [82] Sato T, Pathak H, Orimo Y and Ishikawa K L 2018 J. Chem. Phys. 148 051101
  * [83] Pedersen T B and Koch H 1998 J. Chem. Phys. 108 5194–5204