# Homotopy continuation methods for coupled-cluster theory in quantum
chemistry

Fabian M. Faulstich<EMAIL_ADDRESS>Department of Mathematics,
University of California, Berkeley  Andre Laestadius
<EMAIL_ADDRESS>Department of Computer Science, Oslo Metropolitan
University  Hylleraas Centre for Quantum Molecular Systems, Department of
Chemistry, University of Oslo

###### Abstract

Homotopy methods have proven to be a powerful tool for understanding the
multitude of solutions provided by the coupled-cluster polynomial equations.
This endeavor has been pioneered by quantum chemists that have undertaken both
elaborate numerical as well as mathematical investigations. Recently, from the
perspective of applied mathematics, new interest in these approaches has
emerged using both topological degree theory and algebraically oriented tools.
This article provides an overview of describing the latter development.

##  I Introduction

Coupled-cluster (CC) theory is a widely acclaimed, high-precision wavefunction
approach that is used in quantum chemistry and is of great interest to both
practitioners as well as theoreticians Bartlett and Musiał (2007). The origin
of CC theory dates back to 1958 when Coester proposed to use an exponential
parametrization of the wave function Coester (1958). This parametrization was
independently derived by Hubbard Hubbard (1957) and Hugenholtz Hugenholtz
(1957) in 1957 as an alternative to summing many-body perturbation theory
(MBPT) contributions order by order. A milestone of CC theory is the work by
Čížek from 1966 Čížek (1966). In this work, Čížek discussed the foundational
concepts of second quantization (as applied to many-fermion systems), normal
ordering, contractions, Wick’s theorem, normal-ordered Hamiltonians (which was
a novelty at that time), Goldstone-style diagrammatic techniques, and the
origin of the exponential wave function ansatz. He moreover derived the
connected cluster form of the Schrödinger equation and proposed a general
recipe for how to produce the energy and amplitude equations through
projections of the connected cluster form of the Schrödinger equation on the
reference and excited determinants, which was illustrated using the CC doubles
(CCD) approximation. This work also reports the very first CC computations,
using full and linearized forms of CCD, for nitrogen (treated fully at the _ab
initio_ level) and benzene (treated with a PPP model Hamiltonian). For a more
detailed history of CC theory, several reviews have been written by, e.g.,
Kümmel Kümmel (1991) and Čížek Čížek (1991). Other articles that provide
insight into the history and development of CC theory include those by
Bartlett Bartlett (2005), Paldus Paldus (2005), Arponen Arponen (1991), and
Bishop Bishop (1991).

The CC equations are a set of nonlinear algebraic equations that are typically
solved using (quasi) Newton-type methods Press  _et al._ (2007). Each equation
in the set corresponds to a projection on a specific excitation of the
reference determinant, and the number of equations increases rapidly with the
size of the system and allowed excitations under consideration. Since the CC
equations are a set of nonlinear algebraic equations, there are multiple roots
to the CC equations. The existence of multiple roots to the CC equations can
present challenges for the practical implementation of the method. For
example, the roots can be difficult to adequately converge to, and the
convergence properties of the iterative solution methods can strongly depend
on the employed initial guess. Furthermore, excited states can also be
targeted with equation-of-motion (EOM)-CC Rowe (1968); Monkhorst (1977); Koch
_et al._ (1990); Stanton and Bartlett (1993); Koch  _et al._ (1994), where an
initial “ground-state” CC calculation is the starting point of a non-Hermitian
diagonalization problem.

Although great progress has been made along the fundamental and mathematical
study of the root structure of the CC equations, which include homotopy
methods applied to the CC methodology, progress and wide-spread applications
have been hampered by the high dimensionality and non-linearity of the CC
equations, as well as the steep scaling of algebro-computational methods. The
first study on this topic dates back to 1978, where Živković and Monkhorst
investigated the singularities and multiple solutions of the single-reference
CC equations revealing conditions for reality and the maximum multiplicity of
solutions Živković and Monkhorst (1978). In 1998, Kowalski and Jankowski
revived homotopy methods in connection with CC theory and used them to solve a
CCD system Kowalski and Jankowski (1998a). This was followed by a fruitful
collaboration with Piecuch, who worked on multiple solutions to the single
reference CC and state-universal MRCC equations Paldus  _et al._ (1993);
Piecuch  _et al._ (1990). Piecuch and Kowalski extended the application of
homotopy methods to CC single and doubles-(CCSD)-, CC singles, doubles and
triples-(CCSDT)-, and CC single, doubles, triples and
quadruples-(CCSDTQ)-equations Piecuch and Kowalski (2000) for a 4-electron
system described by a minimum basis set. They also introduced the formalism of
$\beta$-nested equations and proved the Fundamental Theorem of the $\beta$-NE
Formalism, which enabled them to explain the behavior of the curves connecting
multiple solutions of the various CC polynomial systems, i.e., from CCSD to
CCSDT, CCSDT to CCSDTQ, etc. In Kowalski and Piecuch (2000a), Piecuch and
Kowalski used homotopy methods to determine all solutions of nonlinear state-
universal multireference CCSD equations based on the Jeziorski-Monkhorst
ansatz, proving two theorems that provided an explanation for the observed
intruder solution problem. In a sequel work Kowalski and Piecuch (2000b), they
used homotopy methods to obtain all solutions of the generalized Bloch
equation, which is nonlinear even in a CI parametrization. Further articles
utilizing homotopy methods arose in the late ’90s where, amongst other things,
symmetry breaking processes in non-linear CC formulation are explained
Kowalski and Jankowski (1998b); Jankowski and Kowalski (1999a, b, c, d).

In addition to examining the latest developments in homotopy techniques in CC
theory, this article aims to present new findings that result from the
application of homotopy approaches using both topological degree theory and
algebraic geometry tools, from an applied mathematics perspective. These
approaches are essential for expanding the scope of mathematical
investigations beyond the ground state. While local analyses using strong
monotonicity Schneider (2009); Rohwedder (2013); Rohwedder and Schneider
(2013); Laestadius and Kvaal (2018); Faulstich  _et al._ (2018); Laestadius
and Faulstich (2019); Faulstich  _et al._ (2019); Kvaal  _et al._ (2020);
Faulstich  _et al._ (2023) have been useful since Schneider’s seminal work
Schneider (2009), they only provide a limited perspective compared to the
broader study of the CC equations using topological degree theory or algebraic
geometry Csirik and Laestadius (2023); Faulstich and Oster (2022). By adopting
these tools, we can gain a more comprehensive understanding of the
mathematical structure of the CC equations and their solutions, providing new
avenues for further explorations and refinements of the theory.

This article is organized as follows: In Section II we briefly introduce the
CC theory including concepts like the CC variety as well as the equivalence
between a wavefunction being an eigenstate and the CC parametrization
(standard- and EOM-CC). In Section III we recall the concept of homotopy
continuation techniques, highlighting important edge-cases to consider. In
Section IV, we elaborate on the root structure of polynomial systems and the
challenges involved with using (quasi) Newton-type methods. In Section V we
then present upper bounds on the number of roots of the CCSD equations based
on Bézout’s theorem and outline how to establish improved bounds using the
Bernstein-Khovanskii-Kushnirenko theorem. In Section VI, we discuss the
mathematical existence of solution curves connecting CC roots at different
truncation levels. We also present an energy error estimate valid for any
approximate eigenstate of exponential (CC) form. Both these results have
appeared previously in the mathematical literature Csirik and Laestadius
(2023). We conclude in Section VII by summing up the main points as well as
presenting a brief outlook on future work.

##  II Coupled-cluster theory

The underlying idea of coupled-cluster (CC) theory is the exponential
parametrization of the targeted wavefunction $|{\Psi}\rangle$, i.e, for a
given reference determinant $|{\Phi_{0}}\rangle$, we make the ansatz

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

where $\hat{T}$ is the new unknown, called the cluster operator, which we will
define shortly. Using this ansatz, we find

$\mathcal{H}|{\Psi}\rangle=E|{\Psi}\rangle~{}\Leftrightarrow~{}e^{-\hat{T}}\mathcal{H}e^{\hat{T}}|{\Phi_{0}}\rangle=E|{\Phi_{0}}\rangle,$
(2)

where $\mathcal{H}$ is the considered Hamiltonian. Projecting Eq. 2 yields

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

of which the latter is used to compute the cluster operator. The cluster
operator is a linear combination of elementary particle-hole excitation
operators (vide infra), and we shall henceforth highlight the dependence of
$\hat{T}$ to its expansion coefficients by writing $\hat{T}(\bf t)$.
Elementary particle-hole excitation operators 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, highly convenient to label
the elementary particle-hole excitation operators using multi-indices $\mu$
that clarify the projections and wedge products involved (see, e.g., Schneider
(2009)). Employing this notation yields the following expression for the
cluster operator

$\hat{T}({\bf t})=\sum_{\mu}t_{\mu}\hat{X}_{\mu}.$ (4)

As is customary in the quantum chemistry community, we will use upper case
letters to describe cluster operators, and lower case letters to describe
cluster amplitudes, e.g., $\hat{T}({\bf t})$, $\hat{C}({\bf c})$,
$\hat{R}({\bf r})$ or $\hat{S}({\bf s})$. By construction, 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 amplitudes. Hence, using the
standard Galerkin projection approach, we can express the second equation and
the orthogonality constraint in Eq. (3) as

$0=\langle\Phi_{\mu}\mid\hat{\mathcal{H}}({\bf
t})\mid\Phi_{0}\rangle,\quad\forall\mu$ (5)

where $|{\Phi_{\mu}}\rangle=\hat{X}_{\mu}|{\Phi_{0}}\rangle$ and
$\hat{\mathcal{H}}({\bf t})=e^{-\hat{T}({\bf
t})}\hat{\mathcal{H}}e^{\hat{T}({\bf t})}$. Since the cluster operator is
defined by the cluster amplitudes $\bf t$, we can define the CC energy as a
function of $\bf t$, i.e.,

$\mathcal{E}_{\mathrm{CC}}({\bf t})=\langle\Phi_{0}\mid\hat{\mathcal{H}}({\bf
t})\mid\Phi_{0}\rangle.$ (6)

Note that the similarity transformed Hamiltonian $\mathcal{H}({\bf t})$ for
all amplitudes ${\bf t}$ has the same eigenvalues as the original Hamiltonian
$\mathcal{H}$. However, the CC equations in (5) do not arise from
diagonalizing the similarity-transformed Hamiltonian.

We emphasize that the orthogonality conditions in Eq. (5) yields a square
system of polynomial equations, see e.g. Helgaker  _et al._ (2014); Bartlett
and Musiał (2007); Faulstich and Oster (2022). Hence, a key object in CC
theory is the CC variety:

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

where $\mathbb{F}$ is the considered number field (either $\mathbb{R}$ or
$\mathbb{C}$), $K$ is the “system size” (given by the number of correlated
electrons, the size of one-particle basis functions as well as further
selection rules) and $\mathbb{V}$ is the cluster amplitude space, i.e., ${\bf
t}\in\mathbb{V}$ if and only if $\hat{T}({\bf
t})|{\Phi_{0}}\rangle\in\mathcal{V}$. Note that in this work, we always assume
a finite set of one-particle basis functions (i.e., orbitals), in particular,
$\mathbb{V}=\mathbb{F}^{K}$. Note that $\mathbb{F}$ can be $\mathbb{R}$ or
$\mathbb{C}$ determining if we are seeking real or complex valued amplitudes.
We emphasize that although the CC polynomial coefficients are real, the roots
to the polynomial system may not be. Mathematically, the more general theory
describing complex valued solutions is simpler than the theory describing real
valued solutions. Therefore, it is easier to consider $\mathbb{F}=\mathbb{C}$
for mathematical considerations, however, solving the truncated CC equations
for $\mathbb{F}=\mathbb{C}$ may yield complex valued energies, e.g. see
Faulstich and Oster (2022). Moreover, the majority of quantum chemistry
implementations seek real valued solutions, i.e., the CC equations are solved
for $\mathbb{F}=\mathbb{R}$.

All possible determinants that can be generated given the basis set and number
of (correlated) electrons can be represented in an excitation graph,
$G^{\mathrm{full}}$, where the vertices are the determinants
$|{\Phi_{\mu}}\rangle$ and the edges are the operators $\hat{X}_{\mu}$ Csirik,
Mihály A. and Laestadius, Andre (2023). However, we might not always want to
consider all possible determinants that can be generated but rather a subset
$G\subset G^{\mathrm{full}}$. For example, we might just consider excitations
up to a certain order or “rank”, such as CCSD ($\hat{X}_{\mu}$ contains
excitation of at most two electrons), CCSDT ($\hat{X}_{\mu}$ contains
excitations of at most three electrons), etc. Thus, we will sometimes write
$\mathcal{S}=\mathcal{S}(G)$, $\mathbb{V}=\mathbb{V}(G)$, etc. to highlight
that we consider a truncated CC scheme as dictated by $G\subset
G^{\mathrm{full}}$.

In the untruncated case, $G^{\mathrm{full}}$, we note that for given ${\bf
t}_{*}\in\mathcal{S}$ the wavefunction $\Psi=(c_{0}\hat{I}+\hat{C}({\bf
c}))|{\Phi_{0}}\rangle$ is an eigenfunction of $\mathcal{H}$, i.e., there
exists a constant $\mathcal{E}$ such that
$\hat{\mathcal{H}}|{\Psi}\rangle=\mathcal{E}|{\Psi}\rangle$, if and only if
$e^{-\hat{T}({\bf t}_{*})}(c_{0}\hat{I}+\hat{C}({\bf
c}))=r_{0}\hat{I}+\hat{R}({\bf r})$, where

$\left.\begin{aligned} \mathcal{E}_{\mathrm{CC}}({\bf
t_{*}})r_{0}+\langle\Phi_{0}\mid\hat{\mathcal{H}}({\bf t}_{*})\hat{R}({\bf
r})\mid\Phi_{0}\rangle&=\mathcal{E}r_{0}\\\
\Pi_{\mathcal{V}}\hat{\mathcal{H}}({\bf t}_{*})\hat{R}({\bf
r})|{\Phi_{0}}\rangle&=\mathcal{E}\hat{R}({\bf
r})|{\Phi_{0}}\rangle\end{aligned}\right\\}$ (8)

with $\Pi_{\mathcal{V}}$ being the orthogonal projection onto $\mathcal{V}$
and ${c_{0}=r_{0}}$. For a proof of this statement we refer to Lemma 4.1 in
Csirik and Laestadius (2023).

To shed some light on the above description of eigenstates using the CC
framework, we will give a few examples (all familiar to the quantum-chemistry
setting):

  1. (i)

We first note that $\mathcal{E}=\mathcal{E}_{\mathrm{CC}}({\bf t}_{*})$ is
equivalent to ${\bf r}_{0}\neq 0$, $\hat{R}({\bf r})=0$, and then
$|{\Psi}\rangle=e^{\hat{T}({\bf t}_{*})}|{\Phi_{0}}\rangle$.

  2. (ii)

If $\mathcal{E}=\mathcal{E}_{\mathrm{CC}}({\bf t}_{*})$ and there is
$\hat{R}({\bf r})\neq 0$, then $|{\Psi}\rangle=(r_{0}\hat{I}+\hat{R}({\bf
r}))e^{\hat{T}({\bf t}_{*})}|{\Phi_{0}}\rangle$ is another eigenstate, i.e.,
$\mathcal{E}$ is a degenerate energy level.

  3. (iii)

If $\mathcal{E}\neq\mathcal{E}_{\mathrm{CC}}({\bf t}_{*})$ and $\hat{R}({\bf
r})\neq 0$, then $|{\Psi}\rangle=(r_{0}\hat{I}+\hat{R}({\bf
r}))e^{\hat{T}({\bf t}_{*})}|{\Phi_{0}}\rangle$ is an “excited” eigenstate
with respect to the energy level $\mathcal{E}_{\mathrm{CC}}({\bf t}_{*})$).

If we let

$|{\Psi_{(k)}}\rangle=e^{\hat{T}({\bf t}_{(k)})}|{\Phi_{0}}\rangle$ (9)

denote the $k$th eigenstate to $\hat{\mathcal{H}}$, the EOM-CC method Rowe
(1968); Monkhorst (1977); Koch  _et al._ (1990); Stanton and Bartlett (1993);
Koch  _et al._ (1994) gives with $\hat{R}({\bf r}_{(1)})=\hat{I}$ that the
solutions (under certain assumptions) can be written

$|{\Psi_{(k)}}\rangle=\hat{R}({\bf r}_{(k)})e^{\hat{T}({\bf
t}_{(1)})}|{\Phi_{0}}\rangle.$ (10)

We will in the remaining part only consider the case (i) above.

The ground state (intermediately normalized) moreover corresponds to
$|{\Psi_{(1)}}\rangle=e^{\hat{T}({\bf t}_{(1)})}|{\Phi_{0}}\rangle$, where

${\bf t}_{(1)}=\mathop{\mathrm{missing}}{arg\,min}_{{\bf
t}\in\mathcal{S}}~{}\mathcal{E}_{\mathrm{CC}}(\bf t).$ (11)

We emphasize that in the untruncated case, the algebraic variety $\mathcal{S}$
(given in Eq. (7)) describes all solutions to the electronic Schrödinger
equation, given that the eigenstates have non-zero overlap with the reference
determinant $|{\Phi_{0}}\rangle$. Hence, the $k$th eigenstate (that has non-
zero overlap with the reference determinant $|{\Phi_{0}}\rangle$) can be
characterized using the min-max principle, i.e.,

$|{\Psi_{(k)}}\rangle=e^{\hat{T}({\bf t}_{(k)})}|{\Phi_{0}}\rangle$ (12)

where ${\bf t}_{(k)}$ corresponds to the argument that solves the min-max
problem

$\min\\{\max_{{\bf t}\in\mathcal{S}_{k}}~{}\mathcal{E}_{\mathrm{CC}}(\bf
t)~{}|~{}\mathcal{S}_{k}\subset\mathcal{S},\,|\mathcal{S}_{k}|=k\\}.$

In the case of untruncated CC (FCC), amplitudes that solve the FCC equations
describe an eigenstate of the Hamiltonian. Therefore, the cardinality of
$\mathcal{S}$ is equal to the number of eigenstates of $\mathcal{H}$ that are
intermediately normalized, hence, the number of roots is bounded by the number
of Slater determiants. However, when truncations are imposed, the cardinality
of $\mathcal{S}=\mathcal{S}(G)$ 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. To that end, we employ two homotopy continuation
perspectives: The first is to compute $\mathcal{S}$ in its entirety, where a
homotopy is used to connect $\mathcal{S}$ to solutions of a simpler system of
polynomial equations. The second is to characterize the physical solutions in
$\mathcal{S}$ corresponding to truncated CC equations, where a homotopy is
used to connect $\mathcal{S}$ to the FCI solutions (or at least some “higher”
truncation scheme, i.e., less truncated CC equations).

##  III Homotopy continuation

Homotopy continuation methods are well studied mathematically, we refer the
interested reader to 43; 44; 45; 46. Polynomial homotopy continuation is a
numerical method to compute solutions to systems of polynomial equations,

$F\left(x_{1},\ldots,x_{n}\right)=\left[\begin{array}[]{c}f_{1}\left(x_{1},\ldots,x_{n}\right)\\\
\vdots\\\ f_{m}\left(x_{1},\ldots,x_{n}\right)\end{array}\right]=0.$ (13)

Note that in our case, $F$ corresponds to the CC equations. In a general case,
we require $m\geq n$, however, the CC equations are a square system, i.e.,
$m=n$. The underlying idea is straightforward: to solve $F(\mathbf{x})=0$, we
construct an auxiliary system of polynomial equations, called
$G(\mathbf{x})=0$, with known zeroes, together with a homotopy that connects
both systems. More precisely, we define a family of systems
$H(\mathbf{x},\lambda)$ for $\lambda\in\mathbb{R}$ interpolating between $F$
and $G$, i.e., $H(\mathbf{x},0)=F(\mathbf{x})$ and
$H(\mathbf{x},1)=G(\mathbf{x})$. Considering one zero, $\mathbf{y}$, of
$G(\mathbf{x})$ and restricting to $\lambda\in[0,1]$,
$H(\mathbf{x},\lambda)=0$ defines a solution path
$\mathbf{x}(\lambda)\subset\mathbb{C}^{n}$ such that
$H(\mathbf{x}(\lambda),\lambda)=0$ for $\lambda\in[0,1]$ and
$\mathbf{x}(1)=\mathbf{y}$. The path is followed from $\lambda=1$ to
$\lambda=0$ to compute the solution $\mathbf{z}=\mathbf{x}(0)$. This is
equivalent to solving the initial value problem

$\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,\quad\mathbf{x}(1)=\mathbf{y}.$

Which is known as the Davidenko differential equation Davidenko (1953a, b). We
say that $\mathbf{x}(1)=\mathbf{y}$ gets tracked towards $\mathbf{x}(0)$. For
this to work, $\mathbf{x}(\lambda)$ must be a regular zero of
$H(\mathbf{x},\lambda)=0$ for every $\lambda\in(0,1]$. In the case of
nonregular solutions at $\lambda=0$, so-called endgames are employed which are
special numerical methods Morgan  _et al._ (1992).

Figure 1: 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, and
the dotted-dashed line has a unique limit.

When tracking the solutions described by the homotopy, we may encounter
different scenarios, see FIG. 1. As shown in 43, we can see that one path
(solid line) has no finite limit as $\lambda\to 0$, while the other three have
limits. One path (dotted-dashed line) has a unique limit, i.e., the endpoint
of one at $\lambda=0$ is the regular zero of the target system $F(x)$. Two
paths (dashed lines) have the same limit, and their common endpoint is an
isolated zero of $F(x)$ of multiplicity two.

We moreover wish to highlight that great progress has been made along the
lines of 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 highlight PHCpack Verschelde
(1999), Bertini Bates  _et al._ (2006), HOM4PS Chen  _et al._ (2014); Lee  _et
al._ (2008), NAG4M2 Bates  _et al._ (2023), and HomotopyContinuation.jl
Breiding and Timme (2018).

##  IV Root structure of polynomial systems

In this section, we will elaborate on the fundamental and practical importance
of understanding the root structure of the CC equations. We emphasize that it
is hard to describe the root structure of the CC equations in general because
it is a high-dimensional and non-linear system. Although this seems to be a
daunting task, the CC equations show a number of symmetries, which allow some
results regarding the system’s root structure Živković and Monkhorst (1978);
Kowalski and Jankowski (1998a); Piecuch and Kowalski (2000); Kowalski and
Piecuch (2000a, b); Faulstich  _et al._ (2018).

The root structure of a polynomial system is of fundamental importance as it
reveals, e.g., the multiplicity of the roots or whether the roots are real or
complex Hubbard  _et al._ (2001). Having information about the root structure
is also of practical importance when using approximate methods. Most commonly,
(quasi) Newton type methods are employed to find and approximate one root of
the CC equations. However, the convergence behavior of (quasi) Newton type
methods is merely locally well understood; its global convergence behavior can
be highly complicated Hubbard  _et al._ (2001); Schleicher (2002). This can be
illustrated through Newton fractals. Newton fractals are graphical
representations of the iterative process used to find the roots of a given
polynomial system using (quasi) Newton type methods. To create a Newton
fractal, one assigns to each initialization the root to which it converged. In
the case of one polynomial, one can color each point in the complex plane
according to the root to which it converges under the considered (quasi)
Newton type method, see e.g. FIG. 2 where we show the Newton fractal for
$p(z)=z^{3}-1$.

Figure 2: Newton fractal of $p(z)=z^{3}-1$. The white dots correspond to the
roots $x_{1}=1+i\cdot 0$ and $x_{2,3}=-1/2\pm i\sqrt{3}/2$. The different
colored regions, red, blue, and green, correspond to the basins of attraction
of the roots $x_{1}$, $x_{2}$, and $x_{3}$, respectively.

Clearly, FIG. 2 is a simplified perspective, since it shows the root structure
of merely one polynomial. However, it already illustrates the (potential)
issues that the application of (quasi) Newton type methods to the CC equations
can have. Figure 2 shows the three roots of the polynomial $p(z)=z^{3}-1$,
namely $x_{1}=1+i\cdot 0$ and $x_{2,3}=-1/2\pm i\sqrt{3}/2$, as well as the
color-coded basins of attraction that correspond to the respective roots. We
can clearly see that if the initial guess is close to one of the solutions,
(quasi) Newton type methods will stably converge to that solution. However, if
the initial guess is not close to a solution, but lies in one of the “fractal
branches”, the convergence of (quasi) Newton type methods becomes very
unstable Eriksen and Gauss (2019). In fact, the numerical stability can become
so poor that differences at the level of machine precision will change the
converged result. Since the CC equations are so high-dimensional, a bruit
force visualization is already for the smallest systems impractical. Yet,
since the highly intricate global convergence behavior of (quasi) Newton type
methods is rather common, one must be aware of this phenomenon, especially in
a scenario when perturbative initializations are not justifiable. In fact, if
the initial guess lies in a “fractal branch”, the employed (quasi) Newton type
methods may yield a root that is suboptimal, although the used CC ansatz (e.g.
CCSD) is well capable of describing the targeted state. To summarize this
section, using (quasi) Newton type methods to approximate a root of the CC
equations that corresponds to a targeted state, one must ensure to be in the
correct basin of convergence. Otherwise, it is unclear if the obtained CC
result is an approximation to the targeted state.

##  V Homotopy I

In this section, we discuss a path toward numerical methods that aim to find
all solutions to the CC equations. Due to the complicated structure of the
basins of attraction for (quasi) Newton type methods described in the previous
section, employing this approach for finding a CC solution heavily relies on
certain local assumptions of the initial guess, which mathematically reduces
to a perturbative picture. This importance of perturbation theory related to
CC theory is well-known in the computational chemistry community, and we refer
the interested reader to Bartlett and Musiał (2007) for a thorough and text-
book-like presentation. However, this perturbative argument is also well-known
to fail. In particular, we wish to highlight important failure modes presented
in Piecuch  _et al._ (1990) revealing the existence of algebraic branch points
as a perturbation is turned on. The existence of such branch points strongly
indicates that the commonly employed procedure of approximating one root via
(quasi) Newton-type methods starting from the Hartree-Fock reference cannot be
generally applied.  Therefore, in order to reliably extend the CC theory to
the non-perturbative regime, accessing all roots of the CC equations appears
to be inevitable. Here, homotopy continuation methods are one potential way to
achieve this, but the high dimensionality of the equations makes a direct
application of these types of procedures challenging Faulstich and Oster
(2022); Piecuch and Kowalski (2000).

In fact, establishing a good enough bound to the number of roots is already
challenging! This bound is important since it is used to initialize an
auxiliary system $G$, the employed homotopy continuation method will then
track all solutions of $G$. Clearly, when over-estimating the number of roots,
along the path $\lambda\to 0$, spurious roots will collapse, yet, if the
estimated number of roots is too large the produce will simply become
numerically intractable. A trivial bound can be obtained by observing that
Hadamard’s lemma, i.e., the expansion of $\mathcal{H}({\bf t})$, together with
commutator considerations (see e.g., Helgaker  _et al._ (2014); Bartlett and
Musiał (2007)) yields a polynomial system order of less or equal than four.
The corresponding Bézout bound is then

$\mathcal{N}\leq 4^{n_{\mathcal{K}}},$ (14)

where $\mathcal{K}$ describes the number of projective equations. However,
this number grossly overestimates the number of roots, since no further
structure of the CC equations is taken into account. A first dramatic
reduction can be obtained by noticing that the projective equations onto the
singly excited Slater determinants are of the order less or equal to three
(see e.g. Helgaker  _et al._ (2014); Bartlett and Musiał (2007)). Since our
main objective is to incorporate homotopy methods for CCSD, we will for the
remainder of this section only consider this truncation level. The above
consideration then yields

$\mathcal{N}\leq 3^{n_{s}}4^{n_{d}}$ (15)

where $n_{s}$ and $n_{d}$ are the number of projected single and double
equations, respectively Piecuch and Kowalski (2000). Note that Eq. (15) is
specific to the CCSD scheme, and generalizations to this bound for higher-
order CC schemes, such as CCSDT and CCSDTQ are also discussed in Piecuch and
Kowalski (2000).

Although this bound includes the further structure of the CCSD equations, it
is still overestimating the true number of roots. As has been recently shown,
this bound can be significantly improved by rewriting the CCSD equations as a
quadratic system Faulstich and Oster (2022). The key in this reformulation is
to notice that the anti-symmetry property of Slater determinants allows for
factorizing disconnected doubles. To that end, we define the variety
$A\subseteq\mathbb{F}^{n_{s}+2n_{d}}$ (where $\mathbb{F}$ is either
$\mathbb{R}$ or $\mathbb{C}$) over which we seek to optimize the CCSD
equations.

Formally, we define an index map $\iota$ that flattens the tuple $(i,j,a,b)$
and off-sets this compound index by $n_{s}+n_{d}$. Note that since $i$, $j$
are occupied indices, and $a$, $b$ are virtual indices, there are exactly
$n_{d}$ auxiliary indices that are obtained by $\iota$. We then define the
variety $\tilde{\mathcal{S}}$ as

$\tilde{\mathcal{S}}=\\{x\in\mathbb{F}^{n}|x_{k}-x_{i}^{a}x_{j}^{b}+x_{i}^{b}x_{j}^{a}=0,\,\forall
k=\iota(i,j,a,b)\\},$

where we introduced the short-hand notation ${n=n_{s}+2n_{d}}$. On this
variety, every projected equation $f_{\mu}$ takes the form

$f_{\mu}(\Pi_{\mathbb{V}}{\bf x})=\sum_{q,r}h^{\mu}_{q,r}x_{q}x_{r}$ (16)

where $h^{\mu}_{q,r}$ is a matrix (see Faulstich and Oster (2022)),
$\Pi_{\mathbb{V}}$ is the projection onto $\mathbb{V}$ and $\bf
x\in\tilde{\mathcal{S}}$. Note that the l.h.s. in Eq. (16) is a potentially
fourth-order polynomial in $n$ variables, whereas the r.h.s. in Eq. (16) is a
second order polynomial in $n+n_{d}$ variables. This order reduction yields
the improved bound of

$\mathcal{N}\leq 2^{n_{s}+2n_{d}},$ (17)

which is exponentially better than the existing Bézout bound (15).

Although this reduction is already quite significant, we expect that for
special cases further symmetry considerations can be incorporated yielding
even better bounds. However, Bézout type bounds are always worst-case
estimates. A better bound is expected when employing the Bernstein-Khovanskii-
Kushnirenko (BKK) theorem. Bézout type bounds arise from the idea that the
polynomials in the considered system are independent of each other, in which
case the number of roots corresponds to the product of the individual number
of roots; hence, they are worst-case bounds. BKK-type bounds are less
intuitive since the BKK theorem relates the root counting problem for a system
of polynomial equations with the theory of convex bodies. More precisely, the
BKK theorem shows that the generic number of isolated solutions to a system of
(Laurent) polynomial equations equals the mixed volume of the Newton polytopes
of the (Laurent) polynomials. The general challenge when establishing a BKK-
type bound is to compute the volume of convex bodies, i.e., the Newton
polytopes. However, general algorithms for computing the mixed volume are
exponential in the dimension. Since even computing the volume is known to be
#P-hard, a brute force approach seems doomed. That being said, the CC
equations are very structured, and a surrogate system of Newton polytopes can
be established Faulstich and Oster (2022).

##  VI Homotopy II

As already mentioned above, the second form of homotopy approaches we wish to
discuss is the one used to characterize the physical solutions in
$\mathcal{S}=\mathcal{S}(G)$ corresponding to truncated CC equations (as
characterized by the choice of excitation graph $G\subset G^{\mathrm{full}}$).
A homotopy scheme can then be used to connect $\mathcal{S}(G)$ to some
“higher” truncation scheme, in particular the FCC (FCI) regime (i.e., CC
method with $G^{\mathrm{full}}$).

Before continuing further, we note that ${\bf t}\in\mathcal{S}(G)$ is
equivalent to

$\langle\hat{S}({\bf s})\Phi_{0}\mid\hat{\mathcal{H}}({\bf
t})\mid\Phi_{0}\rangle=0$ (18)

for all ${\bf s}\in\mathbb{V}(G)$. Let $\mathbb{V}^{*}$ denote the dual space
of $\mathbb{V}$, we then define the CC mapping $\mathcal{A}({\bf
t}):\mathbb{V}\to\mathbb{V}^{*}$ via

$\langle{\bf s},\mathcal{A}({\bf t})\rangle=\langle\hat{S}({\bf
s})\Phi_{0}\mid\hat{\mathcal{H}}({\bf t})\mid\Phi_{0}\rangle.$ (19)

Here $\langle\cdot,\cdot\rangle$ is also used to denote the dual pairing
between amplitudes and elements in the dual amplitude space. Note that, by
definition, $\mathcal{S}$ is the set of zeros to $\mathcal{A}$. We will assume
for ${\bf t}_{*}\in\mathcal{S}$ that $\det\mathcal{A}^{\prime}({\bf
t}_{*})\neq 0$, i.e., there are no non-trivial right eigenvectors to
$\hat{\mathcal{H}}({\bf t}_{*})$ corresponding to the eigenvalue
$\mathcal{E}_{\mathrm{CC}}({\bf t}_{*})$ (in short, ${\bf t}_{*}$ is a non-
degenerate zero).

Let $\mathbb{V}^{1}$ be the amplitude space corresponding to the excitation
graph $G^{1}$. This can be thought of as the untruncated CC case, i.e., full
graph with all possible determinants (vertices) included. Similarly,
$\mathbb{V}^{0}$ is the amplitude space corresponding to $G^{0}$, which
contains amplitudes with rank $\leq\rho$ ($2\leq\rho<N$). We then decompose
$\mathbb{V}^{1}=\mathbb{V}(G^{1})$ as:
$\mathbb{V}^{1}=\mathbb{V}^{0}\oplus\mathbb{V}^{\perp}$. Here,
$\mathbb{V}^{\perp}$ contains amplitudes with rank rank strictly greater than
$\rho$. It hols $\langle{\bf s}^{0},{\bf s}^{\perp}\rangle=0$ for ${\bf
s}^{0}\in\mathbb{V}^{0}$ and ${\bf s}^{\perp}\in\mathbb{V}^{\perp}$.

###  VI.1 Existence for Kowalski–Piecuch homotopy

We will next discuss (from a more mathematical perspective) how to make sense
of a truncated CC solution by finding a trajectory that connects it to a FCC
(FCI) solution. We define the Kowalski–Piecuch homotopy
$\mathcal{K}_{\mathrm{KP}}:\mathbb{V}^{1}\times[0,1]\to(\mathbb{V}^{1})^{*}$
via the instruction

$\displaystyle\langle\mathcal{K}_{\mathrm{KP}}({\bf t}^{1},\lambda),{\bf
s}^{1}\rangle=\langle\hat{S}({\bf s}^{0})\Phi_{0}\mid\hat{\mathcal{H}}({\bf
t}^{0})\mid\Phi_{0}\rangle$ (20) $\displaystyle\quad+\langle\hat{S}({\bf
s^{\perp}})\Phi_{0}\mid\hat{\mathcal{H}}({\bf t}^{1})\mid\Phi_{0}\rangle$
$\displaystyle\quad+\lambda\langle\hat{S}({\bf
s}^{0})\Phi_{0}\mid\hat{\mathcal{H}}({\bf t}^{0})\mid(e^{\hat{T}({\bf
t}^{\perp})}-\hat{I})\Phi_{0}\rangle$

for all ${\bf t}^{1},{\bf s}^{1}\in\mathbb{V}^{1}$ and $\lambda\in[0,1]$. We
can note that by construction, the Kowalski–Piecuch homotopy satisfies
$\mathcal{K}_{\mathrm{KP}}({\bf t}^{1},1)=\mathcal{A}({\bf t}^{1})$.
Furthermore, $\mathcal{K}_{\mathrm{KP}}({\bf t}^{1}_{**},0)=0$ is equivalent
to

$\displaystyle\langle\hat{S}({\bf s}^{0})\Phi_{0}\mid\hat{\mathcal{H}}({\bf
t}_{**}^{0})\mid\Phi_{0}\rangle$ $\displaystyle=0,$ (21)
$\displaystyle\langle\hat{S}({\bf
s}^{\perp})\Phi_{0}\mid\hat{\mathcal{H}}({\bf t}_{**}^{0}+{\bf
t}_{**}^{\perp})\mid\Phi_{0}\rangle$ $\displaystyle=0.$ (22)

Equation (21) is the usual truncated CC equation (associated with $G^{0}$) and
Eq. (22) is the Kowalski–Piecuch-auxiliary equation.

The Kowalski–Piecuch homotopy has been used to establish an existence result
that connects a truncated solution to a corresponding untruncated (“full”)
solution Csirik and Laestadius (2023):

Let ${\bf t}_{*}^{1}\in\mathbb{V}^{1}$ be a non-degenerate zero of
$\mathcal{A}$. Under technical assumptions (see Theorem 4.34 in Csirik and
Laestadius (2023)), we can find $\varepsilon>0$ such that for any
$\lambda\in[0,1)$ there exists ${\bf
t}_{**}^{1}(\lambda)\in\mathcal{D}_{\varepsilon}$ fulfilling
$\mathcal{K}_{\mathrm{KP}}({\bf t}_{**}^{1}(\lambda),\lambda)=0$, where

$\mathcal{D}_{\varepsilon}=\\{{\bf t}_{*}^{1}+{\bf
r}^{1}\in\mathbb{V}^{1}:\|{\bf r}^{0}\|_{\mathbb{V}}^{2}+\|{\bf
r}^{\perp}\|_{\mathbb{V}}^{2}<\varepsilon\\}.$

In particular, there exists ${\bf t}_{**}^{1}\in\mathcal{D}_{\varepsilon}$
such that $\mathcal{A}({\bf t}_{**}^{0})=0$. (See FIG. 3 for an illustration.)

For the interested reader, the technical assumptions are discussed in Csirik
and Laestadius (2023) (in Theorem 4.34 and directly afterward). In essence,
these assumptions concern the fluctuation potential (and how it couples the
truncated space to the “rest”, which can be controlled by the truncation level
$\rho$) and the size of the ${\bf t_{**}}^{\perp}$ (that cannot be too large).
One of the assumptions can also be interpreted as a perturbative assumption
related to the size of the second derivative of the CC mapping (i.e.,
$\mathcal{A}^{\prime\prime}$) at ${\bf t}_{*}^{1}$.

$\lambda=0$$\lambda=1$$\lambda$$\qquad{\bf t}_{*}^{1}$${\bf
t}_{**}^{1}\qquad{}$$\mathcal{D}_{\varepsilon}$${\bf t}_{**}^{1}(\lambda)$
Figure 3: Visualization of the existence result of the Kowalski–Piecuch
homotopy. Under certain assumptions there is a “tube” in amplitude space that
connects a truncated CC solution to a FCC (FCI) solution. In principle, such
trajectory allows us to select which solutions of a truncated CC calculation
are “physical”.

###  VI.2 Error estimate for Kowalski–Piecuch homotopy

The Fundamental theorem of the formalism of $\beta$-nested equations of
Kowalski and Piecuch Piecuch and Kowalski (2000) can be cast into an a
posteriori error estimate Csirik and Laestadius (2023). Let
$\mathbb{V}^{1}=\mathbb{V}(G^{\mathrm{full}})$. Suppose that ${\bf
t}_{*}^{1}\in\mathbb{V}^{1}$ is a zero of $\mathcal{A}$ and ${\bf
t}_{**}^{1}\in\mathbb{V}^{1}$ a zero of
$\mathcal{K}_{\mathrm{KP}}(\cdot,0)=0$. We then have the energy error estimate
valid for any (approximate) eigenstate (and not just the ground state)
parameterized by the truncated amplitudes ${\bf t}_{**}^{0}$ that solve Eq.
(21):

Let $\rho\geq 2$ and assume that the nonorthogonality condition $\langle
e^{\hat{T}({\bf t}^{0}_{**})}\Phi_{0}\mid e^{\hat{T}({\bf
t}_{*}^{1})}\Phi_{0}\rangle\neq 0$ hold. Then

$|\mathcal{E}_{\mathrm{CC}}({\bf t}_{**}^{0})-\mathcal{E}_{\mathrm{CC}}({\bf
t}_{*}^{1})|\leq C({\bf t}_{**}^{1},{\bf t}_{*}^{1})\,\,\|{\bf
t}_{**}^{\perp}\|_{\mathbb{V}},$ (23)

where ${\bf t}_{**}^{\perp}$ solves Kowalski–Piecuch-auxiliary equation (22).

Here the constant in the rhs. of Eq. (23) is given by

$\displaystyle C({\bf t}_{**}^{1},{\bf t}_{*}^{1})=\tilde{C}({\bf
t}_{**}^{1},{\bf t}_{*}^{1})\,|\langle e^{\hat{T}({\bf
t}^{0}_{**})}\Phi_{0}\mid e^{\hat{T}({\bf t}_{*}^{1})}\Phi_{0}\rangle|^{-1}.$
(24)

The constant $\tilde{C}$ is a multifaceted quantity that comprises several
contributing factors. These factors include a norm equivalence constant, the
similarity transform of the fluctuation operator (where the Hamiltonian is a
sum of the Fock and fluctuation operator), and the size of the orthogonal
space $\mathcal{V}^{\perp}$. We emphasize that the constant becomes smaller in
the perturbative regime of the CC method when the fluctuation potential is
small and $\mathcal{V}^{\perp}$ is less physically relevant. For more details,
we refer the interested reader to Theorem 4.40 in Csirik and Laestadius
(2023).

If the nonorthogonality condition holds and ${\bf t}_{**}^{\perp}=0$, then
${\bf t}_{**}^{0}$ is a FCC (FCI) solution and the energy error estimate
implies that the error is zero, i.e., $\mathcal{E}_{\mathrm{CC}}({\bf
t}_{**}^{0})=\mathcal{E}_{\mathrm{CC}}({\bf t}_{*}^{1})$. Consequently,
$\|{\bf t}_{**}^{\perp}\|_{\mathbb{V}}$ allows us to view Eq.(22) as providing
an _a posteriori_ error estimate for a truncated CC calculation (that gives us
${\bf t}_{**}^{0}$), which could then have practical use if the
Kowalski–Piecuch-auxiliary equation can be solved (at least approximately due
to its FCI complexity).

If $\langle e^{\hat{T}({\bf t}^{0}_{**})}\Phi_{0}\mid e^{\hat{T}({\bf
t}_{*}^{1})}\Phi_{0}\rangle=0$ and ${\bf t}_{*}^{1}$ is assumed to be non-
degenerate, then $e^{\hat{T}({\bf t}^{0}_{**})}|{\Phi_{0}}\rangle$ and
$e^{\hat{T}({\bf t}_{*}^{1})}|{\Phi_{0}}\rangle$ represent different
eigenstates. This means that $e^{\hat{T}({\bf t}^{0}_{**})}|{\Phi_{0}}\rangle$
is an approximation to an eigenstate _different_ from $e^{\hat{T}({\bf
t}_{*}^{1})}|{\Phi_{0}}\rangle$. In this case, it does not make sense to try
to connect these solutions and a comparison between the energies in terms of
an energy error estimate is not meaningful (in agreement with a divergent
constant $C$). Conversely, if the energy difference
$|\mathcal{E}_{\mathrm{CC}}({\bf t}_{**}^{0})-\mathcal{E}_{\mathrm{CC}}({\bf
t}_{*}^{1})|$ diverges for finite cluster amplitudes, then $\langle
e^{\hat{T}({\bf t}^{0}_{**})}\Phi_{0}\mid e^{\hat{T}({\bf
t}_{*}^{1})}\Phi_{0}\rangle\to 0$ since we are trying to connect solutions at
different truncation levels that are orthogonal (without degeneracies this
means solutions belonging to different eigenvalues).

In the beginning of this section we referenced the origin of the energy error
estimate above to the the Fundamental theorem of the formalism of
$\beta$-nested equation due to Kowalski and Piecuch Piecuch and Kowalski
(2000). This theorem offers an explicit formula for the noniterative
correction required to supplement the energy obtained from an (approximate)
SRCC calculations, e.g., on the CCSD level, in order to retrieve the FCI
result (see Eq. (6) in Kowalski and Piecuch (2000c)). The approach of
addressing the many-electron correlation problem, based on this formula, is
called the _method of moments of coupled-cluster equations_ (MMCC) Piecuch and
Kowalski (2000); Kowalski and Piecuch (2000c) In this regard the
Kowalski–Piecuch homotopy has also helped the development of the completely
renormalized (CR) CC and EOM-CC methods and generalizations that exist under
the name CC$(P;Q)$ (see Shen and Piecuch (2012) and references therein).

##  VII Discussion

In this manuscript, we elaborated on recent advances in applying homotopy
methods to CC theory. We introduced the concept of homotopy continuation
methods, which have proven successful in solving polynomial systems in the
wider applied mathematics community. We then delved into the challenges faced
by state-of-the-art approximation techniques when solving the CC equations in
a non-perturbative regime, and we discussed how homotopy continuation methods
can potentially overcome these obstacles.

Despite the promising benefits and previous research by pioneers such as
Živković and Monkhorst Živković and Monkhorst (1978), as well as Kowalski and
Piecuch Kowalski and Jankowski (1998a); Piecuch and Kowalski (2000); Kowalski
and Piecuch (2000a, b), the high dimensionality, non-linearity, and steep
scaling of algebro-computational methods have hindered widespread adoption of
homotopy continuation methods in CC theory. However, significant progress has
been made in scaling algebro-computational methods Verschelde (1999); Bates
_et al._ (2006); Chen  _et al._ (2014); Lee  _et al._ (2008); Breiding and
Timme (2018). We are convinced that in the near future, these approaches will
be adopted and further developed to be applied to CC theory. In this
manuscript, we present two exciting use cases where homotopy methods can make
a significant impact.

The first application involves using homotopy methods to compute all roots of
the CC equations Faulstich and Oster (2022). We outline the advantages of
having access to all roots in CC theory, and how this can extend CC methods
beyond an application in the perturbative regime. The main challenge is that,
despite the great progress in the field, the scaling of algebro-computational
methods is still unfavorable for high-dimensional problems like the CC
equations. However, computational algebraic tools are developed to find the
roots of polynomial systems in the most general sense, in particular, they do
not exploit the structured nature and symmetries of the CC equations. These
provide opportunities to overcome unfavorable scaling, which is the subject of
current investigations. In this article, we outlined the first step for
homotopy approaches, which is an accurate estimate of the number of roots. We
here build upon results by Kowalski and Piecuch Piecuch and Kowalski (2000)
that reduced the estimated number of roots to $3^{n_{s}}4^{n_{d}}$. We
elaborate on how to reduce the order of the CC polynomials to be of merely
second order by extending the CC equations to an algebraic variety
$\tilde{\mathcal{S}}$. The obtained bound is $2^{n_{s}+2n_{d}}$, which
provides a significant improvement to the existing bound.

The second use case involves homotopy methods to continue CC solutions from
one level of accuracy to another, in particular, connecting solutions from a
truncated CC scheme to the FCI regime. This specific use of homotopy is, in
principle, a very important tool to justify what elements of $\mathcal{S}(G)$
constitute physical solutions that approximate eigenstates. The
Kowalski–Piecuch homotopy (with parameter $\lambda\in[0,1]$) allows us to
create trajectories that connect solutions at different levels of accuracy and
we have in this article discussed a mathematical existence result that
guarantees such curves as a function of $\lambda$ (under very technical
assumptions that we have here omitted). Furthermore, building on the
$\beta$-nested equations of Kowalski and Piecuch, an a posteriori energy error
estimate has previously been derived. This result is not restricted to the
ground state. It would be interesting if future work could reveal more about
the role of the Kowalski–Piecuch-auxiliary equations that might be of
practical use. It is noteworthy how these auxiliary equations bear a
resemblance to the tailored CC method Kinoshita  _et al._ (2005); Hino  _et
al._ (2006) that might provide further insights.

##  VIII Acknowledgment

AL acknowledges funding through ERC StG REGAL under agreement No. 101041487
and RCN (the Research Council of Norway) under agreement 287906 (CCerror) and
262695 (Hylleraas Centre for Quantum Molecular Sciences). F.M.F acknowledges
the funding from the Air Force Office of Scientific Research under award
number FA9550-18-1-0095 and from the Simons Targeted Grants in Mathematics and
Physical Sciences on Moiré Materials Magic.

## References

  * Bartlett and Musiał (2007) R. J. Bartlett and M. Musiał, Reviews of Modern Physics 79, 291 (2007).
  * Coester (1958) F. Coester, Nuclear Physics 7, 421 (1958).
  * Hubbard (1957) J. Hubbard, Proc. R. Soc. Lond. A 240, 539 (1957).
  * Hugenholtz (1957) N. Hugenholtz, Physica 23, 533 (1957).
  * Čížek (1966) J. Čížek, The Journal of Chemical Physics 45, 4256 (1966).
  * Kümmel (1991) H. Kümmel, Theoretica chimica acta 80, 81 (1991).
  * Čížek (1991) J. Čížek, Theoretica chimica acta 80, 91 (1991).
  * Bartlett (2005) R. Bartlett, in  _Theory and Applications of Computational Chemistry (The First Forty Years)_ , edited by C. Dykstra, G. Frenking, K. Kim, and G. Scuseria (2005) 1st ed.
  * Paldus (2005) J. Paldus, in  _Theory and Applications of Computational Chemistry (The First Forty Years)_ , edited by C. Dykstra, G. Frenking, K. Kim, and G. Scuseria (2005) 1st ed.
  * Arponen (1991) J. S. Arponen, Theoretica chimica acta 80, 149 (1991).
  * Bishop (1991) R. Bishop, Theoretica chimica acta 80, 95 (1991).
  * Press  _et al._ (2007) W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,  _Numerical recipes 3rd edition: The art of scientific computing_ (Cambridge university press, 2007).
  * Rowe (1968) D. J. Rowe, Rev. Mod. Phys. 40, 153 (1968).
  * Monkhorst (1977) H. J. Monkhorst, International Journal of Quantum Chemistry 12, 421 (1977).
  * Koch  _et al._ (1990) H. Koch, H. J. A. Jensen, P. Jörgensen, and T. Helgaker, The Journal of Chemical Physics 93, 3345 (1990), https://pubs.aip.org/aip/jcp/article-pdf/93/5/3345/11261510/3345_1_online.pdf .
  * Stanton and Bartlett (1993) J. F. Stanton and R. J. Bartlett, The Journal of Chemical Physics 98, 7029 (1993), https://pubs.aip.org/aip/jcp/article-pdf/98/9/7029/11309474/7029_1_online.pdf .
  * Koch  _et al._ (1994) H. Koch, R. Kobayashi, A. Sanchez de Merás, and P. Jo/rgensen, The Journal of Chemical Physics 100, 4393 (1994), https://pubs.aip.org/aip/jcp/article-pdf/100/6/4393/13970649/4393_1_online.pdf .
  * Živković and Monkhorst (1978) T. P. Živković and H. J. Monkhorst, Journal of Mathematical Physics 19, 1007 (1978).
  * Kowalski and Jankowski (1998a) K. Kowalski and K. Jankowski, Physical Review Letters 81, 1195 (1998a).
  * Paldus  _et al._ (1993) J. Paldus, P. Piecuch, L. Pylypow, and B. Jeziorski, Physical Review A 47, 2738 (1993).
  * Piecuch  _et al._ (1990) P. Piecuch, S. Zarrabian, J. Paldus, and J. Čížek, Physical Review B 42, 3351 (1990).
  * Piecuch and Kowalski (2000) P. Piecuch and K. Kowalski, in  _Computational chemistry: reviews of current trends_ (World Scientific, 2000) pp. 1–104.
  * Kowalski and Piecuch (2000a) K. Kowalski and P. Piecuch, Physical Review A 61, 052506 (2000a).
  * Kowalski and Piecuch (2000b) K. Kowalski and P. Piecuch, International Journal of Quantum Chemistry 80, 757 (2000b).
  * Kowalski and Jankowski (1998b) K. Kowalski and K. Jankowski, Chemical physics letters 290, 180 (1998b).
  * Jankowski and Kowalski (1999a) K. Jankowski and K. Kowalski, The Journal of chemical physics 110, 3714 (1999a).
  * Jankowski and Kowalski (1999b) K. Jankowski and K. Kowalski, The Journal of chemical physics 110, 9345 (1999b).
  * Jankowski and Kowalski (1999c) K. Jankowski and K. Kowalski, The Journal of chemical physics 111, 2940 (1999c).
  * Jankowski and Kowalski (1999d) K. Jankowski and K. Kowalski, The Journal of chemical physics 111, 2952 (1999d).
  * Schneider (2009) R. Schneider, Numerische Mathematik 113, 433 (2009).
  * Rohwedder (2013) T. Rohwedder, ESAIM: Mathematical Modelling and Numerical Analysis 47, 421 (2013).
  * Rohwedder and Schneider (2013) T. Rohwedder and R. Schneider, ESAIM: Mathematical Modelling and Numerical Analysis 47, 1553 (2013).
  * Laestadius and Kvaal (2018) A. Laestadius and S. Kvaal, SIAM Journal on Numerical Analysis 56, 660 (2018).
  * Faulstich  _et al._ (2018) F. M. Faulstich, A. Laestadius, S. Kvaal, Ö. Legeza, and R. Schneider, arXiv preprint arXiv:1802.05699 (2018).
  * Laestadius and Faulstich (2019) A. Laestadius and F. M. Faulstich, Molecular Physics 117, 2362 (2019).
  * Faulstich  _et al._ (2019) F. M. Faulstich, M. Máté, A. Laestadius, M. A. Csirik, L. Veis, A. Antalik, J. Brabec, R. Schneider, J. Pittner, S. Kvaal, _et al._ , Journal of chemical theory and computation 15, 2206 (2019).
  * Kvaal  _et al._ (2020) S. Kvaal, A. Laestadius, and T. Bodenstein, Molecular Physics 118, e1810349 (2020), https://doi.org/10.1080/00268976.2020.1810349 .
  * Faulstich  _et al._ (2023) F. M. Faulstich, H. E. Kristiansen, M. A. Csirik, S. Kvaal, T. B. Pedersen, and A. Laestadius, arXiv:2301.11393 (2023).
  * Csirik and Laestadius (2023) M. A. Csirik and A. Laestadius, ESAIM: Mathematical Modelling and Numerical Analysis 57, 545 (2023).
  * Faulstich and Oster (2022) F. M. Faulstich and M. Oster, arXiv:2211.10389 [to appear in SIAM Journal on Applied Algebra and Geometry] (2022).
  * Helgaker  _et al._ (2014) T. Helgaker, P. Jorgensen, and J. Olsen,  _Molecular Electronic-Structure Theory_ (John Wiley & Sons, 2014).
  * Csirik, Mihály A. and Laestadius, Andre (2023) Csirik, Mihály A. and Laestadius, Andre, ESAIM: M2AN 57, 645 (2023).
  * Bates  _et al._ (2023) D. J. Bates, P. Breiding, T. Chen, J. D. Hauenstein, A. Leykin, and F. Sottile, arXiv preprint arXiv:2302.08585 (2023).
  * Garcia and Zangwill (1979) C. B. Garcia and W. I. Zangwill, Mathematical Programming 16, 159 (1979).
  * Morgan (2009) A. Morgan,  _Solving polynomial systems using continuation for engineering and scientific problems_ (SIAM, 2009).
  * Sommese  _et al._ (2005) A. J. Sommese, C. W. Wampler, _et al._ ,  _The Numerical solution of systems of polynomials arising in engineering and science_ (World Scientific, 2005).
  * Davidenko (1953a) D. Davidenko, in  _Dokl. Akad. Nauk SSSR_ , Vol. 88 (1953) pp. 601–602.
  * Davidenko (1953b) D. Davidenko, Ukr. Mat. Zh 5, 196 (1953b).
  * Morgan  _et al._ (1992) A. P. Morgan, A. J. Sommese, and C. W. Wampler, Advances in Applied Mathematics 13, 305 (1992).
  * Verschelde (1999) J. Verschelde, ACM Transactions on Mathematical Software (TOMS) 25, 251 (1999).
  * Bates  _et al._ (2006) D. J. Bates, J. D. Hauenstein, A. J. Sommese, and C. W. Wampler, “Bertini: Software for numerical algebraic geometry,” (2006).
  * Chen  _et al._ (2014) T. Chen, T.-L. Lee, and T.-Y. Li, in  _Mathematical Software–ICMS 2014: 4th International Congress, Seoul, South Korea, August 5-9, 2014. Proceedings 4_ (Springer, 2014) pp. 183–190.
  * Lee  _et al._ (2008) T.-L. Lee, T.-Y. Li, and C.-H. Tsai, Computing 83, 109 (2008).
  * Breiding and Timme (2018) P. Breiding and S. Timme, in  _Mathematical Software–ICMS 2018: 6th International Conference, South Bend, IN, USA, July 24-27, 2018, Proceedings 6_ (Springer, 2018) pp. 458–465.
  * Hubbard  _et al._ (2001) J. Hubbard, D. Schleicher, and S. Sutherland, Inventiones mathematicae 146, 1 (2001).
  * Schleicher (2002) D. Schleicher, Ergodic Theory and Dynamical Systems 22, 935 (2002).
  * Eriksen and Gauss (2019) J. J. Eriksen and J. Gauss, Journal of Chemical Theory and Computation 15, 4873 (2019).
  * Kowalski and Piecuch (2000c) K. Kowalski and P. Piecuch, The Journal of Chemical Physics 113, 18 (2000c), https://pubs.aip.org/aip/jcp/article-pdf/113/1/18/10819898/18_1_online.pdf .
  * Shen and Piecuch (2012) J. Shen and P. Piecuch, The Journal of Chemical Physics 136 (2012), 10.1063/1.3700802, 144104, https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.3700802/15450087/144104_1_online.pdf .
  * Kinoshita  _et al._ (2005) T. Kinoshita, O. Hino, and R. J. Bartlett, The Journal of chemical physics 123, 074106 (2005).
  * Hino  _et al._ (2006) O. Hino, T. Kinoshita, G. K.-L. Chan, and R. J. Bartlett, The Journal of chemical physics 124, 114311 (2006).


# A topological behavior of quantum particles originated from the classicality
of their flow velocity

Tomer Shushi
Center for Quantum Science and Technology
& Department of Business Administration,
Guilford Glazer Faculty of Business and Management,
Ben-Gurion University of the Negev, Beer-Sheva, Israel

###### Abstract

In this letter, we propose a new quantum effect that naturally emerges from
describing the quantum particle as a classical fluid. Following the
hydrodynamical formulation of quantum mechanics for a particle in a finite
convex region, we show how the maximum values of the wavefunction’s amplitude
lie along the boundaries of the region when imposing a vanished quantum
potential, implying a classical flow velocity of the particle. The effect is
obtained for the case of particles in curved space, described by Riemannian
structures. We further show that such an effect cannot be achieved in the
relativistic regime when dealing with quantum particles in flat or curved
spacetime.

Keywords: classicality, curved space, Madelung equations, quantum
hydrodynamics, Riemannian structures

##  1
Introduction

Quantum mechanics is one of the most successful and profound theories in
modern science, offering an incredibly precise framework for understanding the
behavior of particles at the microscopic level. The Schrödinger equation
stands as a cornerstone in quantum mechanics, providing a mathematical
description of the behavior of particles as wavefunctions. However, there are
different formulations of the quantum particle that provide different
interpretations about the nature of the quantum particles. Only a year after
Erwin Schrödinger published his famous paper [1], which proposed the
Schrödinger equation as the description for the dynamics of the quantum
particles, Erwin Madelung published another paper [2], showing that the
Schrödinger equation can be converted into a dual problem that describes a
quantum fluid. Consider the Schrödinger equation of a non-relativistic
particle with mass $m$ and some external potential $V,$
$i\hbar\partial_{t}\psi=-\frac{\hbar^{2}}{2m}\nabla^{2}\psi+V\psi.$ Assuming
that $\psi$ is a smooth function, by taking the polar presentation
$\psi\left(\boldsymbol{x},t\right)=\sqrt{\rho\left(\boldsymbol{x},t\right)}e^{iS\left(\boldsymbol{x},t\right)/\hbar},$
we obtain the Madelung hydrodynamical equations for the quantum particle,

$\partial_{t}\rho+\nabla\cdot\left(\rho\boldsymbol{u}\right)=0,$ (1)

and

$\partial_{t}\boldsymbol{u}+\boldsymbol{u}\cdot\nabla\boldsymbol{u}=-\frac{1}{m}\nabla\left(Q+V\right).$
(2)

Here $\rho$ is the density function of the particle, and
$\boldsymbol{u=}\nabla S/m$ is the flow velocity of the particle, where (1) is
the continuity equation of the fluid, and (2) is the Hamilton-Jacobi equation,
with the addition of the quantum potential

$Q\left(\sqrt{\rho}\right)=-{\frac{\hbar^{2}}{2m}}{\frac{\nabla^{2}\sqrt{\rho}}{\sqrt{\rho}}.}$
(3)

We thus see that by reformulating the Schrödinger equation into a
hydrodynamical reformulation, another term pops up in the description of the
particle, which is the quantum potential. The quantum potential $Q$ gives the
coupling into the hydrodynamical equations (1)-(2), in the sense that the flow
velocity $\boldsymbol{u}$ is coupled to the wavefunction’s amplitude through
the quantum potential. We note that unlike the Schrödinger equation, here, the
description of the particle as a fluid excludes the notion of potentials since
the equations of motion describe forces instead of potentials. The external
force is $F_{\text{ext}}=-\nabla V,$ while the quantum force is defined by
$F_{\text{quant}}=-\nabla Q.$

Since the introduction of quantum hydrodynamics by Madelung, a large amount of
literature has been developed in order to both study the foundations of
quantum mechanics and also to explore novel phenomena within quantum systems
[3-9], and in recent years, quantum hydrodynamics has been gaining more
attention with a growing amount of research in the field. When taking the
classical limit $\hbar\rightarrow 0$ or the case of a massive particle
$m\rightarrow\infty,$ the quantum potential vanishes, and so (1)-(2) describes
a classical fluid, where the flow velocity is independent on the density
function. Playing with the coefficient of $Q$ is a trivial way to achieve this
sort of classicality. However, we can achieve a non-trivial classicality by
imposing a vanished quantum force $F_{\text{quant}}=0,$ which boils down to
the equation

$Q\left(\sqrt{\rho}\right)=C\left(t\right)$ (4)

where $C$ is some real-valued function of time (see, [10]). This framework
shows us that for a suitable shape of the wavefunction, followed by its
density function $\rho,$ we can achieve a description of the particle as a
classical fluid.

In the following, we propose a topological effect that emerges from quantum
particles with a vanished quantum force for particles in Riemannian structures
that describes curved space. We also show that this topological behavior does
not occur for particles in curved spacetime.

##  2 Results

We start with considering a non-relativistic quantum particle in a Riemannian
structure with the metric $ds^{2}=g_{ij}\left(x\right)dx^{i}dx^{j}.$ We assume
that our quantum particle exists in an open subset $M\subset\mathbb{R}^{n}$ of
of the Riemannian space for $n>1$ spatial dimensions. The Schrödinger equation
is then given by

$i\hbar\partial_{t}\psi=-\frac{\hbar^{2}}{2m}\nabla^{2}\psi+V\psi,$ (5)

where $\Delta=\nabla^{2}$ is the Laplace-Beltrami operator. The Hamiltonian
function for the Schrödinger equation is (see, [5])

$H\left(\psi\right)=\frac{\hbar^{2}}{2m}\left\|\nabla\psi\right\|_{L^{2}}^{2}+\left\langle\left\langle
V\psi,\psi\right\rangle\right\rangle_{L^{2}}.$ (6)

Following the transformation $\psi\longmapsto\left(\rho,S\right)$,
$\psi:=\sqrt{\rho}e^{iS/\hbar}$, we can express $H\left(\psi\right)$ in terms
of $\left(\rho,S\right).$ Following
$\frac{\hbar^{2}}{2m}\left\|\nabla\psi\right\|_{L^{2}}^{2}=\frac{\hbar^{2}}{2m}\left\langle\left\langle\nabla\sqrt{\rho},\nabla\sqrt{\rho}\right\rangle\right\rangle_{L^{2}}+\frac{1}{2m}\left\langle\left\langle\rho\nabla
S,\nabla S\right\rangle\right\rangle_{L^{2}},$ the Madelung Hamiltonian is
given by

$H\left(\rho,S\right)=\frac{1}{2m}\left\langle\left\langle\rho\nabla S,\nabla
S\right\rangle\right\rangle_{L^{2}}+\frac{\hbar^{2}}{2m}\left\langle\left\langle\nabla\sqrt{\rho},\nabla\sqrt{\rho}\right\rangle\right\rangle_{L^{2}}+\left\langle\left\langle
V\rho,1\right\rangle\right\rangle_{L^{2}}.$ (7)

The Madelung equations are then obtained using Hamilton’s equations (see,
again [5])

$\partial_{t}S=-\frac{\delta H\left(\rho,S\right)}{\delta\rho},\text{ \
}\partial_{t}\rho=\frac{\delta H\left(\rho,S\right)}{\delta S},$ (8)

leading to

$\partial_{t}\boldsymbol{u}+\boldsymbol{u}\cdot\nabla\boldsymbol{u}+\frac{1}{m}\nabla\left(Q_{g}+V\right)=0,$
(9)

and

$\partial_{t}\rho+\operatorname{div}\left(\rho\boldsymbol{u}\right)=0,$ (10)

where $\boldsymbol{u}:=\nabla S/m$ is the particle’s flow velocity. The
quantum potential is then given by

$Q_{g}\left(\sqrt{\rho}\right)=-\frac{\hbar^{2}}{2m}\frac{\Delta\sqrt{\rho}}{\sqrt{\rho}}=-\frac{\hbar^{2}}{2m}\frac{\frac{1}{\sqrt{g}}\partial_{x^{i}}\left(\sqrt{g}g^{ij}\partial_{x^{j}}\sqrt{\rho}\right)}{\sqrt{\rho}},$
(11)

where $\Delta$ is the Laplace-Beltrami operator, which also depends on the
metric $g.$ Similar to the case of Euclidean space, here, also, the quantum
potential can trivially vanish when taking the limits $\hbar\rightarrow 0$ or
$m\rightarrow\infty.$

From the quantum potential $Q_{g}$ arises the quantum force

$F_{Q_{g}}=-\nabla
Q_{g}=\frac{\hbar^{2}}{2m}\frac{\nabla\frac{1}{\sqrt{g}}\partial_{x^{i}}\left(\sqrt{g}g^{ij}\partial_{x^{j}}\sqrt{\rho}\right)-\frac{1}{2\rho\sqrt{g}}\partial_{x^{i}}\left(\sqrt{g}g^{ij}\partial_{x^{j}}\sqrt{\rho}\right)\nabla\rho}{\sqrt{\rho}}.$
(12)

We produce a classical fluid by canceling the quantum force, with imposing
$F_{Q_{g}}=0,$

$Q_{g}=C\left(t\right),$ (13)

for some real-valued time-dependent function $C\left(t\right)$.
Then, equation (13) boils down to

$g^{ij}\partial_{x^{i}}\partial_{x^{j}}P+\partial_{x^{i}}\left(\sqrt{g}g^{ij}\right)\partial_{x^{j}}P+\frac{2m}{\hbar^{2}}C\left(t\right)P=0.$
(14)

Recalling that $g^{ij}$ is symmetric and positive-definite, and assume that
$g^{ij}$ and $\partial_{x^{i}}\left(\sqrt{g}g^{ij}\right)\partial_{x^{j}}P$
have smooth (bounded) components, when setting $C\left(t\right)\equiv 0,$
equation (13) satisfies the strong maximum principle (smp).

The smp guarantees that the maximum value of the wavefunction’s amplitude $P$
will be on the boundary $\partial M$ of the finite region the particle exists
in, and moreover, these maxima only exist on the boundary. For non-negative
solutions $P_{\text{class}}$ of (14), this also implies that the density
function $\rho^{\ast}=P_{\text{class}}^{2}$ gains its maximum on the boundary,
creating a topological effect of the particle which then moves along the
boundary of the region it exists in.

The following is an illustrative example of the proposed topological effect in
the case of two spatial dimensions for a conformally flat curved space.

Figure 1. The wavefunction’s amplitude $P_{\text{class}}$ for a two-
dimensional, $(x,y),$ quantum particle in a conformally curved space, with
$g^{ij}=\Omega\left(x,y\right)\cdot diag\left(1,1\right),\mu=2.5,$ the
conformal factor
$\Omega\left(x,y\right)=\left(e^{-\left(x-\mu\right)^{2}}+e^{-\left(y-\mu\right)^{2}}\right),$
for a disc-shaped region (a) and different convex-shaped regions (b)-(c).

As can be seen from Figure 1, $P_{\text{class}}$ is concentrated at the
boundaries of the convex closed regions. However, we note that the proposed
effect does not, in general, imply that $P_{\text{class}}$ will mainly exist
on the boundaries, and it can also flow into the closed region.

The desired wavefunction’s amplitude, $P_{\text{class}},$ that gives the
topological effect is, in general, not a solution of the Madelung equations
(9)-(10). We, thus, have to find suitable quantum systems in which (14) can be
satisfied. To do that, we consider the following procedure (see [8]): At the
first stage, we substitute the density function
$\rho^{\ast}=P_{\text{class}}^{2}$ corresponding to $P_{\text{class}}$ into
the continuity equation (10),
$\partial_{t}\rho^{\ast}+\operatorname{div}\left(\rho^{\ast}\boldsymbol{u}\right)=0,$
which allows us to find the flow velocity $\boldsymbol{u}^{\ast}$ in which
this equation is satisfied. The second stage is to substitute both
$\rho^{\ast}$ and $\boldsymbol{u}^{\ast}$ into the second Madelung equation
(9), which then gives the desired external force $F_{\text{exp}}=-\nabla V$
that should be imposed in order to have the desired shape of the
wavefunction’s amplitude, with

$F_{\text{exp}}=m\left(\partial_{t}\boldsymbol{u^{\ast}+u^{\ast}}\cdot\nabla\boldsymbol{u^{\ast}}\right)+\nabla
Q_{g}(P_{\text{class}}).$ (15)

In the hydrodynamical formulation of quantum mechanics,
$\boldsymbol{J=}\rho\boldsymbol{u}$ is defined as the current density of the
function, and so the continuity equation (10),
$\partial_{t}\rho+\operatorname{div}\left(\boldsymbol{J}\right)=0$ essentially
describes the conservation of probability $\rho$ in the system. We note that
in case the current density does not have turbulence behavior, which is
manifested by a zero curl, $\nabla\times\boldsymbol{J=0},$ and assuming that
$\boldsymbol{J}$ goes to zero at the limits
$\left|x_{j}\right|\rightarrow\infty,$ $j=1,2,...,N,$ then we can write
$\boldsymbol{J}$ as the divergence of a scalar function $\phi,$
$\boldsymbol{J=}\nabla\phi.$ This means that we can write the continuity
equation as a Poisson’s equation, $\Delta\phi=-\partial_{t}\rho^{\ast}$. In
the case of flat space, the solution takes an explicit form, leading to the
flow velocity
$\boldsymbol{u=J}/\rho=\frac{1}{P_{\text{class}}^{2}\left(\boldsymbol{x},t\right)}\nabla\int_{\mathbb{R}^{N}}\frac{\partial_{t}P_{\text{class}}^{2}\left(\boldsymbol{x}^{\prime},t\right)}{4\pi\left|\boldsymbol{x}-\boldsymbol{x}^{\prime}\right|}d\boldsymbol{x}^{\prime}.$

In the following, we show that, unlike the case of curved space, when we are
dealing with a quantum particle in curved spacetime, such a classicality of
the quantum particle does not bring the proposed topological behavior. This
feature is directly rooted in the definition of the spacetime metric
$g^{\mu\nu}.$

###  2.1 Classicality of quantum particles in curved spacetime

Consider a relativistic spinless quantum particle in $\left(3+1\right)$ curved
spacetime, modeled by the Klein-Gordon equation in curved spacetime

$-g^{\mu\nu}\partial_{\mu}\partial_{\nu}\Psi+g^{\mu\nu}\Gamma_{\mu\nu}^{\sigma}\partial_{\sigma}\Psi+\frac{m^{2}c^{2}}{\hbar^{2}}\Psi+U\Psi=0$
(16)

where $x=\left(x_{0}=ct,x_{1},x_{2},x_{3}\right),$ and

$\Gamma_{\mu\nu}^{\sigma}=\frac{1}{2}g^{\sigma\beta}\left(\partial_{\mu}g_{\beta\nu}+\partial_{\nu}g_{\beta\mu}-\partial_{\beta}g_{\mu\nu}\right)$
is the christoffel symbol.

To achieve classicality as followed by quantum hydrodynamics, we consider the
polar representation, as before, $\Psi:=\sqrt{\rho}e^{iS/\hbar},$ with noting
that while $\sqrt{\rho}$ describes the wavefunction’s amplitude, it does not
describe the density function of the quantum particle. Then, by substituting
the polar representation into the Klein-Gordon equation (16) we immediately
obtain the corresponding Madelung equations

$g^{\mu\nu}\partial_{\mu}S\cdot\partial_{\nu}S+2m_{0}Q_{g}\left(P\right)+m^{2}c^{2}+\hbar^{2}U=0,$
(17)

and

$\nabla_{\mu}J^{\mu}=0,$ (18)

with the quantum potential

$Q_{g}\left(P\right)=-\frac{\hbar^{2}}{2m_{0}}g^{\mu\nu}\frac{\partial_{\mu}\partial_{\nu}P-\Gamma_{\mu\nu}^{\sigma}\partial_{\sigma}P}{P},$
(19)

where $\rho=P^{2},$ and $J^{\mu}:=P^{2}g^{\mu\nu}\nabla_{\nu}S$ is the four-
current, and $m_{0}>0$ is a constant in units of mass. Then, the spacetime
quantum force is given by

$F_{Q}=-\nabla_{\mu}Q_{g}\left(P\right).$ (20)

To achieve classicality, we impose a vanished quantum force $F_{Q}=0$, which
boils down to

$Q_{g}\left(P\right)+\lambda=0$ (21)

for some real constant $\lambda$ (see, [6]). Now, by setting

$\lambda=m_{eff}c^{2}$

which is in units of energy, similar to the quantum potential, and taking
$m_{0}=m_{eff},$ for some effective mass $m_{eff},$ the equation (21) can then
be converted to

$-g^{\mu\nu}\partial_{\mu}\partial_{\nu}P+g^{\mu\nu}\Gamma_{\mu\nu}^{\sigma}\partial_{\sigma}P+\frac{m_{eff}^{2}c^{2}}{\hbar^{2}}P=0.$
(22)

We, thus, achieve a vanished quantum force when the wavefunction’s amplitude
follows the Klein-Gordon equation in empty curved spacetime with the same
metric as the original Klein-Gordon equation of the quantum particle, but with
some effective mass $m_{eff}$.

Suppose now that our particle is given in a finite region
$M\subset\mathbb{R}^{3+1}$. By setting $m_{eff}=0$, we have

$g^{\mu\nu}\partial_{\mu}\partial_{\nu}P-g^{\mu\nu}\Gamma_{\mu\nu}^{\sigma}\partial_{\sigma}P=0.$
(23)

While $g^{\mu\nu}$ is symmetric, which is one of the requirements for
obtaining the strong maximum principle, the components of $g^{\mu\nu}$ have
different signs, in general,

$sgn\left(g^{\mu\nu}\right)\neq
sgn\left(g^{\mu\prime\nu^{\prime}}\right),\text{ for
}\mu\neq\mu^{\prime},\nu\neq\nu^{\prime},$ (24)

and thus, the strong maximum principle cannot, in principle, be obtained for
such systems. The feature (24) is fundamental in general relativity. In fact,
when (24) is violated, we can have causality problems in the system.

##  3 Discussion

One of the most fundamental aspects of quantum mechanics is the boundary
between the quantum and classical behavior of the quantum particles. While at
the macroscopic level, quantum effects diminish and gradually give way to
classical mechanics, at the microscopic level, quantum mechanics governs the
behavior of the systems. The transition point from quantum to classical
behavior remains a subject of intense investigation and theoretical
exploration. The challenge lies in precisely defining when and how quantum
coherence dissipates, yielding to classical predictability. In this paper, we
have proposed a topological behavior for particles in Riemannian structures,
which naturally emerges when imposing the condition that their quantum
potential ultimately vanishes. We have shown a drastic difference between the
classicality of non-relativistic and relativistic systems, followed by the
proposed topological behavior. While in the non-relativistic regime, the
topological behavior holds, in the relativistic regime, such a topological
behavior does not hold, in principle, followed by the basic nature of the
spacetime metric. We propose to explore this in future research. We can extend
the results into a system containing $n>1$ coupled particles. Studying such
systems can help get new insights into the interplay between the proposed
classicality with a vanished quantum potential and the correspondence limit in
which the system containing many coupled particles transitions into a genuine
macroscopic classical system. We propose to explore it in future research.

## References

  * [1] Schrödinger, E. (1926). An undulatory theory of the mechanics of atoms and molecules. Physical review, 28, 1049.
  * [2] Madelung, E. (1927). Quantum theory in hydrodynamical form. z. Phys, 40, 322.
  * [3] Schönberg, M. (1954). On the hydrodynamical model of the quantum mechanics. Il Nuovo Cimento (1943-1954), 12, 103-133.
  * [4] Tsekov, R., Heifetz, E., & Cohen, E. (2018). Relating quantum mechanics with hydrodynamic turbulence. Europhysics Letters, 122, 40002.
  * [5] Khesin, B., Misiolek, G., & Modin, K. (2018). Geometric hydrodynamics via Madelung transform. Proceedings of the National Academy of Sciences, 115, 6165-6170.
  * [6] Shushi, T. (2023). Classicality of single quantum particles in curved spacetime through the hydrodynamical reformulation of quantum mechanics. Journal of Physics A: Mathematical and Theoretical, 56, 365301.
  * [7] Bonitz, M., Moldabekov, Z. A., & Ramazanov, T. S. (2019). Quantum hydrodynamics for plasmas—Quo vadis?. Physics of Plasmas, 26, 090601.
  * [8] Gay-Balmaz, F., & Tronci, C. (2020). Madelung transform and probability densities in hybrid quantum–classical dynamics. Nonlinearity, 33, 5383.
  * [9] Reddiger, M., & Poirier, B. (2023). Towards a mathematical theory of the Madelung equations: Takabayasi’s quantization condition, quantum quasi-irrotationality, weak formulations, and the Wallstrom phenomenon. Journal of Physics A: Mathematical and Theoretical, 56, 193001.
  * [10] Shushi, T. (2023). An appearance of classical matter from the self-organizing process of quantum systems. Europhysics Letters, 141, 64001.