# Integrating Subsystem Embedding Subalgebras and Coupled Cluster Green’s
Function: A Theoretical Foundation for Quantum Embedding in Excitation
Manifold

Bo Peng<EMAIL_ADDRESS>Physical & Computational Science Division, Pacific
Northwest National Laboratory, Richland WA 99352, USA  Karol Kowalski

###### Abstract

In this study, we introduce a novel approach to coupled-cluster Green’s
function (CCGF) embedding by seamlessly integrating conventional CCGF theory
with the state-of-the-art sub-system embedding sub-algebras coupled cluster
(SES-CC) formalism. This integration focuses primarily on delineating the
characteristics of the sub-system and the corresponding segments of the
Green’s function, defined explicitly by active orbitals. Crucially, our work
involves the adaptation of the SES-CC paradigm, addressing the left eigenvalue
problem through a distinct form of Hamiltonian similarity transformation. This
advancement not only facilitates a comprehensive representation of the
interaction between the embedded sub-system and its surrounding environment
but also paves the way for the quantum mechanical description of multiple
embedded domains, particularly by employing the emergent quantum flow
algorithms. Our theoretical underpinnings further set the stage for a
generalization to multiple embedded sub-systems. This expansion holds
significant promise for the exploration and application of non-equilibrium
quantum systems, enhancing the understanding of system-environment
interactions. In doing so, the research underscores the potential of SES-CC
embedding within the realm of quantum computations and multi-scale
simulations, promising a good balance between accuracy and computational
efficiency.

##  I Introduction

Quantum molecular systems situated in complex environments—often referred to
as open quantum molecular systems—engage in dynamic exchanges with their
surroundings. These exchanges, encompassing charge, spin, particle, and
energy, are fundamental processes critical for a diverse range of scientific
fields. These fields include but are not limited to non-equilibrium chemical
transformations, transient spectroscopies, surface physics, and chemistry in
catalysis, photo-induced charge transfer in molecular optoelectronics, and
photo- and thermally-induced spin transitions in molecular spintronics.
Moreover, understanding decoherence in molecular qubits and noise measurements
in quantum thermodynamics is imperative for advancements in quantum
information science. Accurately capturing the dynamical processes of these
systems opens new avenues for the design and characterization of intricate
quantum molecular systems, with far-reaching implications across these
disciplines Breuer  _et al._ (2002); Weiss (2012); May and KüHn (2011);
Wasielewski  _et al._ (2020).

Despite significant enhancements in computational capabilities over recent
decades, the computational expense associated with describing complex
molecular systems accurately remains a considerable challenge. Strategies to
mitigate these computational demands, while preserving the integrity of high-
fidelity descriptions, range from scaling and extrapolation methods to quantum
embedding (QE) techniques Wesolowski  _et al._ (2015); Sun and Chan (2016).
Traditional scaling and extrapolation methods, while valuable, can prove
prohibitively expensive or even infeasible when applied to systems with
extensive degrees of freedom. In contrast, QE techniques have emerged as
particularly potent tools, demonstrating remarkable efficacy in a variety of
real-world scenarios at the cutting edge of chemical research, including areas
like catalysis, biochemistry, and spectroscopy.

The development of an effective QE method generally begins with the
partitioning of the targeted system into smaller, more manageable subsystems.
This often involves defining an embedded region within a broader embedding
environment. The method then strategically combines advanced, yet
computationally intensive, treatments for the embedded region with more
efficient computational approaches for the environmental degrees of freedom.
This synergy aims to minimize the overall computational burden, making the
comprehensive analysis more practical than if one were to attempt a singular,
high-level computation for the entire system. This principle has guided the
development of numerous QE methodologies, each with unique attributes tailored
to specific types of system partitioning and inter-system interactions Jones
_et al._ (2020); Sushko  _et al._ (2018); Libisch  _et al._ (2014); Seijo and
Barandiarán (1999); Jacob and Neugebauer (2014); Wesolowski  _et al._ (2015);
Sun and Chan (2016).

Recognizing the importance of geometrically sensible system partitioning
(e.g., along chemical bonds or between molecular clusters) and accurate
handling of subsystem interactions, especially in scenarios like charge
delocalization and quantum transport, is essential. In these instances,
interactions are so pivotal that quantum processes span regions significantly
larger than the predefined embedded area. Accordingly, the detailed
understanding and treatment of these interactions become paramount in the
development of reliable and predictive models.

Current QE strategies can broadly classify into two main categories: those
combining quantum mechanics with classical molecular mechanics (QM:MM) and
those integrating various levels of quantum mechanical theories (QM:QM)
Warshel and Levitt (1976); Brunk and Rothlisberger (2015); Svensson  _et al._
(1996); Chung  _et al._ (2015). The QM:QM category encompasses several
sophisticated methods, including but not limited to, fully quantum mechanical
treatment through DFT embedding, partition DFT, potential-functional
embedding, embedded mean-field theory (EMFT), Green’s function embedding,
self-energy embedding, density matrix embedding theory (DMET), and stochastic
embedding DFT Jacob and Neugebauer (2014); Govind  _et al._ (1998); Huang  _et
al._ (2011); Manby  _et al._ (2012); Elliott  _et al._ (2010); Nafziger and
Wasserman (2014); Huang and Carter (2011); Fornace  _et al._ (2015);
Inglesfield (1981); Chibani  _et al._ (2016); Lan  _et al._ (2015); Rusakov
_et al._ (2019); Knizia and Chan (2013); Li  _et al._ (2019). These approaches
leverage various representations of quantum mechanics—such as single-particle
density, single-particle Green’s function, and single-particle density
matrix—for performing density functional embedding, Green’s function
embedding, and density matrix embedding, respectively Sun and Chan (2016);
Jones  _et al._ (2020).

Despite the diverse theories and broad application domains, two predominant
challenges arise:

  * •

Constructing a unified embedding formalism that integrates various quantum
mechanical representations—encompassing wave function, reduced density matrix,
Green’s function, self-energy, among others—rigorously, while minimizing the
embedding scheme’s dependence on system partitioning.

  * •

Describing the coupling between the embedded region and the surrounding
environment in a systematically improvable way. This approach should
seamlessly connect the non-interacting subsystem limit (NSL) and the fully
entangled limit (FEL) of the open system, thereby guiding the balance between
accuracy and computational efficiency for various applications.

In response to these challenges, this paper introduces a quantum embedding
methodology that innovatively applies coupled-cluster subsystem embedding sub-
algebra techniques Kowalski (2018) alongside extensive many-body coupled-
cluster representations of wave functions Arponen (1983) and many-body
coupled-cluster Green’s function theory Nooijen and Snijders (1992, 1993,
1995); Peng  _et al._ (2021). Our approach presents a unique method for the
rigorous partitioning of the Green’s function of open quantum systems,
effectively distinguishing between the low and high energy regimes of the
spectrum. This distinction significantly enhances the precision of electronic
structure representations and quantum dynamics descriptions of open systems,
surpassing the constraints found in traditional embedding algorithms.
Moreover, this methodology is inherently adaptable for multiscale/multiphysics
simulations, promising extensive applicability.

The remainder of the paper is structured into four detailed sections. Section
I will succinctly revisit the Sub-system Embedding Sub-algebras (SES) and its
intricate coupled-cluster formalism (SES-CC), highlighting notable properties.
Section II delves into the generalization of SES-CC, employing the
conventional bi-orthogonal CC formalism, focusing on $\Lambda$ parametrization
and its exponential alternatives. In Section III, we dedicate our efforts to
developing the SES-CC Green’s function formalism, illustrating the systematic
decomposition of conventional CCGF calculations into distinct SES-CCGF
computations for the embedded component and its corresponding environment.
Section IV provides a practical demonstration of the SES-CC Green’s function
approach, utilizing a three-site single impurity Anderson model (SIAM) to
elucidate the workflow of SES-CCGF, particularly demonstrating how embedded
subspace calculations can be conducted with precision. It is crucial to note
that the analyses within this study are premised on the exact wave function
limit formulations of CC and CC Green’s function theories. Consequently, any
simplifications or approximations are thereby a reduction from the equations
derived in this exact limit.

##  II Methodology

###  II.1 A brief review of sub-system embedding subalgebras coupled cluster
(SES-CC) formalism and its properties

In the single reference coupled cluster (CC) approach, the electronic ground-
state wave function $|\Psi\rangle$ of an $N$-electron system is assumed to be
parametrized through an exponential anstaz, $|\Psi\rangle=e^{T}|\Phi\rangle$,
where $|\Phi\rangle$ a single reference. In this exponential ansatz, the
exponent operator $T$ (i.e. the cluster operator) comprises excitation
operators that produce excited configurations from the reference. In the
language of second quantization, $T$ can be expressed as

$T=\sum_{k=1}^{m}\frac{1}{(k!)^{2}}\sum_{\begin{subarray}{c}i_{1},\ldots,i_{k};\\\
a_{1},\ldots,a_{k}\end{subarray}}t^{i_{1}\ldots i_{k}}_{a_{1}\ldots
a_{k}}a_{a_{1}}^{\dagger}\cdots a_{a_{k}}^{\dagger}a_{i_{k}}\cdots
a_{i_{1}}~{}.$ (1)

Here, the indices $i_{1},i_{2},\ldots$ ($a_{1},a_{2},\ldots$) denote occupied
(unoccupied) spin-orbitals in the reference $|\Phi\rangle$, $a_{p}^{\dagger}$
($a_{p}$) is the standard creation (annihilation) operator that add (remove)
electron at spin-orbital $p$, and $t^{i_{1}\ldots i_{k}}_{a_{1}\ldots
a_{k}}$’s denote cluster amplitudes, $m$ is the excitation level ($\leq N$)
that defines the approximate approach within the CC framework (e.g. $m=2$ for
CCSD, Purvis and Bartlett (1982), $m=3$ for CCSDT, Noga and Bartlett (1987,
1988); Scuseria and Schaefer (1988) $m=4$ for CCSDTQ ($m=4$), Kucharski and
Bartlett (1991); Oliphant and Adamowicz (1991)). Once $T$ is defined, we can
plug the exponential ansatz into the Schrödinger equation to get the energy-
dependent form of the CC equation

$\displaystyle(P^{(N)}+Q^{(N)})He^{T}|\Phi\rangle$
$\displaystyle=E_{\text{CC}}^{(N)}(P^{(N)}+Q^{(N)})e^{T}|\Phi\rangle,$ (2)

where $H$ is the electronic Hamiltonian, $E_{\rm CC}^{(N)}$ is the coupled
cluster energy, and $P^{(N)}$ and $Q^{(N)}$ are operators that project Eq. (2)
onto the reference and excited configurations of the $N$-electron system,
respectively,

$\displaystyle P^{(N)}$ $\displaystyle=|\Phi\rangle\langle\Phi|,$
$\displaystyle Q^{(N)}$
$\displaystyle=\sum_{k=1}^{m}~{}\sum_{\begin{subarray}{c}i_{1}<\ldots<i_{k};\\\
a_{1}<\ldots<a_{k}\end{subarray}}|\Phi_{i_{1}\ldots i_{k}}^{a_{1}\ldots
a_{k}}\rangle\langle\Phi_{i_{1}\ldots i_{k}}^{a_{1}\ldots a_{k}}|.$ (3)

Focusing on the excitation operators $T$ in Eq. (1), as discussed in Ref. 29,
these operators constitute a commutative algebra, $\mathfrak{g}^{(N)}$, (see
Fig. 1). Also, we can define a class of active space excitation/de-excitation
subalgebras, denoted as
$\mathfrak{h}\coloneqq\mathfrak{g}^{(N)}(\mathcal{R},\mathcal{S})$, which
facilitate the excitation of electrons from an active subspace $\mathcal{R}$
within the occupied orbital space to an active subspace $\mathcal{S}$ within
the virtual orbital space, encompassing all possible excitations. In the
context of the particle-hole formalism, all subalgebras $\mathfrak{h}$ and the
overarching excitation algebra $\mathfrak{g}^{(N)}$ are commutative. In
particular scenarios where the subalgebras include either all occupied and
selected virtual orbitals ($\mathcal{S}$) or a specific set of occupied
orbitals ($\mathcal{R}$) and all virtual orbitals, we adopt the concise
notations $\mathfrak{g}^{(N)}(\mathcal{S})$ and
$\mathfrak{g}^{(N)}(\mathcal{R})$, respectively. In some instances, the
notation $\mathfrak{g}^{(N)}(x_{\mathcal{R}},y_{\mathcal{S}})$ is employed to
specify the number of orbitals, $x$ and $y$, included in the active subspaces
$\mathcal{R}$ and $\mathcal{S}$, respectively.

The purpose of introducing the active-space excitation/de-excitation
subalgebras lies in the fact that the cluster operator $T$ can be generally
decomposed to the commutative components that belong to $\mathfrak{h}$
(denoted as the internal part, $T_{\rm int}(\mathfrak{h})$) and
$\mathfrak{g}^{(N)}-\mathfrak{h}$ (denoted as an external part, $T_{\rm
ext}(\mathfrak{h})$), respectively, i.e.

$T=T_{\rm int}(\mathfrak{h})+T_{\rm
ext}(\mathfrak{h}),~{}~{}\text{with}~{}~{}[T_{\rm int}(\mathfrak{h}),T_{\rm
ext}(\mathfrak{h})]=0.$ (4)

Accordingly, the corresponding wave function can be rewritten as

$\displaystyle|\Psi\rangle$ $\displaystyle=e^{T}|\Phi\rangle=e^{T_{\rm
ext}+T_{\rm int}}|\Phi\rangle=e^{T_{\rm ext}}|\Psi_{\rm int}\rangle.$ (5)

where we denote $|\Psi_{\rm int}\rangle=e^{T_{\rm int}}|\Phi\rangle$ as the
auxiliary CAS-type wave function corresponding to the internal active
subspace, and drop the explicit $\mathfrak{h}$-dependence in $T_{\rm int}$,
$T_{\rm ext}$, and $|\Psi_{\rm int}\rangle$ for the simplification of notation
(and same below).

The subalgebras $\mathfrak{h}$ can be further refined into one class that
uniquely maps to the CC approximations if (a) $|\Psi(\mathfrak{h})\rangle$
preserves the symmetry properties of $|\Psi\rangle$ and $|\Phi\rangle$, and
(b) the $e^{T_{\rm int}}|\Phi\rangle$ ansatz generates the full configuration
interaction (FCI) expansion for the sub-system defined by the active subspace
corresponding to the $\mathfrak{h}$ subalgebra. In Ref. 29 this is referred as
sub-system embedding subalgebras (SESs).

Figure 1: All the excitations in an $N$-electron system based on a reference
$|\Phi\rangle$ form an excitation algebra $\mathfrak{g}^{N}$ that includes
excitation subalgebras $\mathfrak{g}^{N}(\mathcal{R},\mathcal{S})$’s.

The SES $\mathfrak{h}$ can be employed to derive an alternative representation
of the CC equation that is naturally adapted to the sub-system embedding.
Following the decomposition of $T$, we can further decompose the projection
operator $Q^{(N)}$ as

$\displaystyle Q^{(N)}=Q_{\rm int}^{(N)}+Q_{\rm ext}^{(N)},$ (6)

where $Q_{\rm int}^{(N)}$ and $Q_{\rm ext}^{(N)}$ are constructed from the
excited configurations generated by respectively acting $T_{\rm int}$ and
$T_{\rm ext}$ on the reference. Denote

$\displaystyle\mathcal{P}_{\rm int}^{(N)}=P^{(N)}+Q_{\rm int}^{(N)}$ (7)

(such that $\mathcal{P}_{\rm int}^{(N)}+Q_{\rm ext}^{(N)}=\mathbf{1}^{(N)}$)
and define $\overline{H}_{\rm ext}=e^{-T_{\rm ext}}He^{T_{\rm ext}}$, we can
then re-write CC energy equation as

$\displaystyle\mathcal{P}_{\rm int}^{(N)}\overline{H}_{\rm ext}|\Psi_{\rm
int}\rangle$ $\displaystyle=E_{\rm CC}^{(N)}\mathcal{P}_{\rm
int}^{(N)}|\Psi_{\rm int}\rangle,$ (8)

for the SES $\mathfrak{h}$ with $T_{\rm ext}$ being obtained from the
projections of $(He^{T})_{\rm C}|\Phi\rangle$ on the excited configurations
that are not included in the SES $\mathfrak{h}$

$\displaystyle Q_{\rm ext}^{(N)}\overline{H}_{N}|\Phi\rangle$
$\displaystyle=0,$ (9)

where the normal ordered Hamiltonian is defined as

$\displaystyle\overline{H}_{N}=e^{-T}He^{T}-E_{\rm
CC}^{(N)}=\overline{H}-E_{\rm CC}^{(N)}.$ (10)

Because $Q_{\rm ext}^{\rm(N)}|\Psi_{\rm int}\rangle=0$, the identity between
$\overline{H}_{\rm ext}$ and $|\Psi_{\rm int}\rangle$ in Eq. (8) can be
equivalently replaced by $\mathcal{P}_{\rm int}^{(N)}$, and Eq. (8) can be
equivalently rewritten as

$\displaystyle\mathcal{P}_{\rm int}^{(N)}\overline{H}_{\rm
ext}\mathcal{P}_{\rm int}^{(N)}|\Psi_{\rm int}\rangle$ $\displaystyle=E_{\rm
CC}^{(N)}\mathcal{P}_{\rm int}^{(N)}|\Psi_{\rm int}\rangle.$ (11)

with $\mathcal{P}_{\rm int}^{(N)}|\Psi_{\rm int}\rangle$ denoting the CAS-type
wave function projected to the active space. As can be seen from Eq. (11),
$E_{\rm CC}^{(N)}$ and $|\Psi_{\rm int}\rangle$ (and hence $T_{\rm int}$) can
now be obtained by directly diagonalizing $\overline{H}_{\rm ext}$ in the
$(P^{(N)}+Q_{\rm int}^{(N)})$ space. The generated effective Hamiltonian,

$\displaystyle\overline{H}_{\rm ext}^{\rm eff}=\mathcal{P}_{\rm
int}^{(N)}\overline{H}_{\rm ext}\mathcal{P}_{\rm int}^{(N)},$ (12)

is $\mathfrak{h}$-dependent, and provides a unique and formally exact way to
downfold the full electronic Hamiltonian to the active space associated with
the SES $\mathfrak{h}$. Due to this feature, SES has been proposed to work
with double or generalized unitary CC formalism to provide effective Hermitian
Hamiltonians in the variational quantum approaches, as reported in recent
studies Bauman  _et al._ (2019); Metcalf  _et al._ (2020); Bauman  _et al._
(2020), showing the promising ability of capturing a large part of dynamical
electron correlation usually missing in relatively smaller active space
quantum simulations.

###  II.2 Extension of SES formalism for $\Lambda$ formalism

In contrast to the $T$ operator, the de-excitation $\Lambda$ operator involves
all the de-excitation operators $a_{i_{1}}^{\dagger}\cdots
a_{i_{k}}^{\dagger}a_{a_{k}}\cdots a_{a_{1}}$ defined with respect to the
$|\Phi\rangle$ reference. Several CC formulations such as linear response CC,
analytical CC gradients, and CC Green’s function formulations employ the so-
called bi-orthogonal CC formalism, Arponen (1983); Szalay  _et al._ (1995);
Piecuch and Bartlett (1999); Salter  _et al._ (1987); Stanton and Bartlett
(1993); Monkhorst (1977); Koch and Jørgensen (1990); Nooijen and Snijders
(1992, 1993, 1995); Kowalski  _et al._ (2014); Bhaskaran-Nair  _et al._
(2016); Peng and Kowalski (2016, 2018); Helgaker  _et al._ (2014); Schirmer
and Mertins (2010) which utilizes different parametrization for the bra
($\langle\Psi|$) of the ground-state wave functions

$\displaystyle\langle\Psi|$
$\displaystyle=\langle\Phi|(1+\Lambda)e^{-T}~{}~{},$ (13)

where it is easy to verify that $\langle\Psi|\Psi\rangle=1$ and the de-
excitation $\Lambda$ operator can be determined from solving the left-
eigenvalue problem for the similarity transfromed Hamiltonian $\overline{H}$

$\langle\Phi|(1+\Lambda)\overline{H}(P^{(N)}+Q^{(N)})=E_{\rm
CC}^{(N)}\langle\Phi|(1+\Lambda)(P^{(N)}+Q^{(N)})~{}.$ (14)

In analogy to the “right” eigenvalue problem discussed in the previous
Section, i.e. we can decompose for a given de-excitation SES $\mathfrak{h}$
the exact $\Lambda$ operator from Eq. (14) in the same way as the $T$ operator
(see Eq. (4) and drop off the explicit $\mathfrak{h}$-dependence for the
simplification of notation),

$\Lambda=\Lambda_{\rm int}+\Lambda_{\rm
ext},~{}~{}\text{with}~{}~{}[\Lambda_{\rm int},\Lambda_{\rm ext}]=0.$ (15)

In the exact wave function expansion, the left eigenvector of $\overline{H}$
can be alternatively represented in the extended CC form Arponen (1983),

$\langle\Phi|(1+\Lambda)=\langle\Phi|e^{S}~{},$ (16)

where the de-excitation cluster operator $S$ can also be decomposed into
commuting active and external parts

$S=S_{\rm int}+S_{\rm ext},~{}~{}\text{with}~{}~{}[S_{\rm int},S_{\rm
ext}]=0.$ (17)

In the exact wave function limit both parametrizations are equivalent

$\langle\Phi|(1+\Lambda_{\rm int}+\Lambda_{\rm ext})=\langle\Phi|e^{S_{\rm
int}}e^{S_{\rm ext}}~{}.$ (18)

After projecting onto $Q_{\rm int}^{(N)}$ and $Q_{\rm ext}^{(N)}$ subspaces
one obtains the following identities

$\displaystyle\langle\Phi|(1+\Lambda_{\rm int})\mathcal{P}_{\rm int}^{(N)}$
$\displaystyle=\langle\Phi|e^{S_{\rm int}}\mathcal{P}_{\rm int}^{(N)}~{},$
(19) $\displaystyle\langle\Phi|\Lambda_{\rm ext}Q_{\rm ext}^{(N)}$
$\displaystyle=\langle\Phi|e^{S_{\rm int}}e^{S_{\rm ext}}Q_{\rm
ext}^{(N)}~{}.$ (20)

The cluster analysis of Eqs. (19,20) gives

$\displaystyle\Lambda_{{\rm int},1}$ $\displaystyle=S_{{\rm int},1}~{},$
$\displaystyle\Lambda_{{\rm int},2}$ $\displaystyle=S_{{\rm
int},2}+\frac{1}{2}S_{{\rm int},1}^{2}~{},$ $\displaystyle\Lambda_{{\rm
int},3}$ $\displaystyle=S_{{\rm int},3}+S_{{\rm int},1}S_{{\rm
int},2}+\frac{1}{6}S_{{\rm int},1}^{3}~{},$ (21) $\displaystyle\vdots$

and

$\displaystyle\Lambda_{{\rm ext},1}$ $\displaystyle=S_{{\rm ext},1}~{},$
$\displaystyle\Lambda_{{\rm ext},2}$ $\displaystyle=S_{{\rm
ext},2}+\frac{1}{2}S_{{\rm ext},1}^{2}+S_{{\rm int},1}S_{{\rm ext},1}~{},$
$\displaystyle\Lambda_{{\rm ext},3}$ $\displaystyle=S_{{\rm ext},3}+S_{{\rm
ext},1}S_{{\rm ext},2}+S_{{\rm int},1}S_{{\rm ext},2}$
$\displaystyle~{}~{}+S_{{\rm int},2}S_{{\rm ext},1}+\frac{1}{2}S_{{\rm
int},1}^{2}S_{{\rm ext},1}$ $\displaystyle~{}~{}+\frac{1}{2}S_{{\rm
int},1}S_{{\rm ext},1}^{2}+\frac{1}{6}S_{{\rm ext},1}^{3}~{},$ (22)
$\displaystyle\vdots$

Reverse mapping can be readily derived from the above equations. Introducing
exponential parametrizations (16-20) into Schrödinger equation, one obtains

$\langle\Phi|e^{S_{\rm int}}e^{S_{\rm ext}}\overline{H}e^{-S_{\rm ext}}=E_{\rm
CC}^{(N)}\langle\Phi|e^{S_{\rm int}}~{}.$ (23)

Projecting these equations onto $\mathcal{P}_{\rm int}^{(N)}$ and $Q_{\rm
ext}^{(N)}$ one gets the working equations for $e^{S_{\rm int}}$ and $S_{\rm
ext}$

$\displaystyle\langle\Phi|e^{S_{\rm int}}\overline{\overline{H}}_{\rm
ext}\mathcal{P}_{\rm int}^{(N)}$ $\displaystyle=E_{\rm
CC}^{(N)}\langle\Phi|e^{S_{\rm int}}\mathcal{P}_{\rm int}^{(N)}~{},$ (24)
$\displaystyle\langle\Phi|(1+\Lambda_{\rm int})\overline{\overline{H}}_{\rm
ext}Q_{\rm ext}^{(N)}$ $\displaystyle=0~{},$ (25)

where

$\overline{\overline{H}}_{\rm ext}=e^{S_{\rm ext}}\overline{H}e^{-S_{\rm
ext}}=(e^{S_{\rm ext}}\overline{H})_{C}~{}.$ (26)

One can observe that, similar to Eq. (12), the energy $E_{\rm CC}^{(N)}$ and

$\displaystyle\langle\Psi_{\rm int}|=\langle\Phi|e^{S_{\rm
int}}\mathcal{P}_{\rm int}^{(N)}$ (27)

(and therefore $S_{\rm int}$ and $\Lambda_{\rm int}$) can be obtained by
diagonalizing another effective many-body Hamiltonian

$\displaystyle\overline{\overline{H}}_{\rm ext}^{\rm eff}=\mathcal{P}_{\rm
int}^{(N)}\overline{\overline{H}}_{\rm ext}\mathcal{P}_{\rm int}^{(N)}~{}.$
(28)

Again, the amplitudes defining $\overline{\overline{H}}_{\rm ext}$ and
corresponding to $T$ and $S_{\rm ext}$ operators are decoupled from $S_{\rm
int}$ or $\Lambda_{\rm int}$ amplitudes.

Note that from the above construction $\overline{\overline{H}}_{\rm ext}^{\rm
eff}$ is generally not equal to $\overline{H}_{\rm ext}^{\rm eff}$, thus (i)
$\langle\Psi_{\rm int}|\Psi_{\rm int}\rangle\neq 1$ and (ii)
$\overline{\overline{H}}_{\rm ext}^{\rm eff}$ and $\overline{H}_{\rm ext}^{\rm
eff}$ need to be constructed separately to solve both $T_{\rm int}$ and
$\Lambda_{\rm int}$ (or $S_{\rm int}$). To address these issues, we can
construct the third form of effective Hamiltonian based on the $\Lambda$-CC
energy form (int he exact wave function limit)

$\displaystyle E_{\rm CC}^{(N)}$
$\displaystyle=\langle\Phi|(1+\Lambda)\overline{H}|\Phi\rangle$
$\displaystyle=\langle\Phi|e^{S_{\rm int}}e^{S_{\rm ext}}e^{-T_{\rm
int}}\overline{H}_{\rm ext}e^{T_{\rm int}}|\Phi\rangle$ (29)

Inserting resolution of identity, $1=e^{-T_{\rm int}}e^{T_{\rm int}}$, between
$e^{S_{\rm int}}$ and $e^{S_{\rm ext}}$, and define the effective Hamiltonian

$\displaystyle\widetilde{H}^{\rm eff}_{\rm ext}$
$\displaystyle=\mathcal{P}_{\rm int}^{(N)}W\overline{H}_{\rm
ext}\mathcal{P}_{\rm int}^{(N)},$ (30) $\displaystyle W$
$\displaystyle=e^{T_{\rm int}}e^{S_{\rm ext}}e^{-T_{\rm int}},$ (31)

one gets another working equation for $S_{\rm int}$ and $T_{\rm int}$

$\displaystyle E_{\rm CC}^{(N)}=\langle\Phi|e^{S_{\rm int}}e^{-T_{\rm
int}}\widetilde{H}^{\rm eff}_{\rm ext}e^{T_{\rm int}}|\Phi\rangle,$ (32)

where the left and right eigenvectors of $\widetilde{H}^{\rm eff}_{\rm ext}$
are given by

$\displaystyle\langle\Psi_{\rm int}|=\langle\Phi|e^{S_{\rm int}}e^{-T_{\rm
int}},~{}~{}|\Psi_{\rm int}\rangle=e^{T_{\rm int}}|\Phi\rangle.$ (33)

and $\langle\Psi_{\rm int}|\Psi_{\rm int}\rangle=1$. The construction and use
of the above three effective Hamiltonians are summarized in Tab. 1

Effective Hamiltonian | $E_{\rm CC}^{(N)}$ | $T_{\rm int}$ | $S_{\rm int}$
---|---|---|---
$\overline{H}_{\rm ext}^{\rm eff}=\mathcal{P}_{\rm int}^{(N)}\overline{H}_{\rm ext}\mathcal{P}_{\rm int}^{(N)}$ | $\surd$ | $|\Psi_{\rm int}\rangle=\mathcal{P}_{\rm int}^{(N)}e^{T_{\rm int}}|\Phi\rangle$ | $\bigtimes$
$\overline{\overline{H}}_{\rm ext}^{\rm eff}=\mathcal{P}_{\rm int}^{(N)}\overline{\overline{H}}_{\rm ext}\mathcal{P}_{\rm int}^{(N)}$ | $\surd$ | $\bigtimes$ | $\langle\Psi_{\rm int}|=\langle\Phi|e^{S_{\rm int}}\mathcal{P}_{\rm int}^{(N)}$
$\widetilde{H}^{\rm eff}_{\rm ext}=\mathcal{P}_{\rm int}^{(N)}W\overline{H}_{\rm ext}\mathcal{P}_{\rm int}^{(N)}$ | $\surd$ | $|\Psi_{\rm int}\rangle=\mathcal{P}_{\rm int}^{(N)}e^{T_{\rm int}}|\Phi\rangle$ | $\langle\Psi_{\rm int}|=\langle\Phi|e^{S_{\rm int}}e^{-T_{\rm int}}\mathcal{P}_{\rm int}^{(N)}$
Table 1: The construction of three SES-CC effective Hamiltonians and their use
to get $T_{\rm int}$ and $S_{\rm int}$.

###  II.3 Green’s function embedding based on the sub-system excitation/de-
excitation sub-algebras

In this section we describe a procedure for quantum description of the
composite system based on the integration quantum region (or embedded region)
with the environment using Green’s function formalism. Our strategy is valid
for sub-systems defined by appropriate active spaces. In the spirit of quantum
flow formulation Kowalski (2021), the SES-CC formalism assures a number of
important properties which are ideally suited for multi-center quantum
embedding.

Let’s assume that all active orbitals
$\\{\phi_{M}\\}\in(\mathcal{R},\mathcal{S})$ defining
$g^{(N)}(\mathcal{R},\mathcal{S})$ and
$g^{(N)}(\mathcal{R},\mathcal{S})^{\dagger}$ that determine $T_{\rm int}$ and
$\Lambda_{\rm int}$ operators are defined such that in the non-interacting
limit (NLS) when composite system AB is separated into non-interacting sub-
systems A and B all
$g^{(N)}(\mathcal{R},\mathcal{S})$/$g^{(N)}(\mathcal{R},\mathcal{S})^{\dagger}$
orbitals evolve to orbitals localized on the sub-subsystem A while the
remaining orbitals evolve into orbitals localized on sub-system B.

In this section we will analyze the asymptotic limit of the sub-block
$G_{KL}(\omega)$ of entire CC Green’s function matrix corresponding to
spinorbitals $K$ and $L$ defining
$g^{(N)}(\mathcal{R},\mathcal{S})$/$g^{(N)}(\mathcal{R},\mathcal{S})^{\dagger}$
sub-algebras and which in the NLS evolve into $K_{A}$ and $L_{A}$ spinorbitals
localized on sub-system A, i.e.,

$K\xrightarrow{NSL}K_{A}~{}~{},L\xrightarrow{NSL}L_{A}~{}.$ (34)

Without loss of generality we will focus on the $\omega$-dependent CC Green’s
function $G_{KL}(\omega)$,

$\displaystyle G_{KL}(\omega)=$
$\displaystyle\langle\Phi|(1+\Lambda)\overline{a_{L}^{\dagger}}X_{K}(\omega)|\Phi\rangle+\langle\Phi|(1+\Lambda)\overline{a_{K}}Y_{L}(\omega)|\Phi\rangle~{},$
(35)

where $\overline{a_{K}}$ and $\overline{a_{L}^{\dagger}}$ operators are
defined as Bhaskaran-Nair  _et al._ (2016); Peng and Kowalski (2016)

$\displaystyle\overline{a_{K}}$
$\displaystyle=e^{-T}a_{K}e^{T}=a_{K}+[a_{K},T],$
$\displaystyle\overline{a_{L}^{\dagger}}$
$\displaystyle=e^{-T}a_{L}^{\dagger}e^{T}=a_{L}^{\dagger}+[a_{L}^{\dagger},T]~{},$
(36)

and $X_{K}(\omega)$ and $Y_{L}(\omega)$ satisfy the equations

$\displaystyle
Q^{(N-1)}(\omega_{-i0}+\overline{H}_{N})X_{K}(\omega)|\Phi\rangle$
$\displaystyle=Q^{(N-1)}\overline{a_{K}}|\Phi\rangle~{},$ (37) $\displaystyle
Q^{(N+1)}(\omega_{+i0}-\overline{H}_{N})Y_{L}(\omega)|\Phi\rangle$
$\displaystyle=Q^{(N+1)}\overline{a^{\dagger}_{L}}|\Phi\rangle~{},$ (38)

First, using decomposition (18,19), i.e.,

$\langle\Phi|(1+\Lambda)=\langle\Phi|(1+\Lambda_{\rm int})e^{S_{\rm ext}}~{},$
(39)

we can represent sub-block matrix elements $G_{KL}(\omega)$ as

$\displaystyle G_{KL}(\omega)$ $\displaystyle=\langle\Phi|(1+\Lambda_{\rm
int})e^{S_{\rm ext}}\overline{a_{L}^{\dagger}}X_{K}(\omega)|\Phi\rangle$
$\displaystyle+\langle\Phi|(1+\Lambda_{\rm int})e^{S_{\rm
ext}}\overline{a_{K}}Y_{L}(\omega)|\Phi\rangle~{},$ (40)

Next, if we formally decompose the $X_{K}(\omega)$ and $Y_{L}(\omega)$
operators into their internal and external parts

$\displaystyle X_{K}(\omega)$ $\displaystyle=X_{K,{\rm int}}(\omega)+X_{K,{\rm
ext}}(\omega)~{},$ (41) $\displaystyle Y_{L}(\omega)$ $\displaystyle=Y_{L,{\rm
int}}(\omega)+Y_{L,{\rm ext}}(\omega)~{},$ (42)

with the internal parts, $X_{K,{\rm int}}(\omega)$ and $Y_{L,{\rm
int}}(\omega)$, reconstruct the correlation effects in the embedded sub-system
(in $(N-1)$\- and $(N+1)$-electron Hilbert spaces, respectively) described by
active orbitals whereas the external parts, $X_{K,{\rm ext}}(\omega)$ and
$Y_{L,{\rm ext}}(\omega)$, describe correlation effects in the environment and
effects that simultaneously correlate embedded system and the environment, we
can then decompose the $G_{KL}(\omega)$ as follows

$\displaystyle G_{KL}(\omega)$ $\displaystyle=G^{\rm int}_{KL}(\omega)+G^{\rm
ext}_{KL}(\omega)~{},$ (43) $\displaystyle G^{\rm int}_{KL}(\omega)$
$\displaystyle=\langle\Phi|(1+\Lambda_{\rm
int})\overline{a_{L}^{\dagger}}X_{K,{\rm
int}}|\Phi\rangle+\langle\Phi|(1+\Lambda_{\rm int})\overline{a_{K}}Y_{L,{\rm
int}}|\Phi\rangle~{},$ (44) $\displaystyle G^{\rm ext}_{KL}(\omega)$
$\displaystyle=\langle\Phi|(1+\Lambda_{\rm int})(e^{S_{\rm
ext}}-1)\overline{a_{L}^{\dagger}}(X_{K,{\rm int}}+X_{K,{\rm
ext}})|\Phi\rangle+\langle\Phi|(1+\Lambda_{\rm int})(e^{S_{\rm
ext}}-1)\overline{a_{K}}(Y_{L,{\rm int}}+Y_{L,{\rm ext}})|\Phi\rangle~{}.$
(45)

As shown in Appendices A, B, C which analyze the asymptotic behavior of $T$,
$\Lambda_{\rm int}$, $S_{\rm ext}$, $X_{K}(\omega)$, and $Y_{K}(\omega)$
operators, one can demonstrate (see Appendix C) that $G^{\rm
int}_{KL}(\omega)$ evolves in the NSL to the CC Green’s function for isolated
subsystem A. Therefore, we identify the $G^{\rm int}_{KL}(\omega)$ part
describing embedded subsystem A while the $G^{\rm ext}_{KL}(\omega)$ part
describes the interactions of the embedded system with the surrounding
environment.

Our embedding procedure is based on the construction of the following
representation of the Green’s function,

$\displaystyle\mathbf{G}(\omega)=\left[\begin{array}[]{c|c}\mathbf{G}^{\rm
emb}(\omega)&\mathbf{G}^{\rm env,1}(\omega)\\\ \hline\cr\mathbf{G}^{\rm
env,2}(\omega)&\mathbf{G}^{\rm env,3}(\omega)\end{array}\right],$ (48)

where the embedded block $\mathbf{G}^{\rm emb}(\omega)$ and the environment
blocks ($\mathbf{G}^{\rm env,1}(\omega)$, $\mathbf{G}^{\rm env,2}(\omega)$,
and $\mathbf{G}^{\rm env,3}(\omega)$) are defined by matrix elements
$G_{KL}(\omega)$, $G_{K\overline{p}}(\omega)$, $G_{\overline{p}K}(\omega)$,
and $G_{\overline{p}\overline{q}}(\omega)$, respectively. Here $\overline{p}$
and $\overline{q}$ designate inactive spinorbitals.

For several active-space problems such as the quantum flow equations discussed
for CC formalism in Ref.56, schematic representation of the Green’s function
featuring multiple embedded parts ($\mathbf{G}^{\rm emb(n)}(\omega)$,
$n=1,\cdots,M$, $M\geq 2$ ) is shown below:

$\displaystyle\mathbf{G}(\omega)=\left[\begin{array}[]{c|c|c}\mathbf{G}^{\rm
emb(1)}(\omega)&\mathbf{G}^{\rm env}(\omega)&\cdots\\\
\hline\cr\mathbf{G}^{\rm env}(\omega)&\mathbf{G}^{\rm
emb(2)}(\omega)&\cdots\\\ \hline\cr\vdots&\vdots&\ddots\end{array}\right],$
(52)

The properties discussed in Appendices A, B, C guarantee that in NSL the above
matrix assumes block-form with blocks defined for embedding regions and the
environment, where $X_{K{\rm(i),int(i)}}(\omega)$ and
$X_{K{\rm(i),ext(i)}}(\omega)$, as well as $Y_{L{\rm(i),int(i)}}(\omega)$ and
$Y_{L{\rm(i),ext(i)}}(\omega)$, are defined through

$\displaystyle Q_{\rm
int(i)}^{(N-1)}[\omega_{-i0}+\overline{H}_{N}](X_{K{\rm(i),int(i)}}(\omega)+X_{K{\rm(i),ext(i)}}(\omega))|\Phi\rangle$
$\displaystyle=Q_{\rm
int(i)}^{(N-1)}\overline{a_{K(i)}}|\Phi\rangle~{},(i=1,\ldots,M)$ (53)
$\displaystyle Q_{\rm
int(i)}^{(N+1)}[\omega_{+i0}-\overline{H}_{N}](Y_{L{\rm(i),int(i)}}(\omega)+Y_{L{\rm(i),ext(i)}}(\omega))|\Phi\rangle$
$\displaystyle=Q_{\rm
int(i)}^{(N+1)}\overline{a^{\dagger}_{L(i)}}|\Phi\rangle~{},(i=1,\ldots,M)$
(54)

where the index $i$ goes over all active centers (active spaces). The flow-
type representation of the CC Green’s function embedding in (48) and (52)
offers a flexible selection mechanism for cluster amplitudes included in the
$N$-electron problem. As discussed in Ref. 56 the cluster operator $T$ (the
$\Lambda$ operator is defined in a similar way) is defined as a combination of
unique cluster amplitudes defining the set of internal cluster operators for
each active space involved in the flow. The choice of the active spaces may be
dictated by the nature of processes of interest. For example, if one is
interested in describing transport processes, the active spaces should be
chosen to cover the spatial area relevant to the junction. For other
applications, such as calculating core-level binding energies, the active
spaces in the flow should be chosen to provide a way of correlating core
orbitals with valence orbitals, which is also important in bringing a balance
between ground-state and core-level states. The flow-type representation of
the CC Green’s function also offers a possibility of representing the entire
problem in the form of small-dimensionality computational blocks. Remarkably,
for various active space SES-CC problem either $\overline{H}_{\rm ext}^{\rm
eff}$ or $\overline{\overline{H}}_{\rm ext}^{\rm eff}$ is providing form of
the effective interactions that can be used to calculate the exact energy of
the $N$-electron system (as long as $T_{\rm ext}$ is known) in the active
space. In this context, take $\overline{H}_{\rm ext}^{\rm eff}$ for an
example, by diagonalizing this effective Hamiltonian, we get $E_{\rm
CC}^{(N)}$ and its bi-orthonormal left and right eigenvectors,
$\langle\Psi_{\rm int}|$ and $|\Psi_{\rm int}\rangle$, that bypass the
$\Lambda$-parameterization but satisfy $\langle\Psi_{\rm int}|\Psi_{\rm
int}\rangle=1$. We could then reformulate the sub-block (40) as

$\displaystyle G_{KL}$ $\displaystyle=\langle\Psi_{\rm int}|(P^{(N)}+Q_{\rm
int}^{(N)})W\overline{a_{L}^{\dagger}}_{,\rm ext}Q_{\rm
int}^{(N-1)}(\omega_{-i0^{+}}+\overline{H}_{N,\rm ext})^{-1}Q_{\rm
int}^{(N-1)}\overline{a_{K}}_{,\rm ext}(P^{(N)}+Q_{\rm int}^{(N)})|\Psi_{\rm
int}\rangle$ $\displaystyle+\langle\Psi_{\rm int}|(P^{(N)}+Q_{\rm
int}^{(N)})W\overline{a_{K}}_{,\rm ext}Q_{\rm
int}^{(N+1)}(\omega_{+i0^{+}}-\overline{H}_{N,\rm ext})^{-1}Q_{\rm
int}^{(N+1)}\overline{a_{L}^{\dagger}}_{,\rm ext}(P^{(N)}+Q_{\rm
int}^{(N)})|\Psi_{\rm int}\rangle,$ (55)

where $\overline{H}_{N,\rm ext}=\overline{H}_{\rm
ext}-\langle\overline{H}_{\rm ext}\rangle$ (in this paper we assume for
certain systems, the $\langle\overline{H}_{\rm ext}\rangle$ can be
approximated by $E_{\rm CC}^{(N)}$), and

$\displaystyle\overline{a_{K}}_{,\rm ext}$ $\displaystyle=e^{-T_{\rm
ext}}a_{K}e^{T_{\rm ext}}~{},~{}~{}\overline{a_{L}^{\dagger}}_{,\rm
ext}=e^{-T_{\rm ext}}a_{L}^{\dagger}e^{T_{\rm ext}}~{}.$ (56)

Note that (i) in (55) $\langle\Psi_{\rm int}|$ essentially approximates
$\langle\Phi|(1+\Lambda)e^{-T_{\rm int}}$, such that the projection
$(P^{(N)}+Q_{\rm int}^{(N)})$ can be applied between $\langle\Psi_{\rm int}|$
and $\overline{a_{K}}_{,\rm ext}$ (or $\overline{a_{L}^{\dagger}}_{,\rm
ext}$), (ii) the pole positions associated with the effectiven Hamiltonian are
exactly part of the spectrum of similarity-transformed Hamiltonian that is
associated with the active space, and (iii) the similar reformulation can also
be obtained for the unitary version of CC Ansatz. Note that the unitary CC
effective Hamiltonians are expressed in terms of anti-Hermitian cluster
operators which involve non-terminating expansions in exact limit, and thus
require approximate forms in practice. In the context of double unitary CC
(DUCC) based effective Hamiltonian Bauman  _et al._ (2019); Kowalski and
Bauman (2020); Kowalski (2021); Bauman and Kowalski (2021), approximate
effective Hamiltonian forms involving single- and double-commutator expansions
have been calibrated for reproducing exact ground state energy using small
active spaces, and thus showing potential for working with Green’s function
embedding approaches. Relevant discussions of this topic will be presented in
our future work.

##  III Numerical results and discussion

###  III.1 SES-CCGF workflow

In practical applications approximate schemes of SES-CCGF formalism can be
driven by the truncation of excitation levels used to define these operators
or by using low-order perturbative approximations. In this section, we
demonstrate the SES-CCGF workflow by computing the spectral function of a
single-impurity Anderson model (SIAM) with three electrons and two bath site,
$N_{\rm bath}=2$. The SIAM is a well-established model that describes a system
consisting of a localized impurity with a repulsive interaction between
opposite spin electrons coupled with a non-interacting electronic bath.
Anderson (1961); Newns (1969); Wilson (1975) This model has been successfully
employed to explain the metal-insulator transition through Anderson
localization. Here, we apply SES-CCGF with single and double excitations to
the SIAM to demonstrate how SES-CCSDGF works and how to apply further
approximation.

Figure 2: Electronic structure of the three-site SIAM employed in the present
work. The embedded/active regime including the impurity site and one bath site
is framed. Numbers $0-5$ label the six spin sites. The electron configuration
with one up-spin electron on site ‘0’ and two down-spin electrons on sites ‘3’
and ‘4’ is employed as the reference state.

The SIAM Hamiltonian comprises three parts

$\displaystyle H_{\rm SIAM}=H_{\rm imp.}+H_{\rm bath}+H_{hyb.},$ (57)

where

$\displaystyle H_{\rm
imp.}=\sum_{\sigma}(\epsilon_{c}-\mu)c^{\dagger}_{\sigma}c_{\sigma}+Uc^{\dagger}_{\uparrow}c_{\uparrow}c^{\dagger}_{\downarrow}c_{\downarrow}$
(58)

describes the impurity-site energy $\epsilon_{c}$ and the Coulomb interaction
$U$ between the electrons with opposite spins ($\sigma=\uparrow$ or
$\downarrow$) at the impurity site,

$\displaystyle H_{\rm bath}=\sum_{i=1,\sigma}^{N_{\rm
b}}\epsilon_{d,i}d^{\dagger}_{i,\sigma}d_{i,\sigma}$ (59)

describes the non-interacting bath site with $\epsilon_{d}$’s being the bath-
site energies, and

$\displaystyle H_{\rm hyb.}=\sum_{i=1,\sigma}^{N_{\rm
b}}V_{i}\big{(}c^{\dagger}_{\sigma}d_{i,\sigma}+d^{\dagger}_{i,\sigma}c_{\sigma}\big{)}$
(60)

describes the coupling between the impurity site and the bath levels due to
the hybridization. The electronic structure of the three-site SIAM is
exhibited in Fig. 2, where we select one impurity site and one bath site as
embedded region. The reference state is $|\Phi\rangle=|100110\rangle$ where we
include six spin sites with ‘1’ denoting occupied and ‘0’ denoting virtual
status of the corresponding spin site.

Determinants | $|\Phi_{0}^{1}\rangle$ | $|\Phi_{0}^{2}\rangle$ | $|\Phi_{3}^{5}\rangle$ | $|\Phi_{4}^{5}\rangle$ | $|\Phi_{03}^{15}\rangle$ | $|\Phi_{04}^{15}\rangle$ | $|\Phi_{03}^{25}\rangle$ | $|\Phi_{04}^{25}\rangle$
---|---|---|---|---|---|---|---|---
$t$-amplitudes | -0.628627 | 0.244428 | -0.244428 | -0.628627 | 0.013884 | 0.093685 | -0.003991 | -0.013884
$\lambda$-amplitudes | -0.442227 | 0.165380 | -0.165380 | -0.442227 | 0.075843 | 0.221301 | -0.028853 | -0.075843
$s$-amplitudes | -0.442227 | 0.165380 | -0.165380 | -0.442227 | 0.002707 | 0.025736 | -0.001502 | -0.002707
Table 2: CCSD amplitudes ($t$’s, $\lambda$’s, and $s$’s) of three-site SIAM
model Hamiltonian with $\epsilon_{c}=-0.5$, $\epsilon_{d}=\\{-1.0,1.0\\}$,
$U=1.0$, and $V_{i}=\\{1.0,1.0\\}$. Figure 3: SES-CC workflow to obtain the
ground state energy, effective hamiltonian $\overline{H}^{\rm eff}$ and
cluster amplitude corresponding to the subsystems. When converge, the energy
of each subsystem matches the energy of the entire system.

###  III.2 Ground state energy

Given a three-site SIAM Hamiltonian with $U=1.0$, $\epsilon_{c}=-0.5$,
$\epsilon_{d}=\\{-1.0,1.0\\}$, and $V_{i}=\\{1.0,1.0\\}$, $E_{\rm
0}=-3.7572543$ a.u. Note that since only single and double excitations could
occur in the model system, the CCSD calculation reproduces the true ground
state energy with the computed CC $T$ amplitudes given by

$\displaystyle T$ $\displaystyle=T_{\rm int}+T_{\rm ext}$ $\displaystyle
T_{\rm int}$ $\displaystyle=t_{0}^{1}a_{1}^{\dagger}a_{0}$ $\displaystyle
T_{\rm ext}$
$\displaystyle=t_{0}^{2}a_{2}^{\dagger}a_{0}+t_{3}^{5}a_{5}^{\dagger}a_{3}+t_{4}^{5}a_{5}^{\dagger}a_{4}+t_{03}^{15}a_{1}^{\dagger}a_{5}^{\dagger}a_{3}a_{0}$
$\displaystyle~{}~{}+t_{04}^{15}a_{1}^{\dagger}a_{5}^{\dagger}a_{4}a_{0}+t_{03}^{25}a_{2}^{\dagger}a_{5}^{\dagger}a_{3}a_{0}+t_{04}^{25}a_{2}^{\dagger}a_{5}^{\dagger}a_{4}a_{0}.$
(61)

All the CCSD amplitudes ($t$’s, $\lambda$’s, and $s$’s) are given in Tab. 2

Figure 4: (a) Energy deviation of SES-CCSD ground-state energy with respect to
exact ground-state energy SIAM and (b) double occupancy of the impurity
$\langle n_{c\uparrow}n_{c\downarrow}\rangle$ computed by SES-CCSD as
functions of $U$ for the symmetric SIAM with different numbers of external
subsystems. In the three-site SIAM Hamiltonian, we select $V=1$,
$\epsilon_{c}=-U/2$, $\epsilon_{d}=\\{-1.0,1.0\\}$.

###  III.3 SES-CC effective Hamiltonian

By employing $T_{\text{ext}}$, we can downfold the $2^{6}\times 2^{6}$ three-
site SIAM Hamiltonian to the $2\times 2$ effective Hamiltonian,
$\overline{H}^{\rm eff}_{\rm ext}$, using Eq. (12). In practical calculations,
we propose an iterative computational procedure for performing the SES-CC
embedding calculations, as demonstrated in Fig. 3. This approach is similar to
the recently proposed flow strategy Kowalski (2021); Kowalski and Bauman
(2023) but operates within the conventional nonunitary coupled cluster
framework. As depicted in Fig. 3, the proposed procedure comprises three basic
steps. Initially, the configuration space of the system of interest is
partitioned (potentially with overlaps) into subsystems based on various SESs.
Subsequently, the projected effective Hamiltonians associated with the SESs
are constructed and diagonalized. Finally, the eigenvector of each projected
effective Hamiltonian is mapped back to the corresponding global positions
associated with the coupled cluster amplitude tensor of the full system,
followed by an amplitude synchronization and update. It is important to note
that the second and third steps form a loop to iterate and update the coupled
cluster amplitude. Upon convergence, the eigenvalues of each projected
effective Hamiltonian become consistent and match the coupled cluster energy
of the full system. For current cast, upon converging, we obtained the
$2\times 2$ effective Hamiltonian as follows

$\displaystyle\overline{H}_{\rm ext}^{\rm eff}$ $\displaystyle=(P^{(N)}+Q_{\rm
int}^{(N)})\overline{H}_{\rm ext}(P^{(N)}+Q_{\rm int}^{(N)})$
$\displaystyle=\left(\begin{array}[]{cc}-3.12862715&1.0\\\
1.33811282&-1.62862715\end{array}\right)$ (64)
$\displaystyle=-2.37862715I+1.16905641\sigma_{x}$
$\displaystyle~{}~{}+0.16905641j\sigma_{y}-0.75\sigma_{z}$ (65)

where I is identity matrix, $\sigma_{x,y,z}$ are spin-1/2 Pauli matrices, and
we define

$\displaystyle Q_{\rm int}^{(N)}$
$\displaystyle=|\Phi_{0}^{1}\rangle\langle\Phi_{0}^{1}|,$ $\displaystyle
Q_{\rm ext}^{(N)}$
$\displaystyle=|\Phi_{0}^{2}\rangle\langle\Phi_{0}^{2}|+|\Phi_{3}^{5}\rangle\langle\Phi_{3}^{5}|+|\Phi_{4}^{5}\rangle\langle\Phi_{4}^{5}|+|\Phi_{03}^{15}\rangle\langle\Phi_{03}^{15}|$
$\displaystyle~{}~{}+|\Phi_{04}^{15}\rangle\langle\Phi_{04}^{15}|+|\Phi_{03}^{25}\rangle\langle\Phi_{03}^{25}|+|\Phi_{04}^{25}\rangle\langle\Phi_{04}^{25}|,$
(66)

with the superscripts (subscripts) in the bra/ket labelling virtual (occupied)
orbitals as showed in Fig. 2. Diagonalizing $\overline{H}_{\rm ext}^{\rm eff}$
gives rise to the lowest eigenvalue $E_{\rm CC}^{(N)}$ (i.e. $E_{0}$) and the
corresponding left and right eigenvectors, i.e $\langle\Psi_{\rm int}|$ and
$|\Psi_{\rm int}\rangle$

$\displaystyle|\Psi_{\rm int}\rangle$ $\displaystyle=e^{T_{\rm
int}}|\Phi\rangle=|\Phi\rangle+t_{0}^{1}|\Phi_{0}^{1}\rangle$ (67)

whose normalized form after projecting onto the active space gives

$\displaystyle\left(0.8497802,-0.53419491\right)^{\rm T}.$

Similarly, we have $\langle\Psi_{\rm int}|$ whose normalized form after
projecting onto the active space gives

$\displaystyle\left(0.90848166,-0.42679229\right).$

It’s worth mentioning that, based on the proposed workflow, different levels
of approximation can be introduced through the incrementation of the external
subsystem, which would be particularly useful when the external degrees of
freedom is too larger to deal with while the only limited number of them are
pivotal for the accurate interaction between the internal system of interest
and the external environment. For instance, simple demonstrations are shown in
Figs. 4 for the three-site SIAM, where four external subsystems are included
in an incremental manner to offer varying levels of approximation for energy
and double occupation.

###  III.4 Spectral function

Fig. 5 shows the computed the spectral function from exact diagonalization and
SES-CCSDGF approach with different $U/V$ ratios ranging from weakly to
strongly correlated scenarios. For the three-site model, since the highest
excitations are double, The SES-CCSDGF workflow need to include external
subsystems featuring double excitations from the refrence state in order to
reproduces the exact spectra.

Figure 5: $\omega$-dependent spin-polarized spectral function, $A(\omega)$, of
the three-site SIAM computed from exact diagonalization and CCSD approaches.
In the three-site SIAM Hamiltonian, we select $V=1$, $\epsilon_{c}=-U/2$,
$\epsilon_{d}=\\{-1.0,1.0\\}$.

##  IV Conclusion and outlook

In this paper, we have established a theoretical foundation for a novel and
rigorous sub-system embedding within coupled cluster Green’s function theory.
Distinct from previously reported Green’s function embedding methods Sun and
Chan (2016), our proposed theory relies on partitioning the excitation
manifold through sub-system embedding subalgebras (SES). It operates directly
on the $\omega$-dependent Green’s function formalism within the conventional
non-unitary coupled cluster theoretical framework. Our approach anticipates
offering a meticulous description of system-environment coupling by utilizing
existing coupled cluster machinery. Due to the non-unitary nature of
conventional coupled cluster theory, an essential prerequisite for the
proposed theory—which we have developed in this work—is the generalization of
the SES-CC formalism to address the left eigenvalue problem through a specific
form of Hamiltonian similarity transformation. This embedding approach is
inherently suitable for multiscale/multiphysics simulations, enabling the
iterative construction of multiple embedded sub-systems to balance accuracy
and efficiency optimally.

Future studies could see the proposed theory directly linked to advancements
in quantum computing. Specifically, an effective Hamiltonian—significantly
reduced in dimension compared to the original—could be employed and more
readily encoded within quantum computing frameworks for quantum dynamics
simulations, or for isolating segments of the exact spectrum, assuming there
are ample external subsystems to replicate the desired physics. Concerning the
Green’s function calculations within the proposed theoretical framework, the
necessity arises to evaluate the inverse of the shifted projected Hamiltonian
within the active space. To facilitate larger scale quantum simulations,
further refinements may be derived from the recently suggested connected
moment expansions (see Ref. 63 for a recent review and the references
therein), which have shown superior convergence performance compared to direct
Taylor or geometric expansions. The primary motivation for this simplification
is to employ straightforward and less complex circuits for expectation
measurement, optimizing the number of measurements required through advanced
techniques such as grouping or shadow techniques.

##  V Acknowledgements

B. P. and K. K. were supported by the “Embedding QC into Many-body Frameworks
for Strongly Correlated Molecular and Materials Systems” project, which is
funded by the U. S. Department of Energy, Office of Science, Office of Basic
Energy Sciences (BES), the Division of Chemical Sciences, Geosciences, and
Biosciences under FWP 72689.

##  VI Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  Appendix A Properties of $\Lambda_{\rm int}$ and $S_{\rm ext}$ operators
in the NSL

Let us assume that in the non-interacting sub-systems limit (NSL), the entire
composite system AB dissociates into sub-system A (defined by the active space
generated by SES $\mathfrak{h}$) and the remaining part referred to as the
sub-system B such that

$\displaystyle
H\xrightarrow{NSL}H_{A}+H_{B}~{},~{}|\Phi\rangle\xrightarrow{NSL}|\Phi_{A}\Phi_{B}\rangle~{},$
$\displaystyle E_{\rm CC}^{(N)}\xrightarrow{NSL}E_{A,\rm CC}+E_{B,\rm CC}~{},$
(68)

which leads to
$\overline{H}_{N}\xrightarrow{NSL}\overline{H}_{N,A}+\overline{H}_{N,B}$. In
the further analysis, due to of definition of the internal and external parts
of the operators, we will also assume that

$\displaystyle T\xrightarrow{NSL}T_{A}+T_{B}~{},T_{\rm
int}\xrightarrow{NSL}T_{A}~{},T_{\rm ext}\xrightarrow{NSL}T_{B}~{},$
$\displaystyle\Lambda_{\rm int}\xrightarrow{NSL}\Lambda_{A}~{}.$ (69)

As discussed in Section II.2, the equations for $\Lambda_{\rm int}$ and
$S_{\rm ext}$, in the exact wave function limit, can be cast in the form

$\displaystyle\langle\Phi|(1+\Lambda_{\rm int})[e^{S_{\rm
ext}}\overline{H}_{N}e^{-S_{\rm ext}}](P^{(N)}+Q_{\rm int}^{(N)}+Q_{\rm
ext}^{(N)})=0~{}.$ (70)

Given the fact that $S_{\rm ext}$ cluster amplitudes have to involve by
construction at least one spinorbital index on sub-system B, the general
decomposition of $S_{\rm ext}$ in the NSL takes the form

$S_{\rm ext}\xrightarrow{NSL}S_{B}+S_{AB}$ (71)

(where $S_{AB}$ is defined by the amplitudes indexed by both active and
inactive spin-orbitals), Eq. (70) can then be cast the form:

$\displaystyle\langle\Phi|(1+\Lambda_{A})\left[e^{S_{B}+S_{AB}}\left(\overline{H}_{N,A}+\overline{H}_{N,B}\right)e^{-S_{B}-S_{AB}}\right](P^{(N)}+Q_{A}^{(N)}+Q_{B}^{(N)}+Q_{AB}^{(N)})=0~{},$
(72)

which can be re-written in the form

$\displaystyle\langle\Phi|(1+\Lambda_{A})\left[\left(e^{S_{AB}}\overline{H}_{N,A}\right)_{c}+\left(e^{S_{AB}}\left(e^{S_{B}}\overline{H}_{N,B}\right)_{c}\right)_{c}\right]$
$\displaystyle(P^{(N)}+Q_{A}^{(N)}+Q_{B}^{(N)}+Q_{AB}^{(N)})=0~{},$ (73)

where the subscript ‘$c$’ denotes the connected parts of a given operator
expression. Since the CC equations are satisfied for the subsystem, i.e.,
$Q_{B}\overline{H}_{N,B}|\Phi\rangle$ and all uncontracted lines emanating
from $S_{AB}$ or $\overline{H}_{N,B}$ lead to vanishing expessions when acting
onto reference function $|\Phi\rangle=|\Phi_{A}\Phi_{B}\rangle$ and excited
configurations defining $Q_{A}$, the only non-zero contributions to the
$P+Q_{A}$ projections are given by

$\langle\Phi|(1+\Lambda_{A})\overline{H}_{N,A}(P^{(N)}+Q_{A}^{(N)})=0~{},$
(74)

which represents $\Lambda$-equation for isolated sub-system A. Similar
analysis leads to the form of (73) onto $Q_{B}$,

$\langle\Phi|\left(e^{S_{B}}\overline{H}_{N,B}\right)_{c}Q_{B}=0~{}.$ (75)

At this point, both (74) and (75), as well as their $\Lambda_{A}$ and $S_{B}$
solutions, are independent of the $S_{AB}$ part.

To determine $S_{AB}$ we re-cast Eq. (73) in projections on $Q_{AB}$ in the
form

$\langle\Phi|[\left(\left(e^{S_{AB}}-1\right)\overline{H}_{N,A}\right)_{c}+\left(\left(e^{S_{AB}}-1\right)\left(e^{S_{B}}\overline{H}_{N,B}\right)_{c}\right)_{c}]Q_{AB}=0~{},$
(76)

which represents equations homogeneous in $S_{AB}$ with $S_{AB}=0$ as a
solution. This proves that in the NSL

$S_{\rm ext}\xrightarrow{NSL}S_{B}~{}.$ (77)

##  Appendix B Properties of $X_{K}(\omega)$ and $Y_{L}(\omega)$ operators in
the NSL

In this appendix, we will discuss properties of $X_{K}(\omega)$ and
$Y_{L}(\omega)$ operators in the NLS. Although in previous papers Peng and
Kowalski (2016) the connected character of $X_{K}(\omega)$ and $Y_{L}(\omega)$
operators has been discussed, in this Appendix we will focus on explicit form
of the equations for these operators in the NSL. For the sake of brevity in
our analysis we focus on the $X_{K}(\omega)$. Results of this analysis can be
readily extended to the $Y_{L}(\omega)$ operator. Same as Appendix A, we will
use the asymptotic behavior of the following operators:

$\displaystyle T\xrightarrow{NSL}T_{A}+T_{B}~{},~{}T_{\rm
int}\xrightarrow{NSL}T_{A}~{},~{}T_{\rm ext}\xrightarrow{NSL}T_{B}~{}.$ (78)

Employ $K_{A}$ to denote the active spin-orbital $K$ in the NSL, the equations
for $X_{K_{A}}(\omega)$ take the form

$\displaystyle
Q^{(N-1)}_{A}[(\omega_{-i0}+\overline{H}_{N,A}+\overline{H}_{N,B})(X_{K_{A},A}(\omega)+X_{K_{A},B}(\omega)+X_{K_{A},AB}(\omega))]_{c}|\Phi\rangle$
$\displaystyle=Q^{(N-1)}_{A}\overline{a_{K_{A}}}|\Phi\rangle~{},$ (79)
$\displaystyle
Q^{(N-1)}_{B}[(\omega_{-i0}+\overline{H}_{N,A}+\overline{H}_{N,B})(X_{K_{A},A}(\omega)+X_{K_{A},B}(\omega)+X_{K_{A},AB}(\omega))]_{c}|\Phi\rangle$
$\displaystyle=Q^{(N-1)}_{B}\overline{a_{K_{A}}}|\Phi\rangle~{},$ (80)
$\displaystyle
Q^{(N-1)}_{AB}[(\omega_{-i0}+\overline{H}_{N,A}+\overline{H}_{N,B})(X_{K_{A},A}(\omega)+X_{K_{A},B}(\omega)+X_{K_{A},AB}(\omega))]_{c}|\Phi\rangle$
$\displaystyle=Q^{(N-1)}_{AB}\overline{a_{K_{A}}}|\Phi\rangle~{},$ (81)

where in analogy to analysis in Appendix A

$\displaystyle X_{K,{\rm int}}(\omega)$
$\displaystyle\xrightarrow{NSL}X_{K_{A},A}(\omega)~{},$ (82) $\displaystyle
X_{K,{\rm ext}}(\omega)$
$\displaystyle\xrightarrow{NSL}X_{K_{A},B}(\omega)+X_{K_{A},AB}(\omega)~{}.$
(83)

Here $X_{K_{A},A}(\omega)$, $X_{K_{A},B}(\omega)$, and $X_{K_{A},AB}(\omega)$
are the operators defined by the excitations localized on sub-system A, B, and
the excitations that involve both sub-systems simultaneously, respectively.
The operators $Q_{A}^{(N-1)}$, $Q_{B}^{(N-1)}$, and $Q_{AB}^{(N-1)}$ are
defined in an analogous way in the $(N-1)$-electron Hilbert space. Given the
fact that creation/annihilation operators corresponding to electronic states
on different sub-systems anti-commute, one can write that

$\overline{a_{K_{A}}}=e^{-T_{A}}a_{K_{A}}e^{T_{A}}~{},$ (84)

and the fact that only connected expressions enter $X_{K_{A}}(\omega)$
equations, the Eqs. (79-81) can be simplified to

$\displaystyle
Q^{(N-1)}_{A}[(\omega_{-i0}+\overline{H}_{N,A})X_{K_{A},A}(\omega)+\overline{H}_{N,B}X_{K_{A},AB}(\omega))]_{c}|\Phi\rangle$
$\displaystyle=Q^{(N-1)}_{A}\overline{a_{K_{A}}}|\Phi\rangle~{},$ (85)
$\displaystyle
Q^{(N-1)}_{B}[(\omega_{-i0}+\overline{H}_{N,B})X_{K_{A},B}(\omega)+(\overline{H}_{N,A}+\overline{H}_{N,B})X_{K_{A},AB}(\omega)]_{c}|\Phi\rangle$
$\displaystyle=0~{},$ (86) $\displaystyle
Q^{(N-1)}_{AB}[(\omega_{-i0}+\overline{H}_{N,A}+\overline{H}_{N,B})X_{K_{A},AB}(\omega)]_{c}|\Phi\rangle$
$\displaystyle=0~{}.$ (87)

A close inspection of Eqs. (86-87) shows that both $X_{K_{A},AB}(\omega)$ and
$X_{K_{A},B}(\omega)$ vanish, i.e.,

$\displaystyle X_{K_{A},AB}(\omega)$ $\displaystyle=0~{},$ (88) $\displaystyle
X_{K_{A},B}(\omega)$ $\displaystyle=0~{},$ (89)

and thus

$\displaystyle X_{K,\rm
ext}\xrightarrow[]{NSL}0~{},~{}X_{K}\xrightarrow[]{NSL}X_{K_{A},A}$ (90)

for non-resonant frequencies $\omega$. This implies that the equation for
$X_{K_{A},A}(\omega)$ takes the form

$Q^{(N-1)}_{A}[(\omega_{-i0}+\overline{H}_{N,A})X_{K_{A},A}(\omega)]_{c}|\Phi\rangle=Q^{(N-1)}_{A}\overline{a_{K_{A}}}|\Phi\rangle~{},$
(91)

and is identical to the equation for $X_{K_{A}}(\omega)$’s of the isolated
sub-system A. Analogous results can be derived for the $Y_{L}(\omega)$
operator, which in the NSL is localized only on sub-system A, i.e.,

$\displaystyle Y_{L}(\omega)\xrightarrow{NSL}Y_{L_{A},A}(\omega),$ (92)

where $Y_{L_{A},A}(\omega)$ operator satisfies the equations:

$Q^{(N+1)}_{A}[(\omega_{+i0}-\overline{H}_{N,A})Y_{L_{A},A}(\omega)]_{c}|\Phi\rangle=Q^{(N+1)}_{A}\overline{a^{\dagger}_{L_{A}}}|\Phi\rangle~{},$
(93)

with $\overline{a^{\dagger}_{L_{A}}}=e^{-T_{A}}a^{\dagger}_{L_{A}}e^{T_{A}}$.
This analysis also illustrates the size-intensivity of the $X_{K}(\omega)$ and
$Y_{L}(\omega)$ operators.

##  Appendix C Asymptotic properties of sub-systems CC Green’s function

Using properties of the $T_{\rm int}$, $T_{\rm ext}$, $\Lambda_{\rm int}$,
$S_{\rm ext}$, $X_{K}(\omega)$, and $Y_{L}(\omega)$ operators derived in the
Appendices A and B, we can derive properties of the sub-block of the CC
Green’s function matrix. Following Eqs. (69,77), (90,92), one can derive the
form of $G_{KL}(\omega)$ block in the NSL as:

$\displaystyle G_{KL}(\omega)\xrightarrow{NSL}$
$\displaystyle~{}~{}~{}\langle\Phi|(1+\Lambda_{A})\overline{a^{\dagger}_{L_{A},A}}X_{K_{A},A}(\omega)|\Phi\rangle$
$\displaystyle+\langle\Phi|(1+\Lambda_{A})\overline{a_{K_{A},A}}Y_{L_{A},A}(\omega)|\Phi\rangle~{},$
(94)

which is identical with the algebraic form of the CC Green’s function for
isolated sub-system A.

Similarly, for other blocks in the CC Green’s function matrix, we have

$\displaystyle G_{\overline{p}\overline{q}}(\omega)\xrightarrow{NSL}$
$\displaystyle~{}~{}~{}\langle\Phi|(1+\Lambda_{B})\overline{a^{\dagger}_{\overline{q}_{B}}}_{,B}X_{\overline{p},B}(\omega)|\Phi\rangle$
$\displaystyle+\langle\Phi|(1+\Lambda_{B})\overline{a_{\overline{p}_{B}}}_{,B}Y_{\overline{q},B}(\omega)|\Phi\rangle~{},$
(95) $\displaystyle G_{K\overline{p}}(\omega)\xrightarrow{NSL}$
$\displaystyle~{}~{}0~{},~{}G_{\overline{p}K}(\omega)\xrightarrow{NSL}~{}0.$
(96)

Therefore, in the NSL the $\mathbf{G}$ of the composite system AB simply
becomes a diagonal block matrix

$\displaystyle\mathbf{G}^{\rm\textbf{AB}}(\omega)\xrightarrow{NSL}\left[\begin{array}[]{c|c}\mathbf{G}^{\rm\textbf{A}}(\omega)&\textbf{0}\\\
\hline\cr\textbf{0}&\mathbf{G}^{\rm\textbf{B}}(\omega)\end{array}\right].$
(99)

## References

  * Breuer  _et al._ (2002) H. Breuer, F. Petruccione, and S. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2002).
  * Weiss (2012) U. Weiss,  _Quantum Dissipative Systems_, 4th ed. (WORLD SCIENTIFIC, 2012).
  * May and KüHn (2011) V. May and O. KüHn, eds., in  _Charge and Energy Transfer Dynamics in Molecular Systems_ (John Wiley & Sons, Ltd, 2011).
  * Wasielewski  _et al._ (2020) M. R. Wasielewski, M. D. E. Forbes, N. L. Frank, K. Kowalski, G. D. Scholes, J. Yuen-Zhou, M. A. Baldo, D. E. Freedman, R. H. Goldsmith, T. G. III, M. L. Kirk, J. K. McCusker, J. P. Ogilvie, D. A. Shultz, S. Stoll, and K. B. Whaley, Nat. Rev. Chem. 4, 490 (2020).
  * Wesolowski  _et al._ (2015) T. A. Wesolowski, S. Shedge, and X. Zhou, Chem. Rev. 115, 5891 (2015).
  * Sun and Chan (2016) Q. Sun and G. K.-L. Chan, Acc. Chem. Res. 49, 2705 (2016).
  * Jones  _et al._ (2020) L. O. Jones, M. A. Mosquera, G. C. Schatz, and M. A. Ratner, J. Am. Chem. Soc. 142, 3281 (2020).
  * Sushko  _et al._ (2018) P. V. Sushko, C. Huang, N. Govind, and K. Kowalski, “Embedding methods in materials discovery,” in  _Computational Materials Discovery_, edited by A. R. Oganov, G. Saleh, and A. G. Kvashnin (Royal Society of Chemistry, 2018) Chap. 4, pp. 87–116.
  * Libisch  _et al._ (2014) F. Libisch, C. Huang, and E. A. Carter, Acc. Chem. Res. 47, 2768 (2014).
  * Seijo and Barandiarán (1999) L. Seijo and Z. Barandiarán, “The <italic>ab initio</italic> model potential method: A common strategy for effective core potential and embedded cluster calculations,” in  _Computational Chemistry: Reviews of Current Trends_ (World Scientific, 1999) pp. 55–152.
  * Jacob and Neugebauer (2014) C. R. Jacob and J. Neugebauer, WIREs Comput. Mol. Sci. 4, 325 (2014).
  * Warshel and Levitt (1976) A. Warshel and M. Levitt, J. Mol. Bio. 103, 227 (1976).
  * Brunk and Rothlisberger (2015) E. Brunk and U. Rothlisberger, Chem. Rev. 115, 6217 (2015).
  * Svensson  _et al._ (1996) M. Svensson, S. Humbel, R. D. J. Froese, T. Matsubara, S. Sieber, and K. Morokuma, J. Phys. Chem. 100, 19357 (1996).
  * Chung  _et al._ (2015) L. W. Chung, W. M. C. Sameera, R. Ramozzi, A. J. Page, M. Hatanaka, G. P. Petrova, T. V. Harris, X. Li, Z. Ke, F. Liu, H.-B. Li, L. Ding, and K. Morokuma, Chem. Rev. 115, 5678 (2015).
  * Govind  _et al._ (1998) N. Govind, Y. Wang, A. da Silva, and E. Carter, Chem. Phys. Lett. 295, 129 (1998).
  * Huang  _et al._ (2011) C. Huang, M. Pavone, and E. A. Carter, J. Chem. Phys. 134, 154110 (2011).
  * Manby  _et al._ (2012) F. R. Manby, M. Stella, J. D. Goodpaster, and T. F. Miller, J. Chem. Theory Comput. 8, 2564 (2012).
  * Elliott  _et al._ (2010) P. Elliott, K. Burke, M. H. Cohen, and A. Wasserman, Phys. Rev. A 82, 024501 (2010).
  * Nafziger and Wasserman (2014) J. Nafziger and A. Wasserman, J. Phys. Chem. A 118, 7623 (2014).
  * Huang and Carter (2011) C. Huang and E. A. Carter, J. Chem. Phys. 135, 194104 (2011).
  * Fornace  _et al._ (2015) M. E. Fornace, J. Lee, K. Miyamoto, F. R. Manby, and T. F. Miller, J. Chem. Theory Comput. 11, 568 (2015).
  * Inglesfield (1981) J. E. Inglesfield, J. Phys. C: Solid State Phys. 14, 3795 (1981).
  * Chibani  _et al._ (2016) W. Chibani, X. Ren, M. Scheffler, and P. Rinke, Phys. Rev. B 93, 165106 (2016).
  * Lan  _et al._ (2015) T. N. Lan, A. A. Kananenka, and D. Zgid, J. Chem. Phys. 143, 241102 (2015).
  * Rusakov  _et al._ (2019) A. A. Rusakov, S. Iskakov, L. N. Tran, and D. Zgid, J. Chem. Theory Comput. 15, 229 (2019).
  * Knizia and Chan (2013) G. Knizia and G. K.-L. Chan, J. Chem. Theory Comput. 9, 1428 (2013).
  * Li  _et al._ (2019) W. Li, M. Chen, E. Rabani, R. Baer, and D. Neuhauser, J. Chem. Phys. 151, 174115 (2019).
  * Kowalski (2018) K. Kowalski, J. Chem. Phys. 148, 094104 (2018), https://doi.org/10.1063/1.5010693 .
  * Arponen (1983) J. Arponen, Ann. Phys. 151, 311 (1983).
  * Nooijen and Snijders (1992) M. Nooijen and J. G. Snijders, Int. J. Quantum Chem. 44, 55 (1992).
  * Nooijen and Snijders (1993) M. Nooijen and J. G. Snijders, Int. J. Quantum Chem. 48, 15 (1993).
  * Nooijen and Snijders (1995) M. Nooijen and J. G. Snijders, J. Chem. Phys. 102, 1681 (1995).
  * Peng  _et al._ (2021) B. Peng, N. P. Bauman, S. Gulania, and K. Kowalski (Elsevier, 2021) pp. 23–53.
  * Purvis and Bartlett (1982) G. Purvis and R. Bartlett, J. Chem. Phys. 76, 1910 (1982).
  * Noga and Bartlett (1987) J. Noga and R. J. Bartlett, J. Chem. Phys. 86, 7041 (1987), https://doi.org/10.1063/1.452353 .
  * Noga and Bartlett (1988) J. Noga and R. J. Bartlett, J. Chem. Phys. 89, 3401 (1988), https://doi.org/10.1063/1.455742 .
  * Scuseria and Schaefer (1988) G. E. Scuseria and H. F. Schaefer, Chem. Phys. Lett. 152, 382 (1988).
  * Kucharski and Bartlett (1991) S. Kucharski and R. Bartlett, Theoret. Chim. Acta 80, 387 (1991), https://doi.org/10.1007/BF01117419 .
  * Oliphant and Adamowicz (1991) N. Oliphant and L. Adamowicz, J. Chem. Phys. 95, 6645 (1991), https://doi.org/10.1063/1.461534 .
  * Bauman  _et al._ (2019) N. P. Bauman, E. J. Bylaska, S. Krishnamoorthy, G. H. Low, N. Wiebe, C. E. Granade, M. Roetteler, M. Troyer, and K. Kowalski, J. Chem. Phys. 151, 014107 (2019), https://doi.org/10.1063/1.5094643 .
  * Metcalf  _et al._ (2020) M. Metcalf, N. P. Bauman, K. Kowalski, and W. A. de Jong, J. Chem. Theory Comput. 16, 6165 (2020), pMID: 32915568, https://doi.org/10.1021/acs.jctc.0c00421 .
  * Bauman  _et al._ (2020) N. P. Bauman, J. Chládek, L. Veis, J. Pittner, and K. Kowalski, arXiv (2020), 2011.01985 [quant-ph] .
  * Szalay  _et al._ (1995) P. G. Szalay, M. Nooijen, and R. J. Bartlett, J. Chem. Phys. 103, 281 (1995), https://doi.org/10.1063/1.469641 .
  * Piecuch and Bartlett (1999) P. Piecuch and R. J. Bartlett, Adv. Quantum Chem. 34, 295 (1999).
  * Salter  _et al._ (1987) E. A. Salter, H. Sekino, and R. J. Bartlett, J. Chem. Phys. 87, 502 (1987).
  * Stanton and Bartlett (1993) J. F. Stanton and R. J. Bartlett, J. Chem. Phys. 98, 7029 (1993).
  * Monkhorst (1977) H. J. Monkhorst, Int. J. Quantum Chem. 12, 421 (1977).
  * Koch and Jørgensen (1990) H. Koch and P. Jørgensen, J. Chem. Phys. 93, 3333 (1990).
  * Kowalski  _et al._ (2014) K. Kowalski, K. Bhaskaran-Nair, and W. A. Shelton, J. Chem. Phys. 141, 094102 (2014).
  * Bhaskaran-Nair  _et al._ (2016) K. Bhaskaran-Nair, K. Kowalski, and W. A. Shelton, J. Chem. Phys. 144, 144101 (2016).
  * Peng and Kowalski (2016) B. Peng and K. Kowalski, Phys. Rev. A 94, 062512 (2016).
  * Peng and Kowalski (2018) B. Peng and K. Kowalski, J. Chem. Theory Comput. 14, 4335 (2018).
  * Helgaker  _et al._ (2014) T. Helgaker, P. Jorgensen, and J. Olsen,  _Molecular Electronic-Structure Theory_ (Wiley, 2014).
  * Schirmer and Mertins (2010) J. Schirmer and F. Mertins, Theor. Chem. Acc 125, 145 (2010).
  * Kowalski (2021) K. Kowalski, Phys. Rev. A 104, 032804 (2021).
  * Kowalski and Bauman (2020) K. Kowalski and N. P. Bauman, J. Chem. Phys. 152, 244127 (2020).
  * Bauman and Kowalski (2021) N. P. Bauman and K. Kowalski, “Coupled cluster downfolding theory: towards efficient many-body algorithms for dimensionality reduction of composite quantum systems,” (2021), arXiv:2111.03215 [quant-ph] .
  * Anderson (1961) P. W. Anderson, Phys. Rev. 124, 41 (1961).
  * Newns (1969) D. M. Newns, Phys. Rev. 178, 1123 (1969).
  * Wilson (1975) K. G. Wilson, Rev. Mod. Phys. 47, 773 (1975).
  * Kowalski and Bauman (2023) K. Kowalski and N. P. Bauman, arXiv (2023).
  * Aulicino  _et al._ (2022) J. C. Aulicino, T. Keen, and B. Peng, Int. J. Quantum Chem. 122, e26853 (2022).