# Exploring Ground and Excited States via Single Reference Coupled-Cluster
Theory and Algebraic Geometry

Svala Sverrisdóttir  Department of Mathematics, The University of California,
Berkeley, CA, USA  Fabian M. Faulstich  Department of Mathematics, Rensselaer
Polytechnic Institute, Troy, NY, USA<EMAIL_ADDRESS>
###### Abstract

The exploration of the root structure of coupled cluster equations holds both
foundational and practical significance for computational quantum chemistry.
This study provides insight into the intricate root structures of these non-
linear equations at both the CCD and CCSD level of theory. We utilize
computational techniques from algebraic geometry, specifically the monodromy
and parametric homotopy continuation methods, to calculate the full solution
set. We compare the computed CC roots against various established theoretical
upper bounds, shedding light on the accuracy and efficiency of these bounds.
We hereby focus on the dissociation processes of four-electron systems such as
(H2)2 in both D2h and D∞h configurations, H4 symmetrically distorted on a
circle, and lithium hydride. We moreover investigate the ability of single-
reference coupled cluster solutions to approximate excited state energies. We
find that multiple CC roots describe energies of excited states with high
accuracy. Remarkably, our investigations reveal that for systems like lithium
hydride, CC not only provides high-accuracy approximations to several excited
state energies but also to the states themselves.

##  1 Introduction

Coupled cluster (CC) theory is a high-precision wave function method in
computational quantum chemistry 1. Its origins date back to 1958 when Coester
proposed to use an exponential parametrization of the wave function 2, i.e.

$|{\Psi}\rangle=e^{\hat{T}}|{\Phi_{0}}\rangle,$ (1)

where $\hat{T}$ is the new unknown, and $|{\Phi_{0}}\rangle$ is a chosen
reference Slater determinant (commonly the Hartree-Fock state). This
parametrization was independently derived by Hubbard 3 and Hugenholtz 4 in
1957 as an alternative to summing many-body perturbation theory contributions
order by order. A pivotal moment in the development of CC theory was Čížek’s
seminal work in 1966 5. Čížek presented the first derivation of the CC working
equations together with corresponding simulations. This work established CC
theory within the theoretical framework that we recognize today, including
second quantization applied to many-fermion systems, normal ordering,
contractions, Wick’s theorem, normal-ordered Hamiltonians and Goldstone-style
diagrammatic techniques. There exist several reviews of CC theory, summarizing
different aspects of the method’s success over the past decades, see e.g.
Kümmel 6, Čížek 7, Bartlett 8, Paldus 9, Arponen 10 and Bishop 11.

Employing the exponential ansatz in Eq. (1) yields the following reformulation
of the stationary Schrödinger equation:

$\left\\{\begin{aligned}
\langle{\Phi_{0}}|e^{-\hat{T}}\hat{H}e^{\hat{T}}|{\Phi_{0}}\rangle&=E,\\\
\langle{\Phi}|e^{-\hat{T}}\hat{H}e^{\hat{T}}|{\Phi_{0}}\rangle&=0&\forall|{\Phi}\rangle\perp|{\Phi_{0}}\rangle.\\\
\end{aligned}\right.$ (2)

Employing the standard Galerkin projection, the latter condition yields the
projected CC equations (vide infra), i.e. a set of nonlinear algebraic
equations. Typically, a single solution to this set of equations is computed
using (quasi) Newton-type methods 12. However, as a set of nonlinear algebraic
equations, the CC equations possess multiple roots, raising the following
natural questions:

  * i)

How many solutions to the CC equations exist?

  * ii)

Which solution describes the ground state best?

  * iii)

Do some of the solutions describe excited states?

Especially the last question is controversial, to say the least. On one side,
it can be proven that as long a state is not strictly orthogonal to the chosen
reference state, CC theory can – in its untruncated form – describe that state
exactly. On the other side, excited states commonly have poor overlap with the
Hartree-Fock ground state reference – hindering the convergence of (quasi)
Newton-type methods. Besides such theoretical arguments, it has been
numerically shown that higher-order CC solutions can approximate excited state
energies well while the CC solutions provide poor approximations to the actual
states 13. Conventionally, excited states are computed in a post-single-
reference CC manner. This is done by either constructing and diagonalizing the
similarity-transformed Hamiltonian in equation-of-motion CC 14, 15, 16, 17,
18, 19, 20 or by examining the poles of the linear-response function in
linear-response CC 16, 21, 22, 23. In either case, the excited states provided
by these methods are inherently biased toward the lowest eigenstate upon which
they are built. As a direct consequence, the accuracy of such approaches is
known to be quite poor for certain classes of eigenstates that are
energetically higher than the ground state 24, 25, 26, 27, 28, 29.

Besides the fundamental nature of questions i)–iii), the existence of multiple
roots to the CC equations can present challenges for the practical
implementations of the method. For example, the roots may be difficult to
adequately converge to 5, 30, 31, 32, and the convergence properties of the
iterative solution methods can strongly depend on the employed initial guess
33, 34. See also Ref. 35 for a mathematical exposition of this problem and
Ref. 36 for a potential remedy.

In this work, we elaborate on recent advances coming from algebraic geometry
that provide new ideas on how to bound the number of roots to the CC
equations. These new and significantly improved bounds can then be used to
compute all roots to the CC equations, answering the above-posed questions i)
– iii). In particular, we demonstrate that higher-order solutions of the CC
equations can effectively approximate not only the energies of excited states
but also the eigenstates themselves.

###  1.1 Previous works and perspective

We highlight that the inaugural investigation of the root structure of CC
methods dates back to 1978 when Zivkovic and Monkhorst scrutinized the
singularities and multifarious solutions in single-reference CC equations 37.
Building upon this work, Paldus and colleagues subsequently conducted
mathematical and numerical analyses during the early 1990s, focusing on the
manifold of solutions of both single-reference and state-universal multi-
reference CC equations, while also elucidating their singularities and
analytical characteristics 38, 39. In 1998, Kowalski and Jankowski revived
homotopy methods in connection with the single-reference CC theory and used
them to solve a CCD system 40. Kowalski with Piecuch then extended the
application of homotopy methods to the CCSD, CCSDT, and CCSDTQ 41 for a
4-electron system described in a minimal basis set. From these investigations,
they derived the $\beta$-nested equations and proved the Fundamental Theorem
of the $\beta$-NE Formalism 42, which enabled an explanation of curves
connecting multiple solutions of the various CC polynomial systems, i.e. from
CCSD to CCSDT, CCSDT to CCSDTQ, etc. This Kowalski-Piecuch homotopy was
recently analyzed in-depth using topological degree theory 43, 44. Moreover,
Kowalski with Piecuch used homotopy methods to obtain all solutions of the
generalized Bloch equation 45 and explain symmetry breaking processes in the
CC formulation 46, 47, 48, 49, 50. It is noteworthy that these fundamental
results of CC theory have spurred the development of the method of moment CC
approach 51, 52, 53, 54, which in turn has led to the highly efficient
completely renormalized CC approach 55, 56, 57. This emphasizes the pivotal
role of this research trajectory in driving methodological innovations within
CC theory. A perspective on the recent use of homotopy methods in the context
of quantum chemistry can be found in Ref. 58.

A more algebraic and mathematically driven approach to this problem has
recently gained momentum. Initialized by some of the authors, a first
algebraic perspective was introduced in Ref. 59. This was followed by an in-
depth algebraic investigation leading to the introduction of the CC truncation
varieties 60 – which are the subject of this work. This work moreover
contained open mathematical questions and conjectures, one of which was proven
shortly after 61. More generally, these works are part of a presently
increasing interest in the mathematical understanding of CC theory, driven by
the applied mathematics community. In this context, Schneider performed the
first local analysis of CC theory in 2009 employing a functional analytic
framework 62. Subsequently, this approach was refined 63, 64 and applied to
different CC methods 65, 66, 67. A more versatile framework for analyzing
general CC variants using topological degree theory was later proposed by
Csirik and Laestadius 43, 44. The most recent numerical analysis results
regarding single reference CC were established by Hassan, Maday, and Wang 68,
69.

##  2 Coupled-cluster theory

The starting point of CC theory is the exponential ansatz in Eq. (1), leading
to the CC energy expression

$E=\langle{\Phi_{0}}|e^{-\hat{T}}\hat{H}e^{\hat{T}}|{\Phi_{0}}\rangle,$ (3)

where the cluster operator $\hat{T}$ is determined through the CC equations

$0=\langle{\Phi}|e^{-\hat{T}}\hat{H}e^{\hat{T}}|{\Phi_{0}}\rangle,\qquad\forall|{\Phi}\rangle\perp|{\Phi_{0}}\rangle.$
(4)

The cluster operator is a linear combination of elementary particle-hole
excitation operators, i.e.

$\hat{T}({\bf t})=\sum_{\mu>0}t_{\mu}\hat{X}_{\mu},$ (5)

where we highlight the dependence of $\hat{T}$ on its expansion coefficients
by writing $\hat{T}(\bf t)$. Recall that elementary particle-hole excitation
operators $\hat{X}_{\mu}$ are merely compositions of projections onto a subset
of single-particle basis elements (the occupied orbitals) and wedge products
with another subset of single-particle basis elements (the virtual orbital).
It is therefore convenient to label the elementary particle-hole excitation
operators using multi-indices $\mu>0$ that clarify the projections and wedge
products involved 12, 62. Therefore, acting on the reference state
$|{\Phi_{0}}\rangle$, the elementary particle-hole excitation operators define
the Hilbert space of functions that are $L^{2}$-orthogonal to
$|{\Phi_{0}}\rangle$. We henceforth set
$\mathcal{V}:=\\{|{\Phi_{0}}\rangle\\}^{\perp}$ and note that $\hat{T}({\bf
t})|{\Phi_{0}}\rangle\in\mathcal{V}$ for all CC amplitudes $\mathbf{t}$.
Hence, we can express Eq. (4) in the more common form as a set of projective
equations, i.e.

$0=\langle\Phi_{\mu}|e^{-\hat{T}({\bf t})}\hat{H}e^{\hat{T}({\bf
t})}|\Phi_{0}\rangle,\quad\forall\mu>0,$ (6)

where $|{\Phi_{\mu}}\rangle=\hat{X}_{\mu}|{\Phi_{0}}\rangle$ are the excited
Slater determinants spanning $\mathcal{V}$. In practice, the expansion in Eq.
(5) is commonly limited to e.g. only containing one and two particle-hole
excitations leading to CCSD. In this case, the projective equations (6) are
also restricted to the corresponding set of excited Slater determinants,
yielding a square system of polynomial equations that we seek to solve. Hence,
a key object in CC theory is the set

$\mathcal{S}=\\{{\bf
t}\in\mathbb{F}^{K}~{}|~{}\langle\Phi_{\mu}|e^{-\hat{T}({\bf
t})}\hat{H}e^{\hat{T}({\bf t})}|\Phi_{0}\rangle=0,\quad\forall\mu>0\\},$ (7)

where $\mathbb{F}$ is the considered number field (either $\mathbb{R}$ or
$\mathbb{C}$) and $K$ is the “system size” (determined by the number of
correlated electrons, the size of one-particle basis functions as well as
further selection rules). Note that $\mathbb{F}$ can be $\mathbb{R}$ or
$\mathbb{C}$ depending on whether we are seeking real- or complex-valued CC
amplitudes. We emphasize that although the CC polynomial coefficients are
real, the roots of the polynomial system may not be. Mathematically, the
theory describing complex-valued solutions is simpler and more complete than
the theory describing real-valued solutions. Therefore, we shall consider
$\mathbb{F}=\mathbb{C}$ for mathematical considerations, however, solving the
truncated CC equations over $\mathbb{F}=\mathbb{C}$ may yield unwanted complex
valued energies.

Mathematically, the set $\mathcal{S}$ is an (algebraic) variety. Since the
concept of varieties may be unconventional within the computational chemistry
community, we provide an illustrative example. A variety is a mathematical
object defined by the zeros of a set of polynomial equations. Consider the
multivariate polynomial system

$\left\\{\begin{aligned} p_{1}(x,y,z)&=x^{2}+y^{2}-z^{2},\\\
p_{2}(x,y,z)&=\frac{x+y}{\sqrt{2}}+z-1.\end{aligned}\right.$ (8)

The zeros of $p_{1}$ form a cone, whereas the zeros of $p_{2}$ form a plane;
the set of simultaneous roots defines a variety that corresponds to a
parabolic curve, see Fig. 1. As a curve, this variety is one-dimensional,
embedded in the ambient space which is three-dimensional in this case. A
central quality for this work is the degree of a variety which is defined to
be the maximal number of intersection points of the variety with a general
linear space of complementary dimension. By requiring the linear space to be
general we ensure the intersection is always finite. For the variety in Fig. 1
we see that intersecting the yellow curve with a plane (a linear space of
complimentary dimension) yields at most two intersection points, hence, its
degree is two. We add that over the complex field, the condition of “maximal
number of intersection points” is redundant for the definition of degree since
the number of intersection points is the same for generic complex spaces of
complementary dimension.

Figure 1: Root structure corresponding to Eq. (8). The cone describes the
roots of $p_{1}$, and the plane describes the roots of $p_{2}$. The
simultaneous roots are described by the yellow curve that defines the solution
variety of the multivariate polynomial system.

We conclude this section by noting that in the case of untruncated CC,
amplitudes that solve the full CC equations describe eigenstates of the
Hamiltonian. Therefore, the cardinality of $\mathcal{S}$ is exactly the number
of eigenstates of $\hat{H}$ that are intermediately normalized and the number
of roots to the CC equations is bounded from above by the number of Slater
determinants. However, when truncations are imposed, the cardinality of
$\mathcal{S}$ increases, i.e. truncated CC theory yields (some) unphysical
solutions. Therefore, it becomes less clear what the different elements in
$\mathcal{S}$ describe. Understanding the variety $\mathcal{S}$ and
characterizing (some of) its elements are the subject of this manuscript.

##  3 Truncation varieties

In CC theory one encounters various algebraic equations that define different
varieties. The truncation varieties 60 arise from the exponential
parametrization, i.e.

${\rm exp}:\mathbb{V}\to\mathcal{V},~{}\mathbf{t}\mapsto
e^{\hat{T}(\mathbf{t})}|\Phi_{0}\rangle,$ (9)

where $\mathbb{V}$ is the CC amplitude space. In the untruncated CC method,
the exponential map is bijective. However, when truncations are introduced,
this property no longer holds. Truncations correspond to restrictions of the
exponential map to a subspace, e.g. for CCSD the exponential map is restricted
to $\mathbb{V}_{\\{1,2\\}}\subset\mathbb{V}$ consisting of cluster amplitudes
with only one and two particle-hole excitations. The corresponding truncation
variety, $V_{\\{1,2\\}}$ is then defined as the smallest variety containing
${\rm exp}(\mathbb{V}_{\\{1,2\\}})$. This construction of the truncation
variety can be generalized to any CC truncation: Let $\sigma$ denote the CC
truncation rank, i.e. $\sigma\subset[\\![N]\\!]=\\{1,\dots,N\\}$, where $N$
denotes the number of electrons in the system. Then, the smallest variety
containing ${\rm exp}(\mathbb{V}_{\sigma})$ is the $\sigma$-truncation variety
denoted $V_{\sigma}$.

Remarkably, the truncation varieties can be used to express the projective CC
equations (6) as an eigenvalue problem, i.e.

$[\hat{H}|\Phi\rangle]_{\sigma}=E[|\Phi\rangle]_{\sigma},\quad|\Phi\rangle\in
V_{\sigma},$ (10)

where $[\,\cdot\,]_{\sigma}$ is the projection onto coordinates with
excitation level in $\sigma\cup\\{0\\}$, i.e.

$[|\Psi\rangle]_{\sigma}=\langle\Phi_{\mu}|\Psi\rangle\quad{\rm
for}~{}|\mu|\in\sigma\cup\\{0\\}.$ (11)

This formulation is indeed equivalent to the projective CC equations in Eq.
(6) for CCD, as well as for rank complete CC truncations, i.e. CCSD, CCSDT,
CCSDTQ etc., see Theorem 5.11 in Ref. 60. The number of roots for these CC
variants is bounded by the number of solutions to Eq. (10) for a generic
Hamiltonian matrix $\hat{H}$. We denote this number the CC degree of the
corresponding truncation variety. The CC degree serves as a complexity measure
for solving the CC equations and Bézout’s theorem provides the upper bound

$\mathrm{CCdeg(\sigma)}\leq(\dim(V_{\sigma})+1)\deg(V_{\sigma}),$ (12)

where $\dim(V_{\sigma})$ is the dimension of the truncation variety and
$\deg(V_{\sigma})$ is its degree. The dimension of the truncation variety is
well understood. The reason is that cluster analysis allows us to bijectively
map cluster amplitudes to CI-coefficients using polynomial equations, see e.g.
Refs. 12, 70, 32 or Proposition 2.4 in Ref. 60 for a mathematical proof. Then,
since $V_{\sigma}$ arises as a deformation of the linear space
$\mathbb{V}_{\sigma}$ under a polynomial with polynomial inverse, the
dimension of $V_{\sigma}$ is equal to the dimension of $\mathbb{V}_{\sigma}$.
For example in the case of CCSD, we obtain

$\dim(V_{\sigma})=\dim(\mathbb{V}_{\sigma})=n_{\rm virt}n_{\rm
occ}+\binom{n_{\rm virt}}{2}\binom{n_{\rm occ}}{2}.$ (13)

The degree of ${V}_{\sigma}$ on the other hand is not as easily computable and
requires a more profound understanding of ${V}_{\sigma}$. For special cases,
however, explicit formulas for the degree have been established; for example,
for CCS, the truncation variety is the Grassmanian (see Theorem 3.5 in Ref.
60) which is mathematically well understood. In other cases we rely either on
symbolic methods, provided by e.g. Macaulay2 71, or on numerical methods,
accessible in e.g. HomotopyContinuation.jl 72, to compute the degree.

Despite the complexity of calculating $\deg(V_{\sigma})$ the bound given by
Eq. (12) provides significant insight into bounding the CC degree, since it is
less complex to calculate $\deg(V_{\sigma})$ than to calculate the CC degree
directly. Moreover, this upper bound is notably better than the previously
known upper bounds obtained by applying Bézout’s theorem directly to the
projective CC equations, see Section 6 in Ref. 60.

##  4 Computational methods

In this section, we outline the computational procedure employed to compute
the CC degree and the CC roots for a given CC truncation rank
$\sigma\subset[\\![N]\\!]$. We begin by introducing the CC family, a
parametric polynomial family for Eq. (10), i.e.

$\mathcal{F}_{\rm CC}(\sigma)=\\{[(\hat{G}-E\cdot
I_{K})e^{\hat{T}(\mathbf{t})}]_{\sigma}~{}\big{|}~{}\hat{G}\in\mathbb{F}^{K\times
K}\\},$ (14)

where the energy $E\in\mathbb{C}$ and cluster amplitudes
$\mathbf{t}\in\mathbb{V}_{\sigma}$ are the variables, and the entries of the
matrix $\hat{G}$ are the parameters. We will ultimately set $\hat{G}=\hat{H}$
to be the Hamiltonian of interest, but $\mathcal{F}_{\rm CC}(\sigma)$ is a
much broader family of CC-like systems with $\hat{G}\in\mathbb{F}^{K\times K}$
an arbitrary matrix. The parameter continuation theorem provides an upper
bound to the number of roots of any polynomial system
$F_{\hat{G}}\in\mathcal{F}_{\rm CC}$. Moreover, this bound is tight and
equality is obtained for all polynomial systems with parameters in a certain
dense subspace of $\mathbb{F}^{K\times K}$. We subsequently assume the
validity of the parameter continuation theorem and refer the interested reader
to Section 3 in Ref. 73. The upper bound provided by the parameter
continuation theorem yields the CC degree mentioned in the previous section.

The numerical procedure used in this work consists of two steps: First, we
solve a general system of equations from the CC family, $\mathcal{F}_{\rm
CC}$, using monodromy methods. Second, we use parameter homotopy continuation
to track the solutions from this solved but generic system of equations to
solutions of the CC equations describing a specific chemical system of
interest. Both steps are subsequently outlined in detail.

We emphasize that, at their core, both methods use a homotopy to track new
solutions from known solutions 74, 75, 76, 77. More precisely, we define a
homotopy for a (piecewise smooth) path $\gamma:[0,1]\to\mathbb{F}^{K\times K}$
in the parameter space via

$H(E,\mathbf{t},\lambda)=F_{\gamma(\lambda)}(E,\mathbf{t}),$ (15)

where $H(E,\mathbf{t},\lambda)$ describes the continuous deformation of a
start system $G(E,\mathbf{t})=H(E,\mathbf{t},1)$ to a target system
$F(E,\mathbf{t})=H(E,\mathbf{t},0)$. The individual solution paths to the
system $H(E,\mathbf{t},\lambda)$ are then tracked from $G(E,\mathbf{t})$ to
$F(E,\mathbf{t})$ as $\lambda\to 0$. This is accomplished via an ordinary
differential equation, known as the Davidenko differential equation 78, 79,
i.e.

$\frac{\partial}{\partial\mathbf{x}}H(\mathbf{x},\lambda)\left(\frac{\mathrm{d}}{\mathrm{d}\lambda}\mathbf{x}(\lambda)\right)+\frac{\partial}{\partial\lambda}H(\mathbf{x},\lambda)=0,$
(16)

initialized by a root $\mathbf{x}$ of the start system, i.e. $G(\mathbf{x})=0$
where we set
$\mathbf{x}=(E,\mathbf{t})\in\mathbb{C}\times\mathbb{V}_{\sigma}$.

###  4.1 Monodromy

The monodromy solver is a numerical technique specifically designed for
solving generic polynomial systems within a parameterized family. We emphasize
that the monodromy solver requires the system to be generic, meaning it is not
suitable for solving a specific set of CC equations directly. However, it can
be employed to establish a start system $F_{\hat{G}}\in\mathcal{F}_{\rm CC}$
within the CC family. This start system is then used by the parameter homotopy
continuation method to solve a (physical) target system of interest. Note that
the monodromy method is required only once for a particular system
configuration, i.e. the number of electrons, the number of orbitals, and the
CC truncation level $\sigma\in[\\![N]\\!]$.

We start with a root $\mathbf{x}\in\mathbb{C}\times\mathbb{V}_{\sigma}$ of a
generic system $F_{\hat{G}}\in\mathcal{F}_{\rm CC}$ in the CC family and a
loop $\gamma:[0,1]\to\mathbb{F}^{K\times K}$ in the parameter space, based at
the generic matrix $\hat{G}$, i.e. $\gamma(0)=\gamma(1)=\hat{G}$. Moreover, we
let $\mathbf{x}(\lambda)$ be the solution path of $H(E,\mathbf{t},\lambda)=0$
along $\gamma$, starting at $\mathbf{x}(1)=\mathbf{x}$. Since $\gamma$ is a
loop, the start system in this homotopy is equal to the end system. In
particular, the solution path ends again at a root of $F_{\hat{G}}$. The
monodromy solver exploits the fact that the root $\mathbf{x}(0)$ need not
necessarily be equal to the root $\mathbf{x}(1)$. In fact, let
$(\mathbf{x}_{1},\dots,\mathbf{x}_{M})$ be all roots to $F_{\hat{G}}$, then
the roots $(\mathbf{x}_{1}(0),\dots,\mathbf{x}_{M}(0))$ define a mere
permutation the system’s roots.

We exemplify this procedure with the simple parametric polynomial family:

$\mathcal{F}=\\{x^{3}-6\,x^{2}+11\,z\,x-6~{}\big{|}~{}z\in\mathbb{C}\\},$ (17)

where $x\in\mathbb{C}$ is the variable and $z\in\mathbb{C}$ is the parameter.
For $z=1$, the corresponding polynomial

$x^{3}-6\,x^{2}+11\,x-6=(x-1)(x-2)(x-3)$ (18)

has the roots $x_{1}=1,x_{2}=2,x_{3}=3$. The monodromy solver then tracks
these roots along the path $\gamma:[0,1]\to\mathbb{C},t\mapsto e^{i2\pi t}$.
Note that $\gamma$ is a loop based at $1$. Indeed, this procedure permutes the
roots, see Figure 2.

Figure 2: The solution paths of the roots of $x^{3}-6x^{2}+11zx-6$ for
$z(\lambda)=\gamma(\lambda)$ along the path
$\gamma:[0,1]\to\mathbb{C},t\mapsto e^{i2\pi t}$.

This approach can be effectively utilized to find all solutions to a generic
system of equations. For instance, consider a scenario where the only known
solution for the polynomial in Eq. (18) is $x_{2}=2$. By tracking this
solution along the unit circle, we can determine another solution, $x_{1}=1$.
When we extend this tracking across sufficiently many loops, the method
reveals all possible roots of the system.

This process is central to the monodromy solver. The solver begins with a
single known solution and numerically traces known solutions through multiple
loops within the parameter space. This tracking continues until no additional
solutions are found over a predetermined number of loops, leading to the
conclusion that all possible solutions to the system have been numerically
identified.

The parameters in the CC family $\mathcal{F}_{\rm CC}$ are linear, hence, we
can find a generic start system along with a root by choosing a pair
$(E,\mathbf{t})\in\mathbb{C}\times\mathbb{V}_{\sigma}$ at random, and then
compute one generic solution $\hat{G}_{0}$ fulfilling

$[(\hat{G}_{0}-E\cdot I_{K})e^{\hat{T}(\mathbf{t})}]_{\sigma}=0.$

We then use the monodromy solver to compute all roots to $F_{\hat{G}_{0}}$.

###  4.2 Parametric homotopy continuation

Applying the parametric homotopy continuation method to find the roots of the
CC equations corresponding to a system of interest is now straightforward. The
goal is to solve

$F_{\hat{H}}(E,{\mathbf{t}})=[(\hat{H}-E\cdot
I_{K})e^{\hat{T}(\mathbf{t})}]_{\sigma}=0$ (19)

for a given Hamiltonian $\hat{H}$, arising from a physical electronic system.
By construction, both $\hat{H}$ and $\hat{G}_{0}$ are in the parameter space
of $\mathcal{F}_{\rm CC}$. In particular, we can define a path $\gamma$ in the
parameter space that continuously transforms $\hat{G}_{0}$ into $\hat{H}$.
Since we have computed all roots to the system $F_{\hat{G}_{0}}$, the
Davidenko differential equation (16) allows the tracking of solutions of
$F_{\hat{G}_{0}}$ to $F_{\hat{H}}$ using e.g. Newton’s method.

As $\lambda\to 0$, we may encounter different scenarios when getting close to
$\lambda=0$. Here, the endgame 80 should be employed to determine the target
solutions and their “type”. We illustrate such different scenarios in Figure
3.

Figure 3: Sketch of possible homotopy paths. The solid line shows a path with
no finite limit as $\lambda\to 0$, the dashed lines have the same limit of a
singular root, and the dotted-dashed line has a unique limit of a nonsingular
root.

We conclude this section by highlighting the progress that has been made
within homotopy continuation based software development exploiting parallel
implementations, which can significantly extend their application areas in the
coming years. In particular, we here want to emphasize PHCpack 81, Bertini, 82
HOM4PS, 83, 84 NAG4M2 74 and HomotopyContinuation.jl. 72

##  5 Numerical results

In this section we investigate the root structure and potential energy curves
(PECs) for ground and excited states at the level of CCD and CCSD. We report
PECs corresponding to different CC solutions obtained by the above-outlined
computational procedure. The obtained energies and corresponding eigenstate
approximations are then compared to exact eigenpairs of the respective
Hamiltonian. Besides the computation and analysis of the CC states, we compare
different bounds to the number of CC roots with the computed number of CC
solutions along different bond stretching processes. We restrict our
investigations to four electron systems, such as (H2)2 in D2h and D∞h
configurations, H4 symmetrically disturbed on a circle 85, and LiH. The
reported computations are performed in first quantization and the considered
systems are all at half-filling, i.e. four electrons in eight spin orbitals.
The one-body and two-body integrals for the considered systems are obtained
from the Python-based Simulations of Chemistry Framework (PySCF) 86, 87, 88.
For the algebraic simulations, we used the monodromy and parametric homotopy
continuation methods in HomotopyContinuation.jl. 72

###  5.1 CC bounds for generic four electron systems at half-filling

We emphasize that the theoretical bounds hold for generic Hamiltonians, and
therefore only depend on the number of electrons, the number of spin orbitals,
and the CC truncation level.

##### Coupled cluster doubles:

The roughest bound is obtained by applying Bézout’s theorem directly to the
CCD equations. Note that for four electron systems, the CCD equations are
quadratic. Moreover, in first quantization, there are ${4\choose 2}^{2}=36$
doubly excited Slater determinants. Hence, the Bézout bound directly applied
to the CCD equations yields $2^{36}$ as an upper bound for the number of
roots. The second bound we compare with is the bound in Eq. (12). Again, in
first quantization, there are 36 CCD amplitudes, hence, ${\rm
dim}(V_{\\{2\\}})=36$. Moreover we compute the degree of $V_{\\{2\\}}$ using
Macaulay2 yielding ${\rm deg}(V_{\\{2\\}})=2$. Overall, this yields $74$ as an
upper bound. For a generic Hamiltonian, we can moreover compute the exact
number of roots to the CCD equations using the monodromy method which yields
CCdeg${}_{4,8}(\\{2\\})=73$. This quantifies the improvement provided by the
bound in Eq. (12) over the direct application of Bézout’s theorem.

##### Coupled cluster singles and doubles:

When considering CCSD instead of CCD these numbers dramatically increase.
Frist, we note that the Bézout bound directly applied to the CCSD equation
yields the upper bound $3^{n_{s}}4^{n_{d}}$, where for four electrons in eight
spin orbitals $n_{s}=16$ and $n_{d}=36$. This leads to approximately $2\cdot
10^{29}$ solutions. Second, we note that ${\rm dim}(V_{\\{1,2\\}})=52$ and
numerical calculations yield ${\rm deg}(V_{\\{1,2\\}})=442066$, hence, the
bound provided by Eq. (12) yields approximately $2\cdot 10^{7}$ solutions.
Numerical calculations show that CCdeg${}_{4,8}(\\{1,2\\})=16\,952\,996$. This
clearly illustrates the improvement provided by the bound in Eq. (12) over the
direct application of Bézout’s theorem.

To connect these theoretical bounds to practical systems, we subsequently
calculate the number of CC solutions that result in real-valued energies, as
well as those that accurately predict excited state energies at different bond
lengths for the considered systems. Additionally, we report the number of FCI
solutions. This comparison allows us to assess the accuracy of the bounds
established for generic Hamiltonians.

###  5.2 Dissociation of lithium hydride – CCD

We initiate our energy and state approximation investigations by looking at
the dissociation process of lithium hydride. We consider bond distances
ranging from $1.375$ to $5.95$ bohr. Diagonalizing the Hamiltonian yields the
full spectrum reported in Figure 16(a) (see supplementary material) together
with all roots to the CCD equations. Due to the large number of FCI solutions,
we minimize the distance between FCI solutions and CCD solutions showing a
much clearer picture of potentially well-resolved PECs, see Figure 16(b) in
the supplementary material. This graph allows us to extract the PECs of the
different states that are well-approximated using CCD, see Figure 4(a). We
find that 10 PECs are approximated up to $5\cdot 10^{-2}$ hartree for the
entire dissociation process, however, we also note that for various bond
distances, more CCD energies approximate FCI energies well, see Figure 16(b)
and Table 1. For the 10 PECs, we moreover compute the overlap between the CCD
states and the approximated eigenstates, see Figure 4(b). Remarkably, three of
the four lowest energy states are well approximated both in terms of energy
and states. The highest energy state stands out as well since it is well-
approximated near the equilibrium, both in terms of energy as well as the
eigenstate. However, as we transition away from the equilibrium the level of
approximation deteriorates. Note that poor overlap may be detected when the
eigenspace is higher dimensional, that is if the corresponding eigenvalue has
multiplicity two or higher. Therefore some of the higher energy states
reporting bad overlap might have CCD states that are close to the actual
eigenspace – or even within the eigenspace – since we are only reporting
overlap with one representative from the eigenspace. A detailed investigation
of this is left for future works.

(a)

(b)

Figure 4: (a) PECs of LiH that can be accurately approximated by CCD energies.
(b) Overlaps of the CCD states with the corresponding eigenstate.

Investigating the root structure, we observe that the total number of CCD
roots fluctuates along the bond stretching process, see Table 1. We find that
any bound for a generic Hamiltonian severely overestimates the true number of
roots for this system.

Bond distances | 1.375 | 1.6 | 1.9 | 2.2 | 2.5 | 3.1 | 4.0 | 4.9 | 5.8
---|---|---|---|---|---|---|---|---|---
$\\#$ CCD real | 33 | 40 | 33 | 40 | 34 | 28 | 30 | 33 | 28
$\\#$ CCD approx. | 22 | 25 | 20 | 23 | 18 | 19 | 17 | 18 | 14
Table 1: The number of CCD roots in LiH dissociation for selected bond
distances. By "$\\#$ CCD approx." we denote the energetically relevant CCD
roots. For comparison, we recall that the Bézout bound led to $2^{36}$
solutions, the bound in Eq. (12) led to 74 solutions, and CCdeg4,8({2}) is 73.

###  5.3 Lithium hydride – CCSD

We moreover perform CCSD computations on lithium hydride at the bond distance
2.875 bohr – close to the molecule’s equilibrium. Computing all roots to CCSD
equations yields $15\,954$ roots, $2\,170$ of which yield real-values
energies, and $1\,280$ are real-valued solutions. Visualizing all $15\,954$
CCSD energies shows that the energies tend to cluster around the FCI
solutions, see Figure 5(a). Being able to compare with the true Hamiltonian
spectrum reveals that only 26 of these CCSD solutions yield energies that are
close to an FCI energy up to $10^{-3}$ hartree, see Figure 5(b).

(a)

(b)

Figure 5: (a) FCI solutions (red lines) together with all CCSD energies (b)
CCSD energies that approximate FCI energies up to $10^{-3}$ hartree.

For all 70 eigenstates, we compute the energetically closest CCSD state
together with the overlap with the corresponding states, see Figure 6. We find
that three of the 26 energetically relevant CCSD states have good overlaps
with the targeted states. Therefore there are at least three CCSD states that
are well approximated both in terms of energy and state.

Figure 6: Overlap of the energetically closest CCSD state with the
corresponding eigenstate. The blue dots show the 26 energetically relevant
CCSD states.

###  5.4 Dissociating systems of hydrogen

We now investigate the dissociation of (H2)2 planar model systems in different
geometries 38.

####  5.4.1 (H2)2 dissociation (D2h symmetry)

The bond stretching procedure for (H2)2 in D2h symmetry is sketched in Figure
7. We consider an intra-molecular bond distance of $R=1.4$ bohr, which
corresponds to the equilibrium geometry of H2; this is kept fixed during the
dissociation process. The inter-molecular distance is varied from 1.5 to 4.0
bohr.

Figure 7: Schematic depiction of the dissociation process of (H2)2 in D2h
configuration. The parameter $R=1.4$ bohr.

The full spectrum together with the CCD solutions can be found in Figure 13(a)
in the Supplementary material, and the PECs together with the closest CCD
solutions are reported in Figure 13(b). We extract the PECs that are well-
approximated using CCD and report them in Figure 8(a). Here, we observe that
eight CCD energies approximate FCI PECs well. For these eight CCD states, we
compute the overlap with the corresponding eigenstate, see Figure 8(b). We see
that the ground state is well approximated, both in terms of energy as well as
the state. The following six states (i.e. states two to seven in Figure 8(b))
have intermediate overlaps with the targeted states; the last state has a poor
overlap with the targeted state. This illustrates that CC amplitudes could
yield good energy approximations while providing poor approximations to the
actual eigenstates 13.

(a)

(b)

Figure 8: (a) PECs of (H2)2 in D2h symmetry that are accurately described by
CCD energies. (b) Overlap of the CCD states with the corresponding eigenstate.

In Table 2 we report the number of CCD roots yielding real-valued energies and
energetic relevant CCD roots for selected bond distances. We observe that the
number of CCD roots fluctuates along the bond stretching procedure.

Bond distances | 1.6 | 1.9 | 2.1 | 2.5 | 3.0 | 3.9
---|---|---|---|---|---|---
$\\#$ CCD real | 66 | 45 | 60 | 64 | 55 | 64
$\\#$ CCD approx. | 32 | 26 | 33 | 32 | 33 | 34
Table 2: The number of roots along the (H2)2 dissociation in D2h configuration
for selected bond distances. By "$\\#$ CCD approx." we denote the
energetically relevant CCD roots. For comparison, we recall that the Bézout
bound is $2^{36}$, the bound in Eq. (12) yields 74, and CCdeg4,8({2}) is 73.

####  5.4.2 (H2)2 dissociation (D∞h symmetry)

Next, we investigate the (H2)2 dissociation in D∞h symmetry, see Figure 9 for
a sketch of the dissociation process. The intra-molecular distance is again
set to $R=1.4$ bohr and is kept fixed during the dissociation process. The
inter-molecular distance is varied from 1.5 to 4.0 bohr.

Figure 9: Schematic depiction of the dissociation process of (H2)2 in D∞h
configuration. The parameter $R=1.4$ bohr.

The full spectrum and CCD solutions is presented in Figure 14(a) in the
Supplementary material, with the PECs and their closest CCD solutions detailed
in Figure 14(b). We have identified PECs that are accurately approximated by
CCD, as shown in Figure 10(a), where four CCD energies closely match the FCI
PECs. The overlap of these four CCD states with their respective eigenstates
is analyzed in Figure 10(b). Consistent with observations from the (H2)2 in
D2h configuration, the ground state demonstrates a high degree of accuracy,
while the other three states exhibit intermediate to poor overlap with their
targeted eigenstates.

(a)

(b)

Figure 10: (a) PECs of (H2)2 in D∞h symmetry that are accurately described by
CCD energies. (b) Overlap of the CCD states with the corresponding eigenstate.

Investigating the root structure, we observe that – similar to (H2)2 in D2h
symmetry – the total number of CCD roots fluctuates along the bond stretching
procedure, see Table 3. Note that for (H2)2 in D∞h symmetry, the total number
of CCD roots is lower than the number of CCD roots (H2)2 in D2h symmetry. In
particular, any bound for a generic Hamiltonian severely overestimates the
true number of roots for this system.

Bond distances | 1.6 | 1.9 | 2.1 | 2.5 | 3.0 | 3.9
---|---|---|---|---|---|---
$\\#$ CCD real | 46 | 42 | 40 | 38 | 37 | 35
$\\#$ CCD approx. | 13 | 19 | 12 | 14 | 16 | 21
Table 3: The number of roots along the (H2)2 dissociation in D∞h configuration
for selected bond distances. By "$\\#$ CCD approx." we denote the
energetically relevant CCD roots. For comparison, we recall that the Bézout
bound led to $2^{36}$ solutions, the bound in Eq. (12) led to 74 solutions,
and CCdeg4,8({2}) is 73.

###  5.5 H4 disturbed on a circle

We proceed by investigating a variant of the H4 model consisting of four
hydrogen atoms symmetrically distributed on a circle of radius $R=\sqrt{2}$
bohr 89, see Fig. 11.

Figure 11: Schematic depiction of the H4 model undergoing a symmetric
disturbance on a circle modeled by the angle $\Theta$.

The full spectrum together with the CCD solutions can be found in Figure 15(a)
in the Supplementary material, and the PECs together with the closest CCD
solutions are reported in Figure 15(b). We extract the PECs that are well-
approximated using CCD which are reported in Figure 12(a). Note that it is
well-known that for this system CCD does not perform well in the region close
to 90∘, due to strong degeneracies 38. This results in a dramatic reduction of
the total number of roots, we find three roots at $90^{\circ}$ and six roots
at $89^{\circ}$. For completion, these energies are included in the
supplementary material, however, since the focus of this work is on the root
structure and potentially physically relevant roots accessible by single-
reference CC theory, we focus on cases where single-reference CC theory
provides reasonable approximations and therefore excluded these points in
Figure 12(a). Outside of this challenging region, we find five solutions that
approximate PECs well. For these five solutions, we compute the overlap with
the corresponding eigenstate, see Figure 12(b). Consistent with the previous
H4 systems, this model system shows the effect of CCD states accurately
resolving excited energies but providing poor overlap to the exact
eigenstates. The only state that shows a good approximation to its eigenstate
is the CCD ground state.

(a)

(b)

Figure 12: (a) PECs of H4 symmetrically distributed on a cycle that can be
accurately reconstructed by CCD roots. (b) Overlap of the CCD states with the
corresponding eigenstate.

Similar to (H2)2 in D∞h symmetry, we observe that generic bounds severely
overestimate the number of CCD roots, though this effect is further amplified
for H4 symmetrically distributed on a circle, see Table 4.

Angles in degrees | 45 | 55 | 65 | 75 | 85
---|---|---|---|---|---
$\\#$ CCD real | 35 | 36 | 35 | 36 | 28
$\\#$ CCD approx. | 18 | 18 | 19 | 21 | 22
Table 4: The number of roots of H4 symmetrically distributed on a circle for
selected bond distances. By "$\\#$ CCD approx." we denote the energetically
relevant CCD roots. For comparison, we recall that the Bézout bound led to
$2^{36}$ solutions, the bound in Eq. (12) led to 74 solutions, and
CCdeg4,8({2}) is 73.

##  6 Conclusion

The introduction highlighted critical unanswered questions in the study of the
coupled cluster (CC) equations’ algebraic structures, including the total
number of solutions, the best approximation element, and the feasibility of
approximating excited states with higher-order roots of single reference CC
equations. The ability to compute the full solution spectrum to the CC
equations enabled us to tackle these questions for concrete physical systems.

Most remarkably, our results demonstrate that for lithium hydride in multiple
geometries, higher-order CC solutions not only provide high-accuracy
approximations to the energies but also to the states themselves. For a
broader range of systems including multiple (H2)2 planar model systems in
different geometries, we found that high-order CC solutions accurately predict
several full potential energy curves, though they possibly provide a poor
approximation to the eigenstates. While this effect has been previously
reported, it requires further investigation to fully understand it and
ultimately leverage it for practical applications.

Additionally, we examined the total number of roots and the effectiveness of
various bounds developed for generic Hamiltonians. Our simulations indicate
that while the most recent bounds perform exceptionally well for generic
systems, they significantly overestimate the number of solutions in physical
systems. More precisely, for all considered systems, we computed between
$30$-$60$ real solutions of which only $10$-$30$ approximate physical states,
compared to CC degree of $72$. This excess of solutions gets amplified for
CCSD: Solving the CC equations of a generic Hamiltonian yields $16\,952\,996$
solutions, compared to $1280$ solutions that yield real-valued energies of the
CC equations describing LiH. This overestimation by a factor of over ten
thousand indicates that to further improve the bounds of the CC equations,
specific structures of the Hamiltonian have to be taken into account.

The authors thank Piotr Piecuch for insightful discussions.

## References

  * Bartlett and Musiał 2007 Bartlett, R. J.; Musiał, M. Coupled-cluster theory in quantum chemistry. _Reviews of Modern Physics_ 2007, _79_ , 291
  * Coester 1958 Coester, F. Bound States of a Many-Particle System. _Nuclear Physics_ 1958, _7_ , 421–424
  * Hubbard 1957 Hubbard, J. The description of collective motions in terms of many-body perturbation theory. _Proceedings of the Royal Society A_ 1957, _240_ , 539–560
  * Hugenholtz 1957 Hugenholtz, N. Perturbation approach to the Fermi gas model of heavy nuclei. _Physica_ 1957, _23_ , 533–545
  * Čížek 1966 Čížek, J. On the Correlation Problem in Atomic and Molecular Systems. calculation of wavefunction components in Ursell-type expansion using quantum-field theoretical methods. _The Journal of Chemical Physics_ 1966, _45_ , 4256–4266
  * Kümmel 1991 Kümmel, H. Origins of the coupled cluster method. _Theoretica Chimica Acta_ 1991, _80_ , 81–89
  * Čížek 1991 Čížek, J. Origins of coupled cluster technique for atoms and molecules. _Theoretica Chimica Acta_ 1991, _80_ , 91–94
  * Bartlett 2005 Bartlett, R. In _Theory and Applications of Computational Chemistry (The First Forty Years)_ , 1st ed.; Dykstra, C., Frenking, G., Kim, K., Scuseria, G., Eds.; 2005; pp 1191–1221
  * Paldus 2005 Paldus, J. In _Theory and Applications of Computational Chemistry (The First Forty Years)_ , 1st ed.; Dykstra, C., Frenking, G., Kim, K., Scuseria, G., Eds.; 2005; pp 115–147
  * Arponen 1991 Arponen, J. S. Independent-cluster methods as mappings of quantum theory into classical mechanics. _Theoretica Chimica Acta_ 1991, _80_ , 149–179
  * Bishop 1991 Bishop, R. An overview of coupled cluster theory and its applications in Physics. _Theoretica Chimica Acta_ 1991, _80_ , 95–148
  * Helgaker et al. 2014 Helgaker, T.; Jörgensen, P.; Olsen, J. _Molecular Electronic-Structure Theory_ ; John Wiley & Sons, 2014
  * Piecuch and Kowalski 2000 Piecuch, P.; Kowalski, K. In _Computational Chemistry: Reviews of Current Trends_ ; Leszczynski, J., Ed.; Singapore: World Scientific, 2000; Vol. 5
  * Rowe 1968 Rowe, D. Equations-of-motion method and the extended shell model. _Reviews of Modern Physics_ 1968, _40_ , 153
  * Emrich 1981 Emrich, K. An extension of the coupled cluster formalism to excited states (I). _Nuclear Physics A_ 1981, _351_ , 379–396
  * Sekino and Bartlett 1984 Sekino, H.; Bartlett, R. J. A linear response, coupled-cluster theory for excitation energy. _International Journal of Quantum Chemistry_ 1984, _26_ , 255–265
  * Geertsen et al. 1989 Geertsen, J.; Rittby, M.; Bartlett, R. J. The equation-of-motion coupled-cluster method: Excitation energies of Be and CO. _Chemical Physics Letters_ 1989, _164_ , 57–62
  * Stanton and Bartlett 1993 Stanton, J. F.; Bartlett, R. J. The equation of motion coupled-cluster method. A systematic biorthogonal approach to molecular excitation energies, transition probabilities, and excited state properties. _The Journal of chemical physics_ 1993, _98_ , 7029–7039
  * Comeau and Bartlett 1993 Comeau, D. C.; Bartlett, R. J. The equation-of-motion coupled-cluster method. Applications to open-and closed-shell reference states. _Chemical Physics Letters_ 1993, _207_ , 414–423
  * Watts and Bartlett 1994 Watts, J. D.; Bartlett, R. J. The inclusion of connected triple excitations in the equation-of-motion coupled-cluster method. _The Journal of chemical physics_ 1994, _101_ , 3073–3078
  * Monkhorst 1977 Monkhorst, H. J. Calculation of properties with the coupled-cluster method. _International Journal of Quantum Chemistry_ 1977, _12_ , 421–432
  * Dalgaard and Monkhorst 1983 Dalgaard, E.; Monkhorst, H. J. Some aspects of the time-dependent coupled-cluster approach to dynamic response functions. _Physical Review A_ 1983, _28_ , 1217
  * Koch and Jørgensen 1990 Koch, H.; Jørgensen, P. Coupled cluster response functions. _The Journal of chemical physics_ 1990, _93_ , 3333
  * Watson and Chan 2012 Watson, M. A.; Chan, G. K.-L. Excited states of butadiene to chemical accuracy: Reconciling theory and experiment. _Journal of chemical theory and computation_ 2012, _8_ , 4013–4018
  * Shu and Truhlar 2017 Shu, Y.; Truhlar, D. G. Doubly excited character or static correlation of the reference state in the controversial 21Ag state of trans-butadiene. _Journal of the American Chemical Society_ 2017, _139_ , 13770–13778
  * Barca et al. 2018 Barca, G. M.; Gilbert, A. T.; Gill, P. M. Excitation number: Characterizing multiply excited states. _Journal of chemical theory and computation_ 2018, _14_ , 9–13
  * Loos et al. 2019 Loos, P.-F.; Boggio-Pasqua, M.; Scemama, A.; Caffarel, M.; Jacquemin, D. Reference energies for double excitations. _Journal of chemical theory and computation_ 2019, _15_ , 1939–1956
  * Ravi et al. 2022 Ravi, M.; Perera, A.; Bartlett, R. J., et al. The intermediate state approach for doubly excited dark states in EOM-coupled-cluster theory. _The Journal of Chemical Physics_ 2022, _156_
  * Do Casal et al. 2023 Do Casal, M. T.; Toldo, J. M.; Barbatti, M.; Plasser, F. Classification of doubly excited molecular electronic states. _Chemical Science_ 2023, _14_ , 4012–4026
  * Piecuch and Adamowicz 1994 Piecuch, P.; Adamowicz, L. Solving the single-reference coupled-cluster equations involving highly excited clusters in quasidegenerate situations. _The Journal of chemical physics_ 1994, _100_ , 5857–5869
  * Mihálka and Noga 2023 Mihálka, Z. É.; Noga, J. Exploring alternative approaches to improve the convergence pattern in solving the coupled-cluster equations. _Molecular Physics_ 2023, _121_ , e2140084
  * Shavitt and Bartlett 2009 Shavitt, I.; Bartlett, R. J. _Many-body methods in chemistry and physics: MBPT and coupled-cluster theory_ ; Cambridge university press, 2009
  * Szakács and Surján 2008 Szakács, P.; Surján, P. R. Stability conditions for the coupled cluster equations. _International Journal of Quantum Chemistry_ 2008, _108_ , 2043–2052
  * Adams et al. 1981 Adams, B.; Jankowski, K.; Paldus, J. Symmetry-adapted coupled-pair approach to the many-electron correlation problem. II. Application to the Be atom. _Physical Review A_ 1981, _24_ , 2316
  * Faulstich 2024 Faulstich, F. M. Recent mathematical advances in coupled cluster theory. _arXiv:2401.07383_ 2024,
  * Faulstich et al. 2024 Faulstich, F. M.; Khoo, Y.; Li, K. Augmented Lagrangian method for coupled-cluster. _arXiv:2403.16381_ 2024,
  * Živković and Monkhorst 1978 Živković, T. P.; Monkhorst, H. J. Analytic connection between configuration–interaction and coupled-cluster solutions. _Journal of Mathematical Physics_ 1978, _19_ , 1007–1022
  * Paldus et al. 1993 Paldus, J.; Piecuch, P.; Pylypow, L.; Jeziorski, B. Application of Hilbert-space coupled-cluster theory to simple (H2)2 model systems: Planar models. _Physical Review A_ 1993, _47_ , 2738
  * Piecuch et al. 1990 Piecuch, P.; Zarrabian, S.; Paldus, J.; Čížek, J. Coupled-cluster approaches with an approximate account of triexcitations and the optimized-inner-projection technique. II. Coupled-cluster results for cyclic-polyene model systems. _Physical Review B_ 1990, _42_ , 3351
  * Kowalski and Jankowski 1998 Kowalski, K.; Jankowski, K. Towards complete solutions to systems of nonlinear equations of many-electron theories. _Physical Review Letters_ 1998, _81_ , 1195
  * Piecuch and Kowalski 2000 Piecuch, P.; Kowalski, K. In _Computational Chemistry, Reviews of Current Trends_ ; Leszczynski, J., Ed.; 2000; pp 1–104
  * Kowalski and Piecuch 2000 Kowalski, K.; Piecuch, P. Complete set of solutions of multireference coupled-cluster equations: The state-universal formalism. _Physical Review A_ 2000, _61_ , 052506
  * Csirik, Mihály A. and Laestadius, Andre 2023 Csirik, Mihály A.,; Laestadius, Andre, Coupled-Cluster theory revisited - Part I: Discretization. _ESAIM: Mathematical Modelling and Numerical Analysis_ 2023, _57_ , 645–670
  * Csirik and Laestadius 2023 Csirik, M. A.; Laestadius, A. Coupled-Cluster theory revisited-Part II: Analysis of the single-reference Coupled-Cluster equations. _ESAIM: Mathematical Modelling and Numerical Analysis_ 2023, _57_ , 545–583
  * Kowalski and Piecuch 2000 Kowalski, K.; Piecuch, P. Complete set of solutions of the generalized Bloch equation. _International Journal of Quantum Chemistry_ 2000, _80_ , 757–781
  * Kowalski and Jankowski 1998 Kowalski, K.; Jankowski, K. Full solution to the coupled-cluster equations: the H4 model. _Chemical Physics Letters_ 1998, _290_ , 180–188
  * Jankowski and Kowalski 1999 Jankowski, K.; Kowalski, K. Physical and mathematical content of coupled–cluster equations: Correspondence between coupled–cluster and configuration–interaction solutions. _The Journal of Chemical Physics_ 1999, _110_ , 3714–3729
  * Jankowski and Kowalski 1999 Jankowski, K.; Kowalski, K. Physical and mathematical content of coupled-cluster equations. II. On the origin of irregular solutions and their elimination via symmetry adaptation. _The Journal of Chemical Physics_ 1999, _110_ , 9345–9352
  * Jankowski and Kowalski 1999 Jankowski, K.; Kowalski, K. Physical and mathematical content of coupled-cluster equations. III. Model studies of dissociation processes for various reference states. _The Journal of Chemical Physics_ 1999, _111_ , 2940–2951
  * Jankowski and Kowalski 1999 Jankowski, K.; Kowalski, K. Physical and mathematical content of coupled-cluster equations. IV. Impact of approximations to the cluster operator on the structure of solutions. _The Journal of Chemical Physics_ 1999, _111_ , 2952–2959
  * Kowalski and Piecuch 2000 Kowalski, K.; Piecuch, P. The method of moments of coupled-cluster equations and the renormalized CCSD[T], CCSD(T), CCSD(TQ), and CCSDT(Q) approaches. _The Journal of Chemical Physics_ 2000, _113_ , 18–35
  * Piecuch et al. 2002 Piecuch, P.; Kowalski, K.; Pimienta, I. S.; Mcguire, M. J. Recent advances in electronic structure theory: Method of moments of coupled-cluster equations and renormalized coupled-cluster approaches. _International Reviews in Physical Chemistry_ 2002, _21_ , 527–655
  * Piecuch et al. 2004 Piecuch, P.; Kowalski, K.; Pimienta, I.; Fan, P.-D.; Lodriguito, M.; McGuire, M.; Kucharski, S.; Kuś, T.; Musiał, M. Method of moments of coupled-cluster equations: a new formalism for designing accurate electronic structure methods for ground and excited states. _Theoretical Chemistry Accounts_ 2004, _112_ , 349–393
  * Shen and Piecuch 2012 Shen, J.; Piecuch, P. Biorthogonal moment expansions in coupled-cluster theory: Review of key concepts and merging the renormalized and active-space coupled-cluster methods. _Chemical Physics_ 2012, _401_ , 180–202
  * Kowalski and Piecuch 2000 Kowalski, K.; Piecuch, P. Renormalized CCSD(T) and CCSD(TQ) approaches: Dissociation of the N2 triple bond. _The Journal of Chemical Physics_ 2000, _113_ , 5644–5652
  * Piecuch et al. 2006 Piecuch, P.; Włoch, M.; Gour, J. R.; Kinal, A. Single-reference, size-extensive, non-iterative coupled-cluster approaches to bond breaking and biradicals. _Chemical physics letters_ 2006, _418_ , 467–474
  * Piecuch and Włoch 2005 Piecuch, P.; Włoch, M. Renormalized coupled-cluster methods exploiting left eigenstates of the similarity-transformed Hamiltonian. _The Journal of chemical physics_ 2005, _123_
  * Faulstich and Laestadius 2023 Faulstich, F. M.; Laestadius, A. Homotopy continuation methods for coupled-cluster theory in quantum chemistry. _Molecular Physics_ 2023, e2258599
  * Faulstich and Oster 2024 Faulstich, F. M.; Oster, M. Coupled Cluster Theory: Toward an Algebraic Geometry Formulation. _SIAM Journal on Applied Algebra and Geometry_ 2024, _8_ , 138–188
  * Faulstich et al. 2023 Faulstich, F. M.; Sturmfels, B.; Sverrisdóttir, S. Algebraic varieties in quantum chemistry. _arXiv:2308.05258_ 2023,
  * Borovik et al. 2023 Borovik, V.; Sturmfels, B.; Sverrisdóttir, S. Coupled cluster degree of the grassmannian. _arXiv:2310.15474_ 2023,
  * Schneider 2009 Schneider, R. Analysis of the Projected Coupled Cluster Method in Electronic Structure Calculation. _Numerische Mathematik_ 2009, _113_ , 433–471
  * Rohwedder 2013 Rohwedder, T. The Continuous Coupled Cluster Formulation for the Electronic Schrödinger Equation. _ESAIM: Mathematical Modelling and Numerical Analysis_ 2013, _47_ , 421–447
  * Rohwedder and Schneider 2013 Rohwedder, T.; Schneider, R. Error Estimates for the Coupled Cluster Method. _ESAIM: Mathematical Modelling and Numerical Analysis_ 2013, _47_ , 1553–1582
  * Laestadius and Kvaal 2018 Laestadius, A.; Kvaal, S. Analysis of the extended coupled-cluster method in quantum chemistry. _SIAM Journal on Numerical Analysis_ 2018, _56_ , 660–683
  * Laestadius and Faulstich 2019 Laestadius, A.; Faulstich, F. M. The coupled-cluster formalism–a mathematical perspective. _Molecular Physics_ 2019, _117_ , 2362–2373
  * Faulstich et al. 2019 Faulstich, F. M.; Laestadius, A.; Legeza, Ö.; Schneider, R.; Kvaal, S. Analysis of the tailored coupled-cluster method in quantum chemistry. _SIAM Journal on Numerical Analysis_ 2019, _57_ , 2579–2607
  * Hassan et al. 2023 Hassan, M.; Maday, Y.; Wang, Y. Analysis of the single reference coupled cluster method for electronic structure calculations: the full-coupled cluster equations. _Numerische Mathematik_ 2023, _155_ , 121–173
  * Hassan et al. 2023 Hassan, M.; Maday, Y.; Wang, Y. Analysis of the Single Reference Coupled Cluster Method for Electronic Structure Calculations: The Discrete Coupled Cluster Equations. _arXiv:2311.00637_ 2023,
  * Crawford and Schaefer III 2007 Crawford, T. D.; Schaefer III, H. F. An introduction to coupled cluster theory for computational chemists. _Reviews in computational chemistry_ 2007, _14_ , 33–136
  * 71 Grayson, D. R.; Stillman, M. E. Macaulay2, a software system for research in algebraic geometry. Available at http://www2.macaulay2.com
  * Breiding and Timme 2018 Breiding, P.; Timme, S. HomotopyContinuation. jl: A package for homotopy continuation in Julia. Mathematical Software–ICMS 2018: 6th International Conference, South Bend, IN, USA, July 24-27, 2018, Proceedings 6. 2018; pp 458–465
  * Breiding et al. 2024 Breiding, P.; Kohn, K.; Sturmfels, B. _Metric algebraic geometry_ ; Springer Nature, 2024
  * Bates et al. 2023 Bates, D. J.; Breiding, P.; Chen, T.; Hauenstein, J. D.; Leykin, A.; Sottile, F. Numerical Nonlinear Algebra. _arXiv:2302.08585_ 2023,
  * Garcia and Zangwill 1979 Garcia, C. B.; Zangwill, W. I. Finding all solutions to polynomial systems and other systems of equations. _Mathematical Programming_ 1979, _16_ , 159–176
  * Morgan 2009 Morgan, A. _Solving polynomial systems using continuation for engineering and scientific problems_ ; SIAM, 2009
  * Sommese et al. 2005 Sommese, A. J.; Wampler, C. W., et al. _The Numerical solution of systems of polynomials arising in engineering and science_ ; World Scientific, 2005
  * Davidenko 1953 Davidenko, D. On a new method of numerical solution of systems of nonlinear equations. Proceedings of the USSR Academy of Sciences. 1953; pp 601–602
  * Davidenko 1953 Davidenko, D. On the approximate solution of systems of nonlinear equations. _Ukrainian Mathematical Journal_ 1953, _5_ , 196–206
  * Morgan et al. 1992 Morgan, A. P.; Sommese, A. J.; Wampler, C. W. Computing singular solutions to polynomial systems. _Advances in Applied Mathematics_ 1992, _13_ , 305–327
  * Verschelde 1999 Verschelde, J. Algorithm 795: PHCpack: A general-purpose solver for polynomial systems by homotopy continuation. _ACM Transactions on Mathematical Software (TOMS)_ 1999, _25_ , 251–276
  * Bates et al. 2006 Bates, D. J.; Hauenstein, J. D.; Sommese, A. J.; Wampler, C. W. Bertini: Software for numerical algebraic geometry. 2006
  * Chen et al. 2014 Chen, T.; Lee, T.-L.; Li, T.-Y. Hom4PS-3: a parallel numerical solver for systems of polynomial equations based on polyhedral homotopy continuation methods. Mathematical Software–ICMS 2014: 4th International Congress, Seoul, South Korea, August 5-9, 2014. Proceedings 4. 2014; pp 183–190
  * Lee et al. 2008 Lee, T.-L.; Li, T.-Y.; Tsai, C.-H. HOM4PS-2.0: a software package for solving polynomial systems by the polyhedral homotopy continuation method. _Computing_ 2008, _83_ , 109–133
  * Van Voorhis and Head-Gordon 2000 Van Voorhis, T.; Head-Gordon, M. Benchmark variational coupled cluster doubles results. _J. Chem. Phys._ 2000, _113_ , 8873–8879
  * Sun et al. 2020 Sun, Q.; Zhang, X.; Banerjee, S.; Bao, P.; Barbry, M.; Blunt, N. S.; Bogdanov, N. A.; Booth, G. H.; Chen, J.; Cui, Z.-H., et al. Recent developments in the PySCF program package. _The Journal of chemical physics_ 2020, _153_
  * Sun et al. 2018 Sun, Q.; Berkelbach, T. C.; Blunt, N. S.; Booth, G. H.; Guo, S.; Li, Z.; Liu, J.; McClain, J. D.; Sayfutyarova, E. R.; Sharma, S., et al. PySCF: the Python-based simulations of chemistry framework. _Wiley Interdisciplinary Reviews: Computational Molecular Science_ 2018, _8_ , e1340
  * Sun 2015 Sun, Q. Libcint: An efficient general integral library for g aussian basis functions. _Journal of computational chemistry_ 2015, _36_ , 1664–1671
  * Bulik et al. 2015 Bulik, I. W.; Henderson, T. M.; Scuseria, G. E. Can single-reference coupled cluster theory describe static correlation? _Journal of chemical theory and computation_ 2015, _11_ , 3171–3179

##  7 Dissociating systems of hydrogen

We begin by reporting the full energy spectrum of the Hamiltonian describing
(H2)2 in D2h symmetry together with all CCD roots, see Figure 13(a). Given the
large number of solutions, we minimize over the difference between spectral
values of the Hamiltonian and the found CCD energies. This provides further
insight into energies that are well approximated by CCD, see Figure 13(b).

(a)

(b)

Figure 13: (a) Full spectrum (b) Minimal difference.

Similarly, we proceed for (H2)2 in D∞h symmetry. We reporting the full energy
spectrum of the Hamiltonian describing (H2)2 in D∞h symmetry together with all
CCD roots, see Figure 14(a) Given the large number of solutions, we minimize
over the difference between spectral values of the Hamiltonian and the found
CCD energies. This provides further insight into energies that are well
approximated by CCD, see Figure 14(b).

(a)

(b)

Figure 14: (a) Full spectrum (b) Minimal difference.

##  8 H4 disturbed on a circle

Similarly, we proceed for H4 disturbed on a circle. We report the full energy
spectrum of the Hamiltonian describing H4 disturbed on a circle together with
all CCD roots, see Figure 15(a). Given the large number of solutions, we
minimize over the difference between spectral values of the Hamiltonian and
the found CCD energies. This provides further insight into energies that are
well approximated by CCD, see Figure 15(b).

(a)

(b)

Figure 15: (a) Full spectrum (b) Minimal difference.

##  9 Lithium hydride

We report the full energy spectrum of the Hamiltonian describing LiH together
with all CCD roots, see Figure 16(a) Given the large number of solutions, we
minimize over the difference between spectral values of the Hamiltonian and
the found CCD energies. This provides further insight into energies that are
well approximated by CCD, see Figure 14(b).

(a)

(b)

Figure 16: (a) Full spectrum (b) Minimal difference.


# Synchronous imaging for rapid visualization of complex vibration profiles in
electromechanical microresonators

Y. Linzon<EMAIL_ADDRESS><EMAIL_ADDRESS>Department of Physics
and Optical Engineering, Ort Braude College, PO Box 78, Karmiel 21982, Israel
D. J. Joe  School of Applied and Engineering Physics, Cornell University,
Ithaca, New York 14853, USA  S. Krylov  School of Mechanical Engineering,
Faculty of Engineering, Tel Aviv University, Tel Aviv 69978, Israel  B. Ilic
School of Applied and Engineering Physics, Cornell University, Ithaca, New
York 14853, USA  J. Topolancik  School of Applied and Engineering Physics,
Cornell University, Ithaca, New York 14853, USA  J. M. Parpia  School of
Applied and Engineering Physics, Cornell University, Ithaca, New York 14853,
USA  H. G. Craighead  School of Applied and Engineering Physics, Cornell
University, Ithaca, New York 14853, USA

###### Abstract

Synchronous imaging is used in dynamic space-domain vibration profile studies
of capacitively driven, thin n+ doped poly-silicon microbridges oscillating at
rf frequencies. Fast and high-resolution actuation profile measurements of
micromachined resonators are useful when significant device nonlinearities are
present. For example, bridges under compressive stress near the critical Euler
value often reveal complex dynamics stemming from a state close to the onset
of buckling. This leads to enhanced sensitivity of the vibration modes to
external conditions, such as pressure, temperatures, and chemical composition,
the global behavior of which is conveniently evaluated using synchronous
imaging combined with spectral measurements. We performed an experimental
study of the effects of high drive amplitude and ambient pressure on the
resonant vibration profiles in electrically-driven microbridges near critical
buckling. Numerical analysis of electrostatically driven post-buckled
microbridges supports the richness of complex vibration dynamics that are
possible in such micro-electromechanical devices.

###### pacs:

87.64.-t, 07.10.Cm, 85.85.+j, 05.45.-a

##  I Introduction

Suspended resonant nano- and micro-electro-mechanical systems (NEMS and MEMS)
find use in versatile applications, such as ultra-sensitive mass detectors, rf
filters, and switching devices Craig (2000); Roukes1 (2005). As device
miniaturization advances, optimization of the overall characteristics in high-
frequency MEMS/NEMS resonators becomes increasingly complex and linked with
various mechanical, electrical, thermal and optical parameters of the system
and its environment. This compounds their seemingly superior sensitivity to
environmental conditions, such as the pressure, temperature and chemical
composition of the surrounding gas.

In the characterization of NEMS and MEMS under periodic electrical actuation,
vibration profile (VP) measurements are important in conjunction with
frequency-domain spectral studies Craig (2000); Roukes1 (2005); Pressure_dep
(2011); Carr (1999); Max (2000). While the latter yield important mechanical
properties, the former can be useful in many applications, including
optimization of the excitation parameters, aiding the identification of sites
most effective for localized functionalization to enable sensing, and in
studies of dissipation effects such as intrinsic and pressure-dependent
damping Pressure_dep (2011). Space-domain profiling is crucial in the presence
of significant nonlinearities where boundary conditions become critical RonLif
(2008). For instance, in typical capacitive electrical drive configuration,
the force between the grounded substrate and a device fabricated by patterned
suspended poly-crystalline silicon (polySi) film (serving as an electrode) is
inherently nonlinear with the drive amplitude and film stress Craig (2000);
Roukes1 (2005); Pressure_dep (2011); Carr (1999); Max (2000); RonLif (2008).

With interferometric reflection-mode optical transduction Carr (1999),
thermoacoustic effects can significantly modify the effective device stiffness
or induce autoparametric optical drive Max (2000). We observe all these
effects to be significant in NEMS/MEMS devices defined on films with low
compressive residual stress under applied loads near to the Euler critical
value Buckled_exp (1999). At the critical point, MEMS devices are most
sensitive to changes induced by stress variations in chemically-reactive
coatings SensorReview (2011); Darren (2010). Of additional practical interest
are possible non-uniformities of mechanical and electrical film properties
across the wafer, originating from growth processes and application of
anisotropic etch, which directly affect each circumferentially-clamped
microresonator ($\mu$R) Buckled_exp (1999). Fast space-domain visualization of
resonant VPs serves as a direct means to study the physics of all these
effects on the single device micro-scale during its actuation.

VPs in MEMS are traditionally imaged optically with vibrometric Pressure_dep
(2011), interferometric intMZ (2001); intHeterodyne (2008), or stroboscopic
Strobo (2002); YL (2010) microscopy. Recently, spatiotemporal evaluation of
resonant VPs in high-frequency MEMS $\mu$Rs was demonstrated using resonant
realtime synchronous imaging (RSI) with a pulsed low duty-cycle nanosecond
laser YL (2010). The main feature of stroboscopic imaging is a rapid
production of time-resolved interference pattern movies and static profiles,
as well as the fast evaluation of VPs, thus supplanting scanned motorized
probes that are expensive and inherently slow. This technique is applicable
with mechanical resonant frequencies up to $f_{0}\simeq 1$GHz. In this paper,
we use RSI YL (2010) in averaging mode to rapidly characterize the VPs in
bridge $\mu$Rs close to critical stress as a function of the drive amplitude
and ambient pressure. The effects of high drive nonlinearity and air damping
on the resonator VPs are directly monitored.

##  II Experimental method

$\mu$Rs are fabricated by standard top-down micromachining methods, where
bridges are defined photolithographically on compressively-stressed n+ doped
polySi films, deposited by low pressure chemical vapor deposition over a
sacrificial oxide layer and wet-etch released. Upon release of doubly clamped
bridges, residual stress is relieved through buckling Buckled_exp (1999);
Darren (2010). The devices are driven capacitively with the moving $\mu$R
serving as an electrode and the silicon substrate serving as a bottom ground
electrode. The inset in Fig. 1(a) shows the schematics of a buckled $\mu$R
cross section, as well as definitions we use, and Fig. 1(b) shows SEM images
of our released bridges.

Figure 1: (a) Schematic cross section of the devices studied and definitions
of optical quantities used in the analysis. (b) SEM images of bridges of
dimensions: 25$\times$6$\times$0.12 $\mu$m3 (left, slightly post-buckled), and
20$\times$1$\times$0.14 $\mu$m3 (right, flat), both with $\sim$220 nm static
elevations. Figure 2: (Color online) (a) Schematics of the experimental setup.
(b) Calibration curves for synchronous imaging assuming a device with film
thickness $t$=138 nm and static midpoint elevation $d_{0}$=220 nm. Left:
absolute reflection coefficient $R$. Right: differential reflection $\Delta
R/R_{0}$. With negative values of $\Delta R$, the intensity contrast in the
image is negative.

In Fig. 2(a) a schematic of our RSI configuration is illustrated. A dual
channel pulse source feeds the $\mu$R and optical imaging pulse source
($\lambda_{0}$=661.5 nm) in synchrony. The collimated illumination at a
glancing angle $\theta\simeq 40^{\circ}$ is reflected off the $\mu$R and
collected by an objective lens followed by a $4f$ lens pair, the latter of
which is used for spatial waveform filtering at the Fourier plane with a phase
mask Fourier (1978). The outgoing light is finally imaged on a standard CCD
camera. Changes in the reflection with respect to the static image of the
$\mu$R, due to resonant motion, are monitored as a function of the rf source
frequency $f_{0}$, voltage and phase. The pressure within the chamber is set
with a vacuum pump and venting tubes, and monitored via a Pirani gauge. In
order to calibrate the physical VPs from measured reflection images, an
interferometric analysis is carried out in the out-of-plane direction (shown
in Fig. 1(a)), as detailed below. Application of a 50% duty-cycle to the
imaging pulses (full synchronization with the capacitive drive), high in-phase
sensitivity to _average_ differential actuation amplitudes is attained at the
expense of lost temporal resolution.

For calibration of the physical VPs from measured reflectivity images, a
Fabry-Pèrot interferometer multilayer analysis, as a function of the total
elevation Interference_book (1995), is performed using knowledge of the static
film elevation profile $d_{0}$, thickness $t$, and the refractive indices of
the film (n-doped polySi, $n$=3.916) and substrate $n_{S}$ (single crystal Si,
$n$=3.834). The reflectance coefficients are calculated from the effective
reflectivity matrix, assuming nearly normal incidence:

$\displaystyle M_{total}=M_{2}\cdot
M_{1}=\left(\begin{array}[]{cc}\cos\delta_{2}&\frac{i}{n_{2}}\sin\delta_{2}\\\
in_{2}\sin\delta_{2}&\cos\delta_{2}\\\ \end{array}\right)\cdot$ (3)
$\displaystyle\left(\begin{array}[]{cc}\cos\delta_{1}&\frac{i}{n_{1}}\sin\delta_{1}\\\
in_{1}\sin\delta_{1}&\cos\delta_{1}\\\ \end{array}\right)$ (6)

Figure 3: (Color online) Drive amplitude dependence of VPs in the fundamental
mode of a critically upward buckled resonator. (a) Reference image of the
static bridge. (b) Frequency domain spectra with low ac actuation voltage and
a constant 5 V dc bias. (c) Static height profile of the bridge along the $Y$
direction taken from AFM measurements (0 is defined as the height of the
trench and known film thickness of 140 nm is subtracted on the bridge).
(d),(e) Measured synchronous images at $f_{0}$=4.1 MHz with different ac
amplitudes; (f)-(i) corresponding VPs integrated along $Y$ [in (f),(g),
$X$-profiles], and along $X$ [in (h),(i), $Y$-profiles]. Vertical arrows in
(f),(g) indicate diminished actuation signals with high excitation.

where $\delta_{j}=k_{j}d_{j}$ is the effective phase of layer $j$ and $k_{j}$
is the wave number. Denoting $R_{0}(x,y)$ the reference image of the static
reflection, the contrast signal measured during actuation corresponds to:

$R_{meas}(x,y)=\frac{R(x,y)-R_{0}(x,y)}{R_{0}(x,y)}\equiv\frac{\Delta
R(x,y)}{R_{0}}$ (7)

Under full synchronization of the sampling beam with the drive frequency and
phase, the observed average amplitude $<A>$ at transverse position $(x,y)$ is
then:

$\langle A\rangle(x,y)=d_{1}(x,y)-d_{0}(x,y)$ (8)

With full sampling synchronization and a periodic unipolar square wave
excitation, the observed average amplitude $\langle A\rangle$ at transverse
position $(x,y)$ is $\langle A\rangle=A_{max}/2$. The glancing angle of
illumination (see Fig. 1(a)) introduces an additional scaling factor
$\cos\theta$. The peak VP average amplitude, as measured away from the static
beam position, is then:

$A_{max}(x,y)=\frac{2}{\cos\theta}[d_{1}(x,y)-d_{0}(x,y)]$ (9)

giving rise to a normalization factor of 2.83 in our implementation, with
$\theta=45^{\circ}$. A glancing angle also introduces shadow effects at the
resonator edges, which we can easily eliminate with appropriate phase masks at
the Fourier plane (see Fig. 1(a)). We consider only reflection variations at
the positions of the $\mu$R itself to constitute its real VPs. An example of
the total and differential reflectance curves as function of the total
elevations is shown in Fig. 2(b). Using Eqs. (1)-(4), together with the
calibration curve, the average maximum amplitude profiles are estimated. Here
we will concentrate on characterizations of the fundamental (lowest) resonant
mode.

Figure 4: (Color online) Pressure-dependent studies of VPs in the fundamental
mode of a slightly post-buckled resonator. (a) Reference image of the static
bridge. (b) Frequency domain spectra (inset) and inverse quality factors
(dissipation) as a function of chamber pressure, under continuous 315 mV ac
drive and a 5 V dc bias. (c) Measured interferometric images at $f_{0}$=2.4
MHz as a function of the pressure, and (d) overlaid $Y$-integrated
$X$-profiles of vibration.

##  III Results and Discussion

In Fig. 3 we study a $\mu$R of dimensions (25$\times$6$\times$0.14 $\mu$m3)
and a midpoint elevation of 220 nm under low pressure settings (P$<$1 Torr).
The undriven $\mu$R is almost flat (see Fig. 3(c)) whereas other slightly
longer devices exhibit noticeable static upward buckling, suggesting the
existence of a compressive force whose magnitude is close to critical load.
Highly buckled resonators have been found as hard to drive electrostatically.
Figure 3(a) shows a static optical image of the unactuated device in its
initial reference configuration. In Fig. 3(b), the frequency response under
low-voltage actuation is shown. Even with drive amplitudes as low as 45 mV and
a dc bias of 5 V, we observe the formation of Duffing nonlinearity RonLif
(2008) and significant spectral broadening, with a sudden frequency detuning
between 35 and 45 mV drive voltage. An AFM measurement of the static bridge
height profile, in the transverse ($Y$) direction, is shown in Fig. 3(c), and
the profile is uniform in the axial ($X$) direction, indicating a shell-like
bridge profile. RSI images with intermediate and high ac drive voltages, at a
frequency corresponding to the maximum resonant amplitude, optimal phase and a
dc bias of 5 V, are shown in Figs. 3(d) and 3(e), respectively. Following the
image analysis for amplitude calibration, as detailed in the experimental
section, and integration along the beam width ($Y$), the peak VP amplitude
$X$-profiles are shown in Figs. 3(f) and 3(g), respectively. The peak VP
amplitude $Y$-profiles, integrated and averaged along $X$, are also shown in
Figs. 3(h) and 3(i). With intermediate drive amplitudes the VP shapes are as
shown in Fig. 3(f) and 3(h). With high drive amplitudes, central regions on
the beam appear to undergo diminished displacement at the original frequency
(Fig. 2(g)). However, detuning of the imaging frequency in these cases to
values near multiples of the fundamental mechanical frequency and the same
phase settings show some tiny components of vibration at these locations. We
interpret this observation as resulting from either nonlinear
electromechanical processes inducing transfer of energy to higher harmonics at
locations of high vibration amplitudes on the $\mu$R, or from optical
nonlinearity due to the measured response crossing extreme reflection points.
In any case, positions with diminished signal, such as the one indicated by
the vertical arrows in Figs. 3(f),(g) would clearly not be beneficial to
employ in phase-locked-loop (PLL) sensing applications, using this class of
$\mu$Rs, at this wavelength. Along the $Y$-profiles, slight localization of
the motion at the central region of the bridge is also observed with high
drives.

Figure 4 shows studies using a narrow microbridge of dimensions
(25$\times$1$\times$0.12 $\mu$m3) and a midpoint elevation of 660 nm ($t$=120
nm and $d_{0}$=660 nm) under varying ambient hydrostatic pressure and constant
driving conditions of 1.2 V ac voltage and 5 V dc bias. This bridge is
slightly buckled in the upward direction, as observed in the static optical
reference image of Fig. 4(a). Figure 4(b) shows the dissipation (inverse
quality factor, $Q^{-1}$) of the fundamental resonant mode as a function of
pressure, and corresponding spectra (inset). Different pressure ranges
correspond to well-known dominant dissipation mechanismsPressure_dep (2011);
Darren2 (2009). The total quality factor $Q$ is known to approximately scale
according to Darren2 (2009):

$1/Q=1/Q_{int}+\alpha P$ (10)

Where $Q_{int}$ is the intrinsic (material) quality factor, $\alpha$ is the
coefficient of viscous damping and $P$ is the pressure. In the data
corresponding to Fig. 4(b), a linear fit yields $Q_{int}=154$ and
$\alpha=1.83\times 10^{-4}$ $[\textmd{Torr}^{-1}]$ in this $\mu$R.

In the current experiment we have succeeded in recording RSI images of
sufficient contrast only at pressures below the viscous (gas-dominated)
regime, namely, corresponding to the intrinsic and molecular regimes in Fig.
4(b). It is estimated that the most significant limiting factors are the low
spectral signal-to-noise bandwidth (S/N) at low quality factors (below
$Q\approx$ 20) combined with diminished amplitudes of motion under external
air damping. Figure 4(c) shows RSI images of the $\mu$R as a function of
increasing pressure, with a transition from intrinsic to molecular damping.
Calibrated maximum amplitude $X$-profiles, integrated across the beam width
($Y$), are shown in Fig. 4(d). Increased errors in the VP estimations result
from diminished available S/N, giving rise to less accurate numerical fits. We
consistently find that with increasing pressure, the VPs in this $\mu$R become
suppressed around the regions close to the bridge overhang (see Fig. 4(d)).
This edge suppression effect is not observed in repeated experiments under low
pressure and drive conditions (0.3 V ac voltage and 5 V dc bias), that yield
available signal-to-noise close to the detection limit, with extracted
vibration amplitudes comparable to the highest pressure case shown here and
with more pronounced motion near the overhang.

##  IV Numerical Model

The dynamics of a compressively stressed beam are described by the equation
Nayfeh (2004); Slava (2011):

$\displaystyle EI(\frac{\partial^{4}w}{\partial
x^{4}}-\frac{\partial^{4}w_{0}}{\partial
x^{4}})-[P-\frac{EA}{2L}\int^{L}_{0}((\frac{\partial
w}{dx})^{2}-(\frac{\partial w_{0}}{\partial x})^{2})dx]\times$
$\displaystyle\frac{\partial^{2}w}{\partial x^{2}}+\rho
d\frac{\partial^{2}w}{\partial
t^{2}}=-\frac{\epsilon_{0}bV^{2}}{2(g_{0}+w)^{2}}$

Figure 5: Schematics of the numerical model.

where now $E$ is the Young’s modulus of the beam material, $I=b\times
d^{3}/12$ is the moment of inertia of the beam cross section, $A=b\times d$ is
the sectional area, $L$ is the beam length, $\rho$ is the density, and $b$ and
$d$ are the thickness and width of the beam, respectively. In addition,
$g_{0}$ is distance between the ends of the flat side of the beam and the
electrode (electrostatic gap), $\epsilon_{0}$ is the vacuum permittivity and
$V(t)$ is the time-dependent actuation voltage. In accordance with the
definitions in Fig. 5, the elevation of the beam $w(x)$, as well as the
electrostatic force $f^{e}(x,t)=-\epsilon_{0}bV^{2}/2(g_{0}+w)^{2}$, which is
calculated using the simplest parallel capacitor approximation formula, are
considered positive upwards.

Equation (5) has been reduced to the system of coupled nonlinear ordinary
differential equations by means of the Galerkin decomposition with linear
undamped eigenmodes of a straight beam used as base functions. The equaations
were solved numerically using the ODE45-solver implemented in Matlab. The
details of the formulation and numerical approach used for the analysis are
found in Slava (2011) (see also Nayfeh (2004)).

##  V Numerical Results

Figures 6-8 show numerical solutions of Eq. (5). In all cases, the actuation
voltage contained both ac and dc bias components, such that
$V(t)=V_{dc}+V_{ac}\cos(\omega t)$. Zero initial conditions, corresponding to
the post-buckled configuration of the beam in rest, were used. In all cases,
Young’s modulus $E$=150 GPa and density $\rho$=2300 kg/m3 corresponding to
polySi were used.

Figure 6: (a) Numerical frequency domain resonant response. Midpoint
deflection of the beam is shown. (b) Deflection profiles (difference between
the actual and initial elevations of the beam) averaged over a single period.
In (b) the operation frequency, $\omega$=2.45 MHz, is close to the fundamental
mode resonant frequency. In both simulations, the dimensions of the beam are
25$\times$1$\times$0.12 $\mu$m3; the electrostatic gap is $g_{0}$=660 nm; the
initial elevation of the midpoint above the beam’s ends, due to buckling, is
155 nm; input voltages are $V_{dc}=5$ V and $V_{ac}=1.2$ V, and (a) $Q$=7; (b)
$Q$=700.

Figure 6(a) shows the resonant response of the beam with dimensions
(25$\times$1$\times$0.12 $\mu$m3) and electrostatic gap $g_{0}$=660 nm. The
axial force was chosen such that the midpoint elevation of the beam above its
ends was 155 nm. The driving voltages were $V_{dc}=5$ V and $V_{ac}=1.2$ V,
and the quality factor was $Q$=7. It is observed that the fundamental resonant
frequency is 2.45 MHz, which is close to the experimentally observed value
(Fig. 4(b)). The corresponding resonant displacement profile, averaged over a
single period, is shown in Fig. 6(b). Small initial imperfection of 0.05 in
the initial buckled height, corresponding to an excitation of the second anti-
symmetric mode, was introduced in order to allow non-symmetric mode shapes of
the beam. Calculations show that while the actual beam profiles are dominated
by the fundamental mode, the resonant deflection profiles (i.e., the
differences between the initial buckled shapes and the actual, time dependent
shapes of the vibrating beam) could be more complex.

Figure 7: Dynamic snapshots of the deflection profiles (differences between
actual and initial elevations of the beam) corresponding to different time
sections within a single resonant cycle. The beam dimensions are
25$\times$1$\times$0.12 $\mu$m3; the electrostatic gap is $g_{0}$=220 nm; the
initial elevation of the midpoint above the beam’s ends, due to buckling, is
98 nm; input voltages are $V_{dc}=1$ V and $V_{ac}=350$ mV, and $Q$=1000. The
operation frequency is $\omega$=1.495 MHz. Nine (symmetric and skew-symmetric)
base functions are preserved in the reduced-order model.

A decrease in the initial separation between the beam and the electrode
results in an increased contribution of higher modes in the resonant VPs.
Figure 7 shows time-resolved snapshots of the VPs (relative displacements from
equilibrium) with the same dimensions as in Fig. 6, but with $g_{0}$=220 nm;
Figure 8 shows the same vibration profile as averaged over a single period.
Small initial imperfection of 0.05 in the initial buckled height,
corresponding to a contribution of the second buckled mode, was again used as
an initial condition. Complex displacement profiles are clearly observed in
this case as well.

Figure 8: Deflection profiles (difference between the actual and initial beam
elevations) averaged over a single period, corresponding to the results in
Fig. 7 (same parameters).

##  VI Conclusion

Synchronous imaging has been demonstrated as a robust method for direct and
rapid observations of gradual changes in resonant vibration profiles of
electromechanical microresonators under varying conditions of drive and
ambient pressure. Synchronous imaging can serve as a useful tool for studying
fundamental processes in resonant MEMS/NEMS, as well as for identification of
favorable device positions most suitable for sensitive phase-locked
applications, such as sensors, filters and switches. Numerical analysis of
electrostatically driven post-buckled microbridges supports the richness of
the complex resonant vibrations that are possible in these micro-
electromechanical systems.

##  VII Acknowledgments

This research was funded by the National Science Foundation (grants
DMR-0908634 and DMR-0520404) and Analog Devices. Fabrication was performed at
the Cornell Nanoscale science and technology Facility.

## References

  * Craig (2000) H. G. Craighead, Science 290, 5496 (2000).
  * Roukes1 (2005) K. L. Ekinci and M. L. Roukes, Rev. Sci. Instrum. 76, 061101 (2005).
  * Pressure_dep (2011) R. C. Tung RC, J. W. Lee, H. Sumali, and A. Raman, J. Micromech. Microeng. 21, 025003 (2011); R. A. Bidkar, R. C. Tung, A. A. Alexeenko, H. Sumali, and A. Raman, Appl. Phys. Lett. 94, 163117 (2009); H. Sumali, J. Micromech. Microeng. 17, 2231 (2007).
  * Carr (1999) D. W. Carr, S. Evoy, L. Sekaric, H. G. Craighead, and J. M. Parpia, Appl. Phys. Lett. 75, 920 (1999).
  * Max (2000) M. Zalalutdinov, A. Zehnder, A. Olkhovets, S. Turner, L. Sekaric, B. Ilic, D. Czaplewski, J. M. Parpia, and H. G. Craighead, Appl. Phys. Lett. 79, 695 (2001); B. Ilic, S. Krylov, K. Aubin, R. Reichenbach, and H. G. Craighead, Appl. Phys. Lett. 86, 193114 (2005).
  * RonLif (2008) R. Lifshitz and M. C. Cross, _Review of Nonlinear Dynamics and Complexity_ (Wiley, Meinheim, 2008), Vol. I, pp. 1-52.
  * Buckled_exp (1999) W. Fang, C.-H. Lee, and H.-H. Hu, J. Micromech. Microeng. 9, 236 (1999).
  * SensorReview (2011) A. Boisen, S. Dohn, S. S. Keller, S. Schmid, and M. Tenje, Rep. Prog. Phys. 74, 036101 (2011).
  * Darren (2010) D. R. Southworth, L. M. Bellan, Y. Linzon, H. G. Craighead, and J. M. Parpia, Appl. Phys. Lett. 96, 163503 (2010).
  * intMZ (2001) G. G. Fattinger and P. T. Tikka, Appl. Phys. Lett. 79, 290 (2001).
  * intHeterodyne (2008) K. Kokkonen and M. Kaivola, Appl. Phys. Lett. 92, 063502 (2008).
  * Strobo (2002) C. Rembe and R. S. Muller, J. Microelectromech S. 11, 479 (2002).
  * YL (2010) Y. Linzon, S. Krylov, B. Ilic, D. R. Southworth, R. A. Barton, B. R. Cipriany, J. D. Cross, J. M. Parpia, and H. G. Craighead, Opt. Lett. 15, 2654 (2010).
  * Darren2 (2009) D. R. Southworth, H. G. Craighead, and J. M. Parpia, Appl. Phys. Lett. 94, 213506 (2009).
  * Fourier (1978) J. D. Gaskill, _Linear Systems, Fourier Transforms, and Optics_ (Wiley, New York, 1978).
  * Interference_book (1995) M. Bass, _Handbook of Optics_ , 2nd ed. (McGraw-Hill, San Francisco, 1995), Vol. I, pp. 42.10-42.14.
  * Slava (2011) S. Krylov, B. R. Ilic, and S. Lulinsky, Nonlinear Dyn. 66, 403 (2011).
  * Nayfeh (2004) S. A. Emam and A. H. Nayfeh, Nonlinear Dyn. 35, 1 (2004).