††thanks: Authors contributed equally to the work††thanks: Authors contributed
equally to the work

# A Synergistic Approach towards Optimization of Coupled Cluster Amplitudes by
Exploiting Dynamical Hierarchy

Chayan Patra  Contributed equally to the work  Department of Chemistry,
Indian Institute of Technology Bombay,
Powai, Mumbai, PIN: 400076, India  Valay Agarawal  Contributed equally to the
work  Department of Chemistry,
University of Chicago,
Chicago, IL 60637, USA  Dipanjali Halder  Department of Chemistry,
Indian Institute of Technology Bombay,
Powai, Mumbai, PIN: 400076, India  Anish Chakraborty  Department of Chemistry,
Indian Institute of Technology Bombay,
Powai, Mumbai, PIN: 400076, India  Dibyendu Mondal  Department of Chemistry,
Indian Institute of Technology Bombay,
Powai, Mumbai, PIN: 400076, India  Sonaldeep Halder  Department of Chemistry,
Indian Institute of Technology Bombay,
Powai, Mumbai, PIN: 400076, India  Rahul Maitra<EMAIL_ADDRESS>Department of Chemistry,
Indian Institute of Technology Bombay,
Powai, Mumbai, PIN: 400076, India

###### Abstract

The coupled cluster iteration scheme for determining the cluster amplitudes
involves a set of nonlinearly coupled difference equations. In the space
spanned by the amplitudes, the set of equations are analysed as a multivariate
time-discrete map where the concept of time appears in an implicit manner.
With the observation that the cluster amplitudes have difference in their
relaxation timescales with respect to the distributions of their magnitudes,
the coupled cluster iteration dynamics are considered as a synergistic motion
of coexisting slow and fast relaxing modes, manifesting a dynamical
hierarchical structure. With the identification of the highly damped auxiliary
amplitudes, their time variation can be neglected compared to the principal
amplitudes which take much longer time to reach the fixed points. We
analytically establish the adiabatic approximation where each of these
auxiliary amplitudes are expressed as unique parametric functions of the
collective principal amplitudes, allowing us to study the optimization with
the latter taken as the independent degrees of freedom. Such decoupling of the
amplitudes significantly reduces the computational scaling without sacrificing
the accuracy in the ground state energy as demonstrated by a number of
challenging molecular applications. A road-map to treat higher order post-
adiabatic effects is also discussed.

##  I Introduction

Coupled Cluster Theory (CC)[1, 2, 3, 4, 5] is well established electronic
structure methodology for accurately solving molecular energetics and
properties. In CC method, a correlated wavefunction is generated by the action
of an exponential wave-operator involving rank-one ($T_{1}$), rank-two
($T_{2}$), …, rank-$n$ ($T_{n}$) cluster operators, that act on a suitably
chosen reference determinant:
$|\psi_{CC}\rangle=e^{T_{1}+T_{2}+T_{3}+...}|\phi_{0}\rangle$. For closed
shell cases, the reference determinant $\phi_{0}$ is often taken as the
Hartree-Fock determinant and the amplitudes ($t$) associated with the cluster
operators are the unknown quantities that are self-consistently optimized. The
cluster amplitudes associated with $T_{k}$ are determined by projecting a
similarity transformed effective hamiltonian
$\hat{H}_{eff}=e^{-T}\hat{H}e^{T}$ against the $k-$tuply excited determinant
that, in principle, folds in the effects of excited determinants such that the
correlated ground state energy can be determined as an expectation value of
$\hat{H}_{eff}$ with respect to the reference determinant.

CC theory is size-extensive and size-consistent at any level of truncation in
the rank of the cluster operators. In this manuscript, without any loss of
generality, we would consider the case where only rank-two cluster operators
are taken. The resulting CC theory with doubles (CCD) involves coupled
optimization of $n_{o}^{2}n_{v}^{2}$ amplitudes, where $n_{o}$ and $n_{v}$ are
the number of occupied (hole, to be denoted as $i,j,k,...$) and unoccupied
(virtual, denoted as $a,b,c,...$) orbitals. Due to the exponential structure
of the waveoperator, any optimization strategy involves nonlinearly coupled
set of equations among $n_{o}^{2}n_{v}^{2}$ amplitudes. This is often done by
iterative minimization of a residue vector
$R_{ijab}=\bra{\phi_{ij}^{ab}}e^{-T}\hat{H}e^{T}\ket{\phi_{0}}$, where
$\ket{\phi_{ij}^{ab}}$ is the excited determinant. Thus, at the fixed point of
the nonlinear optimization process, as $R_{ij}^{ab}\rightarrow 0,\forall
i,j,a,b$, $T\rightarrow T_{CCD}$. Note that due to the excitation structure of
the cluster operators, the cluster operators are only allowed to contract with
the hamiltonian and thus, this optimization procedure scales as
$n_{o}^{2}n_{v}^{4}$ at the worst.

Due to the nonlinear structure of the working equations, it is alluring to
visualize the iterative scheme as a time-discrete multivariate map where
“time” enters in an implicit manner with each iterative step being embedded as
one discrete-time step. Such iteration scheme undergoes chaotic dynamics when
the system of equations is perturbed with an input perturbation. The stability
of the equations were, for the first time, demonstrated by Szakács and
Surján[6, 7] which is later extended to Lippmann-Schwinger equation[8]. Some
of the present authors led the concept one step further to demonstrate that
under the influence of single source of perturbation, the nonlinear iteration
scheme essentially behaves like a multivariate time-discrete map of one-
parameter family[9] which obeys the universality of the Feigenbaum
dynamics[10, 11]. The authors argued the existence of a set of dominant
collective modes that macroscopically govern the optimization process.
Contrarily, there exists large number of other recessive variables
(amplitudes) that essentially evolve synergistically as dictated by the
dominant ones. The present authors, for the first time, exploited the
synchronization among the cluster amplitudes during the optimization
trajectory via machine learning that resulted tremendous reduction in the
degrees of freedom and $40\%-50\%$ savings in computational time to achieve
sub-microHartree ($\mu E_{h}$) [12, 13]. As a matter of fact, such a
synchronous evolution of all the variables during a (discrete) time evolution
is a central theme of Synergetics which gives us a prescription to write the
entire dynamics through the collective dominant modes only. Here we briefly
explain the concept of nonlinear dynamics and synergetics and how one may
apply the concepts of dimensionality reduction for the CC optimization
strategy.

In the regions near the fixed point equilibrium, the linear stability analysis
enables us to classify the variables into two different sets: the auxiliary
(recessive) and principal (dominant) modes[14, 15]. The auxiliary modes are
the ones that relax much faster than the principal modes[16], and in the
characteristic timescale of the principal modes, one can neglect the time
variation of the auxiliary modes (known as Adiabatic Approximation)[17]. This
is conceptually the generalised version of the Born-Oppenheimer
approximation[18]. For most of the practical cases the following conditions
hold[19]:

  1. 1.

the number of principal modes are much less than than the number of auxiliary
modes

  2. 2.

the magnitude of the principal modes are significantly larger than the
magnitude of the auxiliary modes.

Given that the system of variables can be grouped into the auxiliary and
principal modes, the “Slaving Principle” enables us to express each of the
auxiliary modes as unique function of the principal modes. This thus allows us
to confine our attention only to the dynamics of the principal modes. The
auxiliary modes unanimously follow the “orders” of the principal modes and
therefore, this is known as the master-slave dynamics. With the time evolution
of the principal modes, the auxiliary modes also get updated according to
their parametric dependence on the principal modes. However, there exists
circular causal relationship such that the updated information of the
auxiliary modes gets coupled back to the dynamical equations of the principal
modes, which may be referred to as the feedback coupling[19]. The whole
dynamical system thus evolves in the interdependent way. As the auxiliary
modes at each step merely chip in as a parametric function of the principal
modes, the effective degrees of freedom of the whole system is significantly
reduced. In one of our previous publications, we developed a first-principle
based CC optimization strategy, termed as the Adiabatically Decoupled Coupled
Cluster (ADCC), where we have shown that it is possible to distinguish the
various amplitudes at their MP2 level into principal and auxiliary modes.
Based on certain assumptions which are conceptually close in spirit to the
adiabatic approximation, the authors were able to reproduce very accurate
results by optimizing the variables in a reduced subspace[20].

The main purpose of this article is to show that some of these features of
nonlinear dynamics for the systems close to classical critical points can also
be applied in numerically accurate manner (although not mathematically exact)
to reduce the computational scaling where nonlinear iterative optimization is
involved. We will also demonstrate how the “adiabatic approximation” comes out
quite naturally in a more general and rigorous mathematical way, given that
some of the cluster amplitudes have significantly longer time scale of
relaxation to converge. This will also elucidate the scope to further improve
the ADCC results by incorporating some of the higher order post-adiabatic
terms[21, 22]. Towards this, we will first present mathematical preliminaries,
taking general class of dynamical systems that motivate us towards the
development of the adiabatically decoupled CC iterative formulation. We will
mainly focus on the underlying concepts of master and slave variables, which
are related by the slaving principle, and will show how the CC iterative
scheme can be adapted within the master-slave dynamics framework for
dimensionality reduction. We will show the performance of our model as a
function of the dimension of the master modes to justify that CC theory in
principle can be optimized in significantly reduced dimension than that
dictated by the size of the basis set.

##  II Adiabatically Decoupled Scheme for Coupled Cluster Amplitude
Optimization:

###  II.1 General Mathematical Preliminaries towards Dimensionality Reduction:
Concepts of Adiabatic Approximation and Slaving Principle

As mentioned in the introduction, for specific multivariate dynamical systems,
one may divide the entire variable space into auxiliary and principal modes
having large difference in their magnitude and characteristic timescale of
relaxation. In order to have an easy readability, in the following, we first
discuss the mathematical preliminaries that allow us to express the dynamics
of the whole system solely in terms of the principal modes. We will consider a
general dynamical system with continuous time variation. The corresponding
discrete analogue pertinent to the CC iterative scheme will be developed in
the next subsection.

Let us consider a multivariate dynamical system $\\{\mathbf{\hat{q}}\\}$
evolving according to the following general form of an equation:

$\mathbf{\dot{\underline{\hat{q}}}}=L\mathbf{\underline{\hat{q}}}+\mathbf{N}(\\{\mathbf{\hat{q}}\\})$
(1)

where, $L$ is a coefficient matrix independent of the dynamical variables and
$\mathbf{N}$ contains all the nonlinear couplings among them. One may perform
a linear stability analysis[23, 24, 25] to obtain a set of eigenvalues
$\\{\lambda\\}$ of the stability matrix that can be classified into the
following two categories:[22, 19]

$\begin{split}&\\{{\lambda}_{u}\\}\gtrless 0;\hskip 14.22636pt\mbox{(total
number }N_{u})\\\ &\\{{\lambda}_{s}\\}<0;\hskip 17.07164pt\mbox{(total number
}N_{s})\end{split}$ (2)

with the condition:

$\mid\lambda_{s_{i}}\mid>\mid\lambda_{u_{I}}\mid$ (3)

In the usual scenario, $N_{s}>>N_{u}$, where the individual elements from the
sets of eigenvectors of $\lambda_{u}$ ($\\{u\\}$) and $\lambda_{s}$
($\\{s\\}$) are denoted by $u_{I}$ and $s_{i}$. Positive values of
$\lambda_{u}$ indicates the dynamics is non-equilibrium and the corresponding
eigenvectors are the unstable modes which move away from the fixed point as
time progresses[24]. This justifies the subscript $u$ and $s$ referring to
unstable and stable modes, respectively. On the other hand, keeping in mind
the CC iteration dynamics, we are only interested here in the negative
$\lambda_{u}$ cases, referring to an overall equilibrium system that
eventually reaches one of the fixed points of the dynamics. This indicates
that all the $\lambda$s can be considered as damping factors. Hence, we drop
the terminology of unstable and stable modes from here on, rather we keep the
same subscripts but prefer to refer them in a more general way as the
principal and auxiliary modes. In the transient period, the variables with
smaller and larger amplitudes generally have larger and smaller damping
factors respectively, allowing us to establish a demarcation between auxiliary
and principal amplitudes. The nonlinearly coupled equations of motion for
these two sets of amplitudes can be written as:

$\dot{u_{I}}=\lambda_{u_{I}}u_{I}+Q_{I}({u,s});\hskip 14.22636pt\mbox{where
}I=1,2,...,N_{u}$ (4)

and

$\dot{s_{i}}=\lambda_{s_{i}}s_{i}+P_{i}({u,s});\hskip 14.22636pt\mbox{where
}i=1,2,..,N_{s}$ (5)

Here $Q_{I}$, $P_{i}$ contains all the information of nonlinearity and inter-
mode coupling. In a narrow region around a fixed point (i.e.
$\\{\dot{u_{I}}\\},\\{\dot{s_{i}}\\}=0$) where linearization can be applied,
one may analyze the dynamics by neglecting the small contributions from the
nonlinear terms in $Q_{I}$ and $P_{i}$. In such a scenario, when the system is
perturbed slightly from the fixed point equilibrium, the whole system moves
under the influence of the principal modes $\\{u_{I}\\}$ in their
characteristic timescale, whereas, the auxiliary modes $\\{s_{i}\\}$ decay
back to the fixed point. In other words, $s_{i}$ can be considered as heavily
damped modes whereas $u_{I}$ is relatively under-damped.

By exploiting the characteristic relaxation timescale, we wish to show that
$s_{i}$ amplitudes can be expressed by means of $u_{I}$ only, such that one
can also eliminate $\\{s_{i}\\}$ from Eq. (4) by direct substitution. The most
general solution to Eq. (5) can be obtained as

$\begin{split}&\Big{(}\frac{d}{dt}-\lambda_{s_{i}}\Big{)}s_{i}=P({u,s})\\\
\implies&s_{i}=\Big{(}\frac{d}{dt}-\lambda_{s_{i}}\Big{)}^{-1}P(\\{u,s\\})\\\
\end{split}$ (6)

The inverse differential operator can be solved with an integral form that
allows us to write[22]:

$s_{i}(t)=\int_{-\infty}^{t}e^{\lambda_{s_{i}}(t-\tau)}P(\\{u,s\\})d{\tau}$
(7)

Assuming the integral is finite and well-behaved[26, 27], one may perform an
integration by parts to further simplify Eq. (7)

$s_{i}(t)=-\frac{1}{\lambda_{s_{i}}}P({\\{u,s\\}})-\frac{1}{(-\lambda_{s_{i}})}\int_{-\infty}^{t}e^{\lambda_{s_{i}}(t-\tau)}(\frac{d}{d\tau}P)d{\tau}$
(8)

The second term in right hand side onward are usually very small. Hence,
neglecting the small higher order terms, Eq. (8) reads

$s_{i}(t)=-\frac{1}{\lambda_{s_{i}}}P(\\{u,s\\})$ (9)

A careful observation leads to the interesting fact that we can avoid all the
mathematical jargon and still get the same result as that of Eq. (9) if we
simply set $\dot{s_{i}}=0$ in Eq. (5). This simply implies that in the
characteristic time scale of the principal modes, the time-variation of the
auxiliary modes can be neglected. This is to be referred as the adiabatic
approximation.

Generally in all practical cases, the amplitude of the principal modes are
extremely large in comparison to auxiliary modes. Thus all the contributions
from $s_{i}$ can be neglected compared to $u_{I}$ such that the right hand
side of Eq. (9) can be written as a function of $u_{I}$ alone, which
simplifies the expression as:

$s_{i}(t)=-\frac{1}{\lambda_{s_{i}}}P(\\{u\\})$ (10)

As we have the expression in Eq. (10) for the auxiliary modes as a function of
the principal modes, we can now imitate the system dynamics in a reduced
subspace spanned only by the principal modes, which are much fewer in number
compared to the auxiliary modes. In the context of CC iterative approach, our
goal is to accurately determine $s$ as a function of $\\{u\\}$ only, which
would be fed back to the dynamics of $\\{u\\}$ to get an updated set of
principal modes. This loop, referred to as the circular causality loop, would
continue till the principal modes converge below a pre-defined threshold.

###  II.2 Coupled Cluster iterative optimization from the synergistic
viewpoint:

We now turn our attention to the optimization of CC theory using Jacobi
iterative scheme. We would consider the iterative scheme as a time-discrete
map where each iteration is embedded as one time step. Furthermore, since the
iterative optimization leads to a set of fixed points, it can be considered as
an equilibrium system from the perspective of nonlinear science, and hence the
demarcation of the auxiliary and principal modes (as shown in Eq. (2) for non-
equilibrium systems) does not strictly hold according to the linear stability
analysis[9]. In general, for equilibrium systems (in the perspective of
nonlinear dynamics), all the values of $\lambda$ may be negative; however, for
all practical purposes, one may bypass the linear stability criteria and
classify the amplitudes based on their distribution of magnitudes into
auxiliary and principal modes having shorter or longer characteristic
relaxation time to reach their fixed point solutions. We will generally refer
to these amplitudes as auxiliary and principal amplitudes, respectively.

With the approximations in mind, let us now consider the amplitudes updating
Jacobi equation

$\begin{split}\Delta
t_{\mu}&=t_{\mu}^{(k+1)}-t_{\mu}^{k}=\frac{R_{\mu}}{D_{\mu}}\\\
&=\frac{1}{D_{\mu}}(H_{\mu}+(\contraction{}{H}{}{T}HT)_{\mu}+\frac{1}{2}(\contraction{}{H}{}{T}\contraction[2ex]{}{H\;}{T}{}H\;TT)_{\mu})\end{split}$
(11)

Here $T$ is the cluster operator, $t_{\mu}$ denotes the associated amplitudes
with hole-particle excitation structure ‘$\mu$’ and $R_{\mu}$ is the
corresponding residue. The hamiltonian $H$ contains the one and two electron
terms. It was observed from a number of numerical examples[20] that the
amplitudes with large magnitude (at the MP2 level) take substantially more
number of iterations to reach their converged values (fixed points) than the
ones with smaller magnitudes. Based on the relative magnitude of these
amplitudes, the entire amplitude space can be subdivided into a Large
Amplitude Subset (LAS, spanned by $\\{T_{L}\\}$ with dimension $n_{L}$) and
Small Amplitude Subset (SAS, spanned by $\\{T_{S}\\}$, with dimension
$n_{S}$). In our time-discrete iterative map, their amplitudes $\\{t_{L}\\}$
plays the role of the principal modes and $\\{t_{S}\\}$ plays the role of the
auxiliary modes. Note also that $n_{S}>>n_{L}$ and
$\mid{t_{L}}\mid>\mid{t_{S}}\mid$. They will thus be referred to as the
principal and auxiliary amplitudes, respectively. In conjunction to the
dynamics of the principal and auxiliary modes shown in Eqs. (4) and (5), the
same for the principal and auxiliary amplitudes can be written as:

$\begin{split}\Delta
t_{L_{I}}&=\frac{1}{D_{L_{I}}}{(\contraction{}{H^{d}_{L_{I}}}{}{T_{L_{I}}}H^{d}_{L_{I}}T_{L_{I}})}_{L_{I}}+Q_{L_{I}}(\\{T_{L},T_{S}\\});\forall
I=1,2,...,n_{L}\end{split}$ (12)

and

$\begin{split}\Delta
t_{S_{i}}&=\frac{1}{D_{S_{i}}}{(\contraction{}{H^{d}_{S_{i}}}{}{T_{S_{i}}}H^{d}_{S_{i}}T_{S_{i}})}_{S_{i}}+P_{S_{i}}(\\{T_{L},T_{S}\\});\forall
i=1,2,...,n_{S}\end{split}$ (13)

where, $L_{I}$ and $S_{i}$ are composite hole-particle indices for LAS and SAS
elements, respectively. $D_{\mu}$ is the corresponding orbital energy
difference for the general hole-particle index $\mu$ . Here $Q_{L_{I}}$ and
$P_{S_{i}}$ contain all the information of nonlinearity for the dynamics of
the principal and auxiliary amplitudes, and they can be expanded as:

$\begin{split}&Q_{L_{I}}(\\{T_{L},T_{S}\\})=\frac{1}{D_{L_{I}}}\Big{[}\Big{(}H_{L_{I}}+{{(\contraction{}{H}{}{T_{L_{J}}}HT_{L_{J}})}_{L_{I}}}+{{(\contraction{}{H}{}{T_{S_{i}}}HT_{S_{i}})}_{L_{I}}}\Big{)}+\\\
&\frac{1}{2}\Big{(}{(\contraction{}{H}{}{T_{L_{J}}}\contraction[2ex]{}{H\;}{T_{L_{K}}}{}H\;T_{L_{J}}T_{L_{K}})_{L_{I}}}+{(\contraction{}{H}{}{T_{L_{J}}}\contraction[2ex]{}{H\;}{T_{S_{i}}}{}H\;T_{L_{J}}T_{S_{i}})_{L_{I}}}+{(\contraction{}{H}{}{T_{S_{i}}}\contraction[2ex]{}{H\;}{T_{S_{j}}}{}H\;T_{S_{i}}T_{S_{j}})_{L_{I}}})\Big{)}\Big{]}\end{split}$
(14)

and

$\begin{split}&P_{S_{i}}(\\{T_{L},T_{S}\\})=\frac{1}{D_{S_{i}}}\Big{[}\Big{(}H_{S_{i}}+{{(\contraction{}{H}{}{T_{L_{I}}}HT_{L_{I}})}_{S_{i}}}+{{(\contraction{}{H}{}{T_{S_{j}}}HT_{S_{j}})}_{S_{i}}}\Big{)}+\\\
&\frac{1}{2}\Big{(}{(\contraction{}{H}{}{T_{L_{I}}}\contraction[2ex]{}{H\;}{T_{L_{J}}}{}H\;T_{L_{I}}T_{L_{J}})_{S_{i}}}+{(\contraction{}{H}{}{T_{L_{I}}}\contraction[2ex]{}{H\;}{T_{S_{j}}}{}H\;T_{L_{I}}T_{S_{j}})_{S_{i}}}+{(\contraction{}{H}{}{T_{S_{j}}}\contraction[2ex]{}{H\;}{T_{S_{k}}}{}H\;T_{S_{j}}T_{S_{k}})_{S_{i}}})\Big{)}\Big{]}\end{split}$
(15)

Figure 1: Representative diagrams contributing to the leading order of
computational scaling. Diagrams of the same topology are shown for forward
mapping (a) and backward mapping (b) with different interpretation.

Here the diagonal terms, $H^{d}$, in Eqs. (12) and (13) needs some
explanation. The diagonal part of the hamiltonian matrix $H^{d}_{\mu}$ are the
composite tensorial terms which upon contraction with $T_{\mu}$ generate a
structure with hole-particle composite structure as $\mu$. The diagonal part
is extracted out in the expressions of $\Delta t_{\mu}$ in Eqs. (12) and (13)
to show their one-to-one correspondence with Eqs. (4) and (5). The explicit
expression for $H^{d}_{\mu}$ is given by

$H_{\mu}^{d}=(1+P(i,j)P(a,b))\Big{(}-f_{ii}+f_{aa}+\frac{1}{2}v_{ab}^{ab}+2v_{ia}^{ai}\\\
-(1+\delta_{ij}\delta_{ab})v_{ia}^{ia}+\frac{1}{2}v_{ij}^{ij}-v_{ib}^{ib}\Big{)}\\\
\mu=\\{ijab\\}\in\mbox{{SAS} or {LAS}}$ (16)

Here $i,j$ refer to the occupied (hole) orbital indices while $a,b$ refer to
the unoccupied (particle) orbital indices. The canonical Fock operator matrix
elements and the two electron integrals are represented by f and v
respectively, and $P$ is the permutation operator.

We are now in a position to develop the discrete-time analogue of Eq. (7) for
the CC optimization. Without going into the details of the equation, for which
we refer to Haken[28, 15], the most general solution for the auxiliary
amplitudes, $t_{S}$, can be written as:

$t_{S_{i}}=\sum_{m=-\infty}^{l}(1+\lambda_{S_{i}})^{l-m}P_{S_{i}}$ (17)

where,

${\lambda_{S_{i}}=\frac{H^{d}_{S_{i}}}{D_{S_{i}}}}$ (18)

Note that $H^{d}_{S_{i}}$ and $D_{S_{i}}$ have opposite sign (and hence CC
equations form a convergent series) and thus $\lambda_{S_{i}}$ is negative.
The expression in the right hand side of Eq. (18) can be expanded using
summation by parts, a mathematical trick similar to integration by parts of
Eq. (8), to obtain:

$\begin{split}&t_{S_{i}}=P_{S_{i}}\sum_{m=-\infty}^{l}(1+\lambda_{S_{i}})^{l-m}\\\
&-\sum_{m=-\infty}^{l}(1+\lambda_{S_{i}})^{l+1-m}\sum_{m^{\prime}=-\infty}^{m-1}(1+\lambda_{S_{i}})^{m-1-m^{\prime}}\Delta
P_{S_{i}}\end{split}$ (19)

Here, “$\Delta$” is the time-discrete version of the time derivative operator
“$\nicefrac{{d}}{{dt}}$”. Since the auxiliary modes have
$\lambda_{S_{i}}<0,\forall S_{i}$, one can simplify the intricate summations
above using the identity:

$\sum_{m=-\infty}^{l}(1+x)^{l-m}=-\frac{1}{x};\hskip 2.84526pt\mid x+1\mid<1$
(20)

The simplification of Eq. (19) leads to

$t_{S_{i}}=\underbracket{-\frac{P_{S_{i}}(\\{t_{L},t_{S}\\})}{\lambda_{S_{i}}}}_{\text{\clap{adiabatic}}}\overbracket{-\frac{\Delta
P_{S_{i}}(\\{t_{L},t_{S}\\})}{\lambda_{S_{i}}^{2}}}^{\text{\clap{post-
adiabatic}}}$ (21)

We will briefly explain the significance of these two terms that appear on the
right hand side of Eq. (21). As mentioned previously, one may neglect the time
variation of the auxiliary modes (auxiliary amplitudes) in the characteristic
timescale of the principal modes (principal amplitudes). By setting up the
condition $\Delta{t_{S_{i}}}=0$ to neglect the time variation of the auxiliary
amplitudes, one obtains the first term on the right hand side of Eq. (21).
This thus can be interpreted as the adiabatic approximation. Note that in one
of our earlier publications[20], the adiabatic approximation was introduced
simply by setting the decoupling condition. Here we further corroborate the
concept by introducing a discrete-time dependent picture and timescale
decoupling in the amplitude space. The second term on the right hand side of
Eq. (21) depends on the discrete-time variation of the (part of the) residue
$P_{S_{i}}$ and depends on $\lambda_{S_{i}}^{2}$. This thus may be interpreted
as post- adiabatic correction and will not be considered in this article any
further. We will introduce the post-adiabatically corrected optimization
scheme in our future publication.

In what follows, we will now proceed exactly the same way that allowed us to
derive Eq. (10) from Eq. (9). Noting the fact[20] the auxiliary amplitudes
$t_{S}$ are significantly small compared to the principal amplitudes, $t_{L}$,
one may neglect $t_{S}$ from the adiabatically decoupled expression of
$t_{S_{i}}$:

$t_{S_{i}(ad)}=-\frac{P_{s_{i}}(\\{t_{S},t_{L}\\})}{\lambda_{S_{i}}}\xrightarrow{{\mid
t_{S}\mid}\approx 0}-\frac{P_{s_{i}}(\\{t_{L}\\})}{\lambda_{S_{i}}}$ (22)

or in the long hand notation by writing $P_{S_{i}}(\\{t_{L}\\})$ explicitly:

$t_{S_{i}(ad)}=-\frac{1}{H^{d}_{S_{i}}}\Big{[}H_{S_{i}}+{{(\contraction{}{H}{}{T_{L_{I}}}HT_{L_{I}})}_{S_{i}}}+\frac{1}{2}\Big{(}{(\contraction{}{H}{}{T_{L_{I}}}\contraction[2ex]{}{H\;}{T_{L_{J}}}{}H\;T_{L_{I}}T_{L_{J}})_{S_{i}}}\Big{)}\Big{]}$
(23)

where all the terms that involve $T_{S}$ are set to zero. We note that setting
the adiabatic condition $\Delta{t_{S_{i}}}=0$ does not imply that the
auxiliary amplitudes are frozen to given pre-computed values. Rather, they are
updated through parametric functional dependence on the principal amplitudes.
We now turn our attention to the computation of the determination of the
principal amplitudes via feedback coupling.

Figure 2: A schematic representation of the circular algorithm where the
auxiliary amplitudes, $t_{S}$ are determined through sole function of the
principal amplitudes, $t_{L}$, which are in turn updated through the back
coupling of both sets of amplitudes. At each of these steps, the computational
scaling goes as $n_{L}n_{v}^{2}<<n_{o}^{2}n_{v}^{4}$.

#### Feedback coupling and the determination of the principal amplitudes:

As previously mentioned, there exists circular causality relationship among
the principal and the auxiliary amplitudes. This implied that the updated
information of the auxiliary amplitudes (obtained from Eq. (23)) must get fed
back to the updating equations for $T_{L}$. Noting the fact that
$T_{s_{i}(ad)}=T_{S_{i}}(\\{T_{L}\\})$, one may explicitly write the principal
amplitude update equations as:

$\begin{split}&\Delta
t_{L_{I}}=\frac{1}{D_{L_{I}}}\Big{[}\Big{(}H_{L_{I}}+{{(\contraction{}{H}{}{T_{L_{J}}}HT_{L_{J}})}_{L_{I}}}+{{(\contraction{}{H}{}{T_{S_{i}(ad)}}HT_{S_{i}(ad)})}_{L_{I}}}\Big{)}+\\\
&\frac{1}{2}\Big{(}{(\contraction{}{H}{}{T_{L_{J}}}\contraction[2ex]{}{H\;}{T_{L_{K}}}{}H\;T_{L_{J}}T_{L_{K}})_{L_{I}}}+{(\contraction{}{H}{}{T_{L_{J}}}\contraction[2ex]{}{H\;}{T_{S_{i}(ad)}}{}H\;T_{L_{J}}T_{S_{i}(ad)})_{L_{I}}}+\\\
&{(\contraction{}{H}{}{T_{S_{i}(ad)}}\contraction[2ex]{}{H\;}{T_{S_{j}(ad)}}{}H\;T_{S_{i}(ad)}T_{S_{j}(ad)})_{L_{I}}})\Big{)}\Big{]}\end{split}$
(24)

which can be written as:

$\begin{split}&\Delta
t_{L_{I}}=\frac{1}{D_{L_{I}}}\Big{(}H_{L_{I}}+{{(\contraction{}{H}{}{T}HT)}_{L_{I}}}+\frac{1}{2}{(\contraction{}{H}{}{T}\contraction[2ex]{}{H\;}{T}{}H\;TT)_{L_{I}}}\Big{)}\end{split}$
(25)

where

$T=\sum_{I}T_{L_{I}}\oplus\sum_{i}T_{S_{i}(ad)}(\\{T_{L}\\})$ (26)

This implies that both the principal and auxiliary amplitudes couple to the
equations of the former. Since the contribution of the auxiliary amplitudes is
significantly smaller, one may truncate the right hand side of Eq. (24) by
keeping only those terms containing $T_{S}$ which are overall linear.
Depending on which terms are retained, we have developed two schemes. In
scheme-I, all the terms in the right hand side of Eq. (24) are retained. This
implies that scheme-I takes care of complete feedback coupling. In scheme-II,
only terms I-IV in the right hand side of Eq. (24) are retained. In other
words, the feedback of the auxiliary amplitudes in the equation of the
principal amplitudes is taken up to overall linear order.

#### Analysis of computational scaling:

Figure 3: Comparative study of the accuracy of ADCCD scheme-I (black line) and
scheme-II (red line) against canonical CCD as functions of LAS size for four
different molecules. Note that the horizontal blue line denotes the canonical
CCD energy scaled to zero and the blue shaded region denotes $\pm 0.5mE_{h}$
region. ADCCD scheme-I is slightly more accurate over scheme-II, particularly
with smaller LAS, as expected. ADCCD scheme-I achieves $\sim 0.01mE_{h}$
accuracy with only about 30 % of the total nonzero amplitudes taken as
independent (principal) ones. The inset shows the accuracy of the combined set
of cluster amplitudes as compared to canonical CCD and the former is shown to
be as accurate as the latter with LAS taken to be 40$\%$ where one obtains
accuracy up to $\mu E_{h}$ precision.

We now briefly discuss about the computational scaling associated with the
determination of the auxiliary amplitudes (Eq. (23)) and principal amplitudes
(Eq. (24)). For detailed analysis of the computational scaling, we refer to
one of our earlier publications[20]. Here, for our analysis of the scaling, we
would only consider the most expensive linear diagram that appears in CC
theory, namely the one with all particle contraction. Such a representative
diagram is shown in Fig. 1. Fig. 1(a) represents its interpretation in the
context of the forward mapping via Eq. (23) where $t_{L}$ amplitudes uniquely
determine the $t_{S}$ amplitudes. Diagrammatically, only the $T_{L}$ operators
(filled squared box, of dimension $n_{L}$) contracts with the hamiltonian
(filled circles) to generate $T_{S}$. Here the index quartet of $T_{L}$,
$(ijcd)$, necessarily belongs to one of the $n_{L}$ elements of LAS. On the
other hand, the uncontracted particle indices $a,b$ can be arbitrary with the
constraint that the uncontracted index quartet $(ijab)$ should necessarily
belong to one of the SAS elements. This automatically renders the leading
computational scaling for the forward mapping to be $n_{L}n_{v}^{2}$.

We now turn our attention to the scaling of the backward mapping (Eq. (24)).
For this purpose, we resort to the diagram of the same topology as we
discussed for the forward mapping, and interpret it in a different way as
shown in Fig. 1(b). In this case, all the amplitudes are allowed to couple
with the hamiltonian, however, the resulting structure should have the
corresponding amplitude in LAS. This implies that the index quartet $(ijcd)$
can be any one of the principal or auxiliary amplitudes. However, the
uncontracted index quartet $(ijab)$ should necessarily belong to one of the
$n_{L}$ elements of $T_{L}$. The restriction on the uncontracted index to the
ones belonging to LAS automatically ensures that $i$ and $j$ cannot take any
arbitrary hole orbital level. The scaling for the construction of this diagram
is again $n_{L}n_{v}^{2}$ at the worst.

Thus the overall scaling of the scheme never exceeds $n_{L}n_{v}^{2}$ for the
ADCC scheme, whereas, for CCD, the same scales as $n_{o}^{2}n_{v}^{4}$. Since
$n_{L}$ is only a small fraction (vide infra) compared to
$n_{o}^{2}n_{v}^{2}$, there is a clear computational advantage of our scheme
offers over the conventional one.

In summary, the overall algorithm moves in a circular manner where at each
step, the auxiliary amplitudes are determined via Eq. (23) as functions of
$\\{t_{L}\\}$ alone. On the other hand, the principal amplitudes, $t_{L}$ are
updated through complete or partial feedback coupling (for scheme-I and
scheme-II, respectively) via Eq. (24). A schematic representation of the
algorithm is shown in Fig. 2. Note that at each step, the scaling goes as
$n_{L}n_{v}^{2}$. With $n_{L}<<n_{o}^{2}n_{v}^{2}$, our scheme offers
significant computational savings over the conventional CC algorithm.

##  III Results and Discussion :

In this section, we will discuss the efficiency and accuracy of the ADCC
formalism and will compare it against the canonical coupled cluster method. As
a proof-of-concept, without any loss of generality, we will restrict the
cluster amplitudes to only doubles (CCD). All the results are obtained with
our in-house codes and the convergence threshold was set to be $10^{-6}$.
Furthermore, no DIIS optimization is used to accelerate the convergence, even
though inclusion of it is fairly straightforward.

Figure 4: Accuracy of ADCCD scheme-I for various size of LAS as function of
the cumulative cost of the ADCCD calculation obtained from the highest scaling
term. The blue line shows the cumulative number of matrix operation for the
construction of the diagram containing all particle contraction for canonical
CCD.

As our objective is to simulate the canonical CCD optimization in a reduced
subspace with order of magnitude less computational scaling, it is imperative
to take the same as our reference. In Fig. 3, we have plotted the difference
in energy obtained with ADCCD (both scheme-I and scheme-II) and canonical CCD
as a function of the size of LAS. Note that the percentage of LAS (N) in the
horizontal axis signifies that the largest $N\%$ of total nonzero amplitudes
(at the MP2 level) are taken as independent parameters. Thus, a larger LAS
signifies more number of principal amplitudes taken as the independent
variable, against which all the remaining amplitudes (the auxiliary ones) are
mapped. Note that the horizontal blue line (y=0) denotes the exact CCD results
and the light-blue shaded region indicates a $0.5$ milli-Hartree ($mE_{h}$)
energy band vertically each sides of the reference energy line that provides a
measure of the accuracy. Four different molecules were chosen for our
application which have distinct variety of electronic complexity: (a)
symmetric water (equilibrium geometry), (b) asymmetric water, (c) ammonia at
equilibrium geometry and (d) formaldehyde[29]. As we expand LAS, the energy
calculated from ADCCD (both scheme-I and scheme-II) converge monotonically to
the corresponding CCD values. Understandably, ADCCD scheme-I is somewhat
better than ADCCD scheme-II due to the fact that the former takes into account
the complete feedback coupling of the auxiliary amplitudes towards the
principal amplitudes, while the later does so in an approximate manner. For
most of the cases, ADCCD Scheme-I enters the shaded region in between
$20\%$-$25\%$ of LAS and acquires a remarkable accuracy of $\sim 0.1mE_{h}$ at
around $30\%$ in almost all the cases under consideration, while ADCCD scheme-
II does so around $25\%$-$30\%$ and $40\%$ of the LAS percentage,
respectively. This implies that with $70\%$ reduction in the number of
independent amplitudes, ADCCD scheme-I can achieve an accuracy within tens of
$\mu E_{h}$ as compared to the canonical CCD, while ADCCD scheme-II does so
with $60\%$ reduction in the number of amplitudes. We will soon analyse how
these tremendous reduction in the number of independent parameters transcends
into reduction in the number of matrix operations needed to achieve the
accuracy.

In the inset of each subfigure of Fig. 3, we have plotted the overlap between
the adiabatically determined cluster amplitudes and the canonical cluster
amplitudes. Towards this, we have mapped the amplitudes as normalized column
vectors $\ket{t_{CCD}}$ and $\ket{t_{ADCCD}}$. Note that for ADCC schemes,
both the principal and auxiliary amplitudes are included in $\ket{t_{ADCCD}}$
where their index-ordering were kept consistent. The accuracy of the ADCCD
cluster amplitudes are measured by the quantity
$(1-\bra{t_{ADCCD}}\ket{t_{CCD}})$. The inset plot shows the monotonic
approach of the overlap value to one for scheme-I with increase in the size of
LAS, signifying that the ADCCD cluster amplitudes are of the similar accuracy
to those obtained by the canonical CCD scheme.

Having discussed the precision of the optimization algorithm, we now comment
on the number of matrix operations needed. While canonical CCD requires
$\mathcal{O}(n_{o}^{2}n_{v}^{4})$ matrix operations each iterative cycle, our
scheme requires $\mathcal{O}(n_{L}n_{v}^{2})$ matrix operations. In Fig.4, we
have demonstrated that one can achieve precision of the order of tens of $\mu
E_{h}$ with LAS taken to be $30\%$ while it takes at least an order fewer
number of cumulative matrix operations. One may note that the ADCC scheme, at
each iterative cycle, requires one order of magnitude less memory compared to
canonical CCD, and thus enabling large scale computations without
significantly sacrificing the accuracy.

##  IV Conclusion and Future Outlook :

In this article we have rigorously established the fact that the discrete-time
CC iteration dynamics can be considered as a collective motion of coexisting
slow and fast relaxing modes, manifesting a dynamical hierarchical structure.
Although the concepts of adiabatic elimination and slaving principle in
synergetics and nonlinear dynamics are applicable for non-equilibrium systems
in most of the existing literatures, borrowing these ideas we have
successfully developed an equilibrium dynamical formalism (in the amplitude
space) to study the complex transient behaviour of coupled cluster iteration
before convergence. Our scheme shows order of magnitude reduction in
computational efforts to obtain molecular ground state energy with
sub-$mE_{h}$ to tens of $\mu E_{h}$ precision to the conventional CC schemes.
The results are, furthermore, systematically improvable with inclusion of more
number of amplitudes in the LAS. However while doing so we have neglected some
higher order post-adiabatic terms that opens up scopes for further improvement
of the adiabatic results[21]. Also, so far we have hived off LAS from the full
cluster amplitude space purely based on the observation that the cluster
amplitudes have a difference in relaxation timescales with respect to their
magnitude distribution. There exist some information theoretical techniques
[30, 31, 32, 33, 34, 35] and entropy maximization approaches [36, 37] to
accurately determine the principal amplitudes. Although most of these
principles are mathematically heavy and applicable for thermodynamical systems
far from equilibrium, we are looking for devising our own methods based on
these concepts to construct LAS in a more mathematically accurate manner. Some
stochastic adiabatic elimination principles[38, 39, 40, 41, 42, 43, 44] can
also be applied to introduce stochasticity or external noise in CC to study
thermal properties. We can utilise some studies that have shown how the stable
and unstable manifolds can suddenly change or collide among themselves in the
phase portrait under certain conditions to show some interesting phenomena[45,
46]. In the immediate future, we will try to incorporate the higher order
post-adiabatic correction terms to investigate further the intricate
mathematical structures of the theory and to improve the results on top of
ADCC.

##  V Acknowledgments

RM acknowledges the financial support from Industrial Research and Consultancy
Centre, IIT Bombay, and Science and Engineering Research Board, Government of
India.

## Data Availability

The data is available upon reasonable request to the corresponding author.

## Keywords:

Electronic structure. Ab initio calculations. Quantum Chemistry. Nonlinear
Dynamics. Synergetics.

## References

  * C̆íz̆ek [1966] J. C̆íz̆ek, “On the correlation problem in atomic and molecular systems. calculation of wavefunction components in ursell-type expansion using quantum-field theoretical methods,” J. Chem. Phys. 45, 4256–4266 (1966).
  * C̆íz̆ek [1969] J. C̆íz̆ek, “On the use of the cluster expansion and the technique of diagrams in calculations of correlation effects in atoms and molecules,” Adv. Chem. Phys. 14, 35–89 (1969).
  * Čížek and Paldus [1971] J. Čížek and J. Paldus, “Correlation problems in atomic and molecular systems iii. rederivation of the coupled-pair many-electron theory using the traditional quantum chemical methods,” Int. J. Quantum Chem. 5, 359–379 (1971).
  * Bartlett and Musiał [2007] R. J. Bartlett and M. Musiał, “Coupled-cluster theory in quantum chemistry,” Reviews of Modern Physics 79, 291 (2007).
  * Crawford and Schaefer [2000] T. D. Crawford and H. F. Schaefer, “An introduction to coupled cluster theory for computational chemists,” Reviews in computational chemistry 14, 33–136 (2000).
  * Szakács and Surján [2008a] P. Szakács and P. R. Surján, “Stability conditions for the coupled cluster equations,” Int. J. Quantum Chem. 108, 2043–2052 (2008a).
  * Szakács and Surján [2008b] P. Szakács and P. R. Surján, “Iterative solution of bloch-type equations: stability conditions and chaotic behavior,” Journal of mathematical chemistry 43, 314–327 (2008b).
  * Surján, Simon, and Szabados [2022] P. R. Surján, K. Simon, and Á. Szabados, “Stability analysis of the lippmann–schwinger equation,” Molecular Physics , e2091053 (2022).
  * Agarawal, Chakraborty, and Maitra [2020] V. Agarawal, A. Chakraborty, and R. Maitra, “Stability analysis of a double similarity transformed coupled cluster theory,” The Journal of Chemical Physics 153, 084113 (2020).
  * Feigenbaum [1978] M. J. Feigenbaum, “Quantitative universality for a class of nonlinear transformations,” Journal of statistical physics 19, 25–52 (1978).
  * Feigenbaum [1979] M. J. Feigenbaum, “The universal metric properties of nonlinear transformations,” Journal of Statistical Physics 21, 669–706 (1979).
  * Agarawal  _et al._ [2021] V. Agarawal, S. Roy, A. Chakraborty, and R. Maitra, “Accelerating coupled cluster calculations with nonlinear dynamics and supervised machine learning,” The Journal of Chemical Physics 154, 044110 (2021).
  * Agarawal  _et al._ [2022] V. Agarawal, S. Roy, K. K. Shrawankar, M. Ghogale, S. Bharathi, A. Yadav, and R. Maitra, “A hybrid coupled cluster–machine learning algorithm: Development of various regression models and benchmark applications,” The Journal of Chemical Physics 156, 014109 (2022).
  * Haken [1989a] H. Haken, “Synergetics: an overview,” Rep. Prog. Phys. 52, 515–553 (1989a).
  * Haken and Wunderlin [1982] H. Haken and A. Wunderlin, “Slaving principle for stochastic differential equations with additive and multiplicative noise and for discrete noisy maps,” Z. Phys. B 47, 179–187 (1982).
  * Wunderlin [1987] A. Wunderlin, “On the slaving principle,” in  _Lasers and Synergetics_ (Springer, 1987) pp. 140–147.
  * Haken [1983a] H. Haken, “Nonlinear equations. the slaving principle,” in  _Advanced Synergetics: Instability Hierarchies of Self-Organizing Systems and Devices_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 1983) pp. 187–221.
  * Born [1927] M. Born, “Born-oppenheimer approximation,” Ann. Phys 84, 457–484 (1927).
  * Synergetics [1976] H. H. Synergetics, “Introduction and advanced topics,” Springer, Berlin, Heidelberg, New York Tokio 202, 204 (1976).
  * Agarawal, Patra, and Maitra [2021] V. Agarawal, C. Patra, and R. Maitra, “An approximate coupled cluster theory via nonlinear dynamics and synergetics: The adiabatic decoupling conditions,” The Journal of Chemical Physics 155, 124115 (2021).
  * Haken [1975] H. Haken, “Higher order corrections to generalized ginzburg-landau equations of non-equilibrium systems,” Zeitschrift für Physik B Condensed Matter 22, 69–72 (1975).
  * Wunderlin and Haken [1981] A. Wunderlin and H. Haken, “Generalized ginzburg-landau equations, slaving principle and center manifold theorem,” Zeitschrift für Physik B Condensed Matter 44, 135–141 (1981).
  * Strogatz [2018] S. H. Strogatz,  _Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering_ (CRC press, 2018).
  * Alligood, Sauer, and Yorke [1996] K. T. Alligood, T. D. Sauer, and J. A. Yorke, “Two-dimensional maps,” Chaos: An Introduction to Dynamical Systems , 43–104 (1996).
  * Lam [2003] L. Lam,  _Introduction to nonlinear physics_ (Springer Science & Business Media, 2003).
  * Schmidt [1985] G. Schmidt, “Haken, h., advanced synergetics. instabilities of self-organizing systems and devices. berlin-heidelberg-new york-tokyo, springer-verlag 1983. xv. 356 s., 105 abb., dm 98,—. us $38.90. isbn 3-540-12 162-5 (springer series in synergetics 20),” Zeitschrift Angewandte Mathematik und Mechanik 65, 394–394 (1985).
  * Haken [1983b] H. Haken, “Nonlinear equations. the slaving principle,” in  _Advanced Synergetics_ (Springer, 1983) pp. 187–221.
  * Haken [1997] H. Haken, “Discrete dynamics of complex systems,” Discrete Dynamics in Nature and Society 1, 1–8 (1997).
  * Johnson  _et al._ [2006] R. D. Johnson _et al._ , “Nist computational chemistry comparison and benchmark database,” http://srdata. nist. gov/cccbdb (2006).
  * Haken and Portugali [2016] H. Haken and J. Portugali, “Information and self-organization,” (2016).
  * Haken [1985a] H. Haken, “Application of the maximum information entropy principle to selforganizing systems,” Zeitschrift für Physik B Condensed Matter 61, 335–338 (1985a).
  * Haken [1987] H. Haken, “Self-organization and information,” Physica Scripta 35, 247 (1987).
  * Haken [2006] H. Haken,  _Information and self-organization: A macroscopic approach to complex systems_ (Springer Science & Business Media, 2006).
  * Haken and Portugali [2021] H. Haken and J. Portugali, “Information and self-organization ii: Steady state and phase transition,” Entropy 23, 707 (2021).
  * Haken [1985b] H. Haken, “Information, information gain, and efficiency of self-organizing systems close to instability points,” Zeitschrift für Physik B Condensed Matter 61, 329–334 (1985b).
  * Haken [1986] H. Haken, “The maximum entropy principle for non-equilibrium phase transitions: Determination of order parameters, slaved modes, and emerging patterns,” Zeitschrift für Physik B Condensed Matter 63, 487–491 (1986).
  * Haken [1989b] H. Haken, “Synergetics: an overview,” Reports on Progress in Physics 52, 515 (1989b).
  * Schöner and Haken [1986] G. Schöner and H. Haken, “The slaving principle for stratonovich stochastic differential equations,” Zeitschrift für Physik B Condensed Matter 63, 493–504 (1986).
  * Schöner and Haken [1987] G. Schöner and H. Haken, “A systematic elimination procedure for ito stochastic differential equations and the adiabatic approximation,” Zeitschrift für Physik B Condensed Matter 68, 89–103 (1987).
  * Aoki and Kaneko [2013] H. Aoki and K. Kaneko, “Slow stochastic switching by collective chaos of fast elements,” Physical review letters 111, 144102 (2013).
  * Wu and Cao [1990] D.-j. Wu and L. Cao, “Relations between different representations of stochastic adiabatic solutions,” Zeitschrift für Physik B Condensed Matter 81, 451–455 (1990).
  * Constable, McKane, and Rogers [2013] G. W. Constable, A. J. McKane, and T. Rogers, “Stochastic dynamics on slow manifolds,” Journal of Physics A: Mathematical and Theoretical 46, 295002 (2013).
  * Li, Cao, and Wu [1992] J.-H. Li, L. Cao, and D.-J. Wu, “Relations of systematic adiabatic approximation and stochastic adiabatic approximation,” Communications in Theoretical Physics 18, 403 (1992).
  * Wu, He, and Huang [2000] S. Wu, K. He, and Z. Huang, “Suppressing complexity via the slaving principle,” Physical Review E 62, 4417 (2000).
  * Grebogi, Ott, and Yorke [1983] C. Grebogi, E. Ott, and J. A. Yorke, “Crises, sudden changes in chaotic attractors, and transient chaos,” Physica D: Nonlinear Phenomena 7, 181–200 (1983).
  * Takatsuka [2018] K. Takatsuka, “Adiabatic and nonadiabatic dynamics in classical mechanics for coupled fast and slow modes: sudden transition caused by the fast mode against the slaving principle,” Molecular Physics 116, 2556–2570 (2018).