# Quantum computer formulation of the FKP-operator eigenvalue problem for probabilistic learning on manifolds

Christian Soize, Loïc Joubert-Doriol, Artur Izmaylov

► To cite this version:

Christian Soize, Loïc Joubert-Doriol, Artur Izmaylov. Quantum computer formulation of the FKP-operator eigenvalue problem for probabilistic learning on manifolds. Computer Methods in Applied Mechanics and Engineering, 2025, 443, pp.118080. 10.1016/j.cma.2025.118080. hal-05078073

HAL Id: hal-05078073

https://univ-eiffel.hal.science/hal-05078073v1

Submitted on 8 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Quantum computer formulation of the FKP-operator eigenvalue problem for probabilistic learning on manifolds

Christian Soize<sup>a,*</sup>, Loïc Joubert-Doriol<sup>a</sup>, Artur F. Izmaylov<sup>b</sup>

<sup>a</sup>Université Gustave Eiffel, MSME UMR 8208, 5 bd Descartes, 77454 Marne-la-Vallée, France

<sup>b</sup>University of Toronto Scarborough, Department of Physical and Environmental Sciences, 1265 Military Trail Toronto, Canada

## Abstract

We present a quantum computing formulation to address a challenging problem in the development of probabilistic learning on manifolds (PLoM). It involves solving the spectral problem of the high-dimensional Fokker-Planck (FKP) operator, which remains beyond the reach of classical computing. Our ultimate goal is to develop an efficient approach for practical computations on quantum computers. For now, we focus on an adapted formulation tailored to quantum computing. The methodological aspects covered in this work include the construction of the FKP equation, where the invariant probability measure is derived from a training dataset, and the formulation of the eigenvalue problem for the FKP operator. The eigen equation is transformed into a Schrödinger equation with a potential $V$, a non-algebraic function that is neither simple nor a polynomial representation. To address this, we propose a methodology for constructing a multivariate polynomial approximation of $V$, leveraging polynomial chaos expansion within the Gaussian Sobolev space. This approach preserves the algebraic properties of the potential and adapts it for quantum algorithms. The quantum computing formulation employs a finite basis representation, incorporating second quantization with creation and annihilation operators. Explicit formulas for the Laplacian and potential are derived and mapped onto qubits using Pauli matrix expressions. Additionally, we outline the design of quantum circuits and the implementation of measurements to construct and observe specific quantum states. Information is extracted through quantum measurements, with eigenstates constructed and overlap measurements evaluated using universal quantum gates.

**Keywords:** Quantum computing formulation, Fokker-Planck operator, eigenvalue problem, second quantization, probabilistic learning on manifolds, PLoM

## 1. Introduction

We present a quantum computing methodology to address a challenging problem recently encountered in the development of probabilistic learning on manifolds [1]. This approach aims to construct statistical surrogate models from small training datasets, which are particularly useful for uncertainty quantification, nonconvex optimization under uncertainty, and statistical inverse problems. The main challenge lies in solving the spectral problem associated with the high-dimensional Fokker-Planck (FKP) operator, a task currently beyond the capabilities of classical computing. In this paper, we explore a quantum computing-based methodology for calculating the eigenfunctions associated with the smallest eigenvalues of the FKP operator. The ultimate goal is to develop an efficient approach that, in the near future, could enable effective computation on a quantum computer. For now, we focus on an adapted formulation tailored to quantum computing. It is worth noting that this problem may also be relevant to frameworks and applications beyond those that motivated the developments presented in this work.

<footnote>*Corresponding author: C. Soize, christian.soize@univ-eiffel.fr</footnote>

<footnote>Email addresses: christian.soize@univ-eiffel.fr (Christian Soize), loic.joubert-doriol@univ-eiffel.fr (Loïc Joubert-Doriol), artur.izmaylov@utoronto.ca (Artur F. Izmaylov)</footnote>
*Quantum Computing Applications in Computational Mechanics and Engineering.* Since research on quantum information and quantum computation began in the 1980s, numerous papers have been published, covering all aspects necessary for developing quantum computing techniques to address large-scale simulations in engineering science, including data science. This effort requires the development of suitable quantum computing algorithms, software advancements, quantum error correction methods, and the definition of hardware and infrastructure requirements for executing large-scale numerical simulations. These challenges are particularly complex and demand interdisciplinary collaboration among quantum physicists, mathematicians, computer scientists, and engineers. Given the vast literature in this multidisciplinary field, this short introduction does not aim to review all existing works. However, the engineering science community has recently shown interest in leveraging quantum computers for numerical simulations and algorithm development. For example, quantum computing applications have been explored in computational mechanics [2, 3, 4], materials science [5, 6, 7], finite element methods [8], and engineering simulation and optimization [9, 10].

*Objectives and organization of the paper.* The aim of this paper is to propose a method to solve the spectral problem of the high-dimensional Fokker-Planck (FKP) operator, which remains beyond the capabilities of classical computing. It is known (see, for example, [11, 12]) that the FKP operator can be reformulated as a Schrödinger-type operator exhibiting a potential $V$. This reformulation bears a strong resemblance to the vibrational spectral problem in molecular systems, where quantum computing has recently gained traction as an appealing alternative to classical computers [13, 14]. We follow a similar path in the current paper, also aligning with recent efforts to implement the dynamical solution of the FKP problem on quantum computers [15]. Our approach presents one path to solve the spectral problem on quantum computers, here the Variational Quantum Eigensolver [16]. Nonetheless, there are other known approaches such as the Quantum Phase Estimation [17] (similar to the Fourier transform of the time-dependent Fokker-Planck Equation). The FKP equation under consideration is constructed within the framework of probabilistic learning on manifolds (PLoM) [1], as we will explain in Section 2. The potential $V$, associated with a probability density function related to the training dataset, is a non-algebraic function that is neither simple nor a polynomial representation, and therefore does not satisfy the algebraic prerequisites of the quantum algorithm we will consider. In the first part of this paper, we present in Sections 3 and 4 the construction of a representation of the potential of the Schrödinger operator. To achieve this, we propose a methodology for constructing a multivariate polynomial approximation [13, 14] of $V$, leveraging polynomial chaos expansion within the Gaussian Sobolev space [18]. This approach preserves the algebraic properties of the potential and adapts it for quantum algorithms. It enables the construction of a polynomial chaos expansion in Gaussian space with good convergence properties, preserving the algebraic characteristics that control the spectrum of the operator. Section 5 addresses the numerical validation of the computation of the polynomial chaos expansion coefficients of the potential. The second part of the paper, Sections 6 to 8, focuses on the quantum computing formulation. This formulation employs a finite basis representation (spectral method) that incorporates second quantization with creation and annihilation operators. Explicit formulas for the Laplacian and potential are derived and mapped onto qubits using Pauli matrix expressions. The design of quantum circuits and the implementation of measurements to construct and observe specific quantum states are presented in Section 7. Finally, Section 8 addresses the construction of eigenstates and the extraction of overlaps through quantum measurements, evaluated using universal quantum gates.

*Key steps in the proposed methodology.* Hereinafter, we present methodological aspects. Although presented in the context of the FKP-operator eigenvalue problems, the proposed methodology illustrates well the main steps required to use quantum computers.

(i) Construction of the FKP equation whose invariant probability measure, associated with the steady-state solution, is derived from a given training dataset used within the framework of probabilistic learning on manifolds.

(ii) Formulation of the eigenvalue problem related to the FKP operator and transforming it into a Schrödinger-type operator exhibiting a potential $y \mapsto V(y)$ on $\mathbb{R}^v$ with high dimension $v$.

(iii) Construction of a polynomial chaos expansion (PCE) of the potential $V$ in the Gaussian Sobolev space [18], and introduction of a convergence analysis for the truncated PCE representation of $V$. Transformation of the truncated PCE into a polynomial in $y^m$, in which $m$ is a multi-index of dimension $v$. This is the final representation of $V$ that is perfectly adapted for developing a quantum algorithm based on a second quantization.
(iv) Development of a quantum computing formulation based on the following three steps: (a) Defining a finite basis representation, for which the eigenvalue problem can be rewritten in a second quantized form. (b) Developing the second quantization, introducing the creation and annihilation operators in order to rewrite the eigenvalue problem in the Fock space, allowing us to obtain explicit formulas for the Laplacian operator $\nabla^2$ and for the multivariate monomials $y^m$. Deducing a representation of the Fokker-Planck operator as a linear combination with known real coefficients of products of powers of the creation and annihilation operators. (c) Mapping onto a system of qubits by expressing the products of creation and annihilation operators with the Pauli matrices.

(v) Design of the quantum circuits and measurements in order to construct and measure specific states on the quantum computer. (a) The quantum circuits are constructed with gates. To construct the eigenstates, a set of universal quantum gates is used from which any unitary transformation can be reconstructed. The successive application of these gates on the set of qubits forms the quantum circuits. (b) The extraction of information from the quantum states is performed via quantum measurements of certain physical observables, which are probabilistic.

(vi) Finally, we perform the eigenstate construction and overlap extraction.

## 1.1. Convention for the variables, vectors, and matrices

*x*, *η*: lower-case Latin or Greek letters are deterministic real variables.

*x*, *η*: boldface lower-case Latin or Greek letters are deterministic vectors.

*X*: upper-case Latin letters are real-valued random variables.

*X*: boldface upper-case Latin letters are vector-valued random variables.

*[x]*: lower-case Latin letters between brackets are deterministic matrices.

*[X]*: boldface upper-case letters between brackets are matrix-valued random variables.

## 1.2. Convention used for random variables

In this paper, for any finite integer $m \ge 1$, the Euclidean space $\mathbb{R}^m$ is equipped with the $\sigma$-algebra $\mathcal{B}_{\mathbb{R}^m}$. If $Y$ is a $\mathbb{R}^m$-valued random variable defined on the probability space $(\Theta, \mathcal{T}, \mathcal{P})$, $Y$ is a mapping $\theta \mapsto Y(\theta)$ from $\Theta$ into $\mathbb{R}^m$, measurable from $(\Theta, \mathcal{T})$ into $(\mathbb{R}^m, \mathcal{B}_{\mathbb{R}^m})$, and $Y(\theta)$ is a realization (sample) of $Y$ for $\theta \in \Theta$. The probability distribution of $Y$ is the probability measure $P_Y(dy)$ on the measurable set $(\mathbb{R}^m, \mathcal{B}_{\mathbb{R}^m})$ (we will simply say on $\mathbb{R}^m$). The Lebesgue measure on $\mathbb{R}^m$ is noted $dy$ and when $P_Y(dy)$ is written as $p_Y(y)dy$, $p_Y$ is the probability density function (pdf) on $\mathbb{R}^m$ of $P_Y(dy)$ with respect to $dy$. Finally, $E$ denotes the mathematical expectation operator that is such that $E\{Y\} = \int_{\mathbb{R}^m} y P_Y(dy)$.

# 2. Probabilistic learning on manifolds with a given probability measure

## 2.1. Preamble related to the context of the probabilistic learning on manifolds

Probabilistic Learning on Manifolds (PLoM) is a tool in computational statistics, introduced in 2016 [19], which can be viewed as a tool for scientific machine learning. The PLoM approach has specifically been developed for small training-dataset cases [19, 20]. The method avoids the scattering of learned realizations associated with the probability distribution to preserve its concentration in the neighborhood of the random manifold defined by the parameterized computational model. Several extensions have been proposed to account for implicit constraints induced by physics, computational models, and measurements [21, 22, 23, 24], to reduce the stochastic dimension using a statistical partition approach [25], and to update the prior probability distribution with a target dataset, whose points are, for instance, experimental realizations of the system observations [26]. Consequently, PLoM, constrained by a stochastic computational model and statistical moments or samples/realizations, allows for performing probabilistic learning inference and constructing predictive statistical surrogate models for large parameterized stochastic computational models.

PLoM approach relies on projecting, along the data axis, a matrix-valued Itô equation linked to a stochastic dissipative Hamiltonian system, serving as the MCMC generator for the probability measure estimated via Gaussian Kernel Density Estimation (GKDE) on the training dataset. The projection basis is the diffusion maps (DMAPS) basis, associated with a time-independent isotropic kernel introduced in [27, 28], referred to as the reduced-order DMAPS basis in the PLoM context. Since 2016, all published PLoM applications have used the isotropic kernel,
yielding quality results even for heterogeneous data and complex systems in various dimensions. Improving statistical surrogates for stochastic manifolds with conditional statistics via PLoM using a projection basis from a transient anisotropic (time-dependent) kernel has recently been proposed [1].

The first algorithmic step of PLoM [19, 20] involves performing a principal component analysis (PCA) on the non-Gaussian random vector $X$, which is generally in high dimension. This vector $X$ results from concatenating the random vector $Q$, representing the quantities of interest (QoI), with the random vector $W$, representing the random control variables of the considered stochastic system. Additionally, this system depends on random latent/uncontrolled variables. The vector-valued random QoI, $Q = F(W)$, is connected to $W$ by an unknown random nonlinear mapping. Denoting by $C_w$ the support of the probability measure of $W$, the random graph $\{(w, F(w)) | w \in C_w\}$ forms a stochastic manifold due to the latent random variables in the system implying the randomness of mapping $F$. The PCA is performed on the training dataset of $X$. The new coordinates resulting from the PCA form a random vector $H$ with a reduced dimension. The training dataset of $H$ is obtained through the PCA projection of the training dataset of $X$. By construction, $H$ is a non-Gaussian normalized random vector, which is centered and has a covariance matrix equal to the identity matrix. Under these conditions, the available information is only the training dataset of $H$. Using this dataset, the probability density function $p_H$ of $H$ is constructed by using the modification [29] of the GKDE method [30]. Thus, the only available information is represented by the probability density function $p_H$ of $H$. Given the normalization of $H$, the challenges in constructing the learned dataset of $H$, with or without constraints, involve preserving the concentration of the learned probability measure for $H$ (induced by the random manifold) and ensuring the “good quality” of the learned joint probability measure of the components of $H$, not just the marginal probability measure of each component. The quality of the statistical surrogate model of $Q$ given $W$, based on conditional statistics, is directly linked to the learned joint probability measure of the components of $H$. Here, to simplify the presentation, we will directly begin on $H$ rather than starting from $X = (Q, W)$.

The method presented in [1] starts from the FKP equation, whose stationary solution is $p_H$. The spectral problem of this FKP operator is introduced, and its eigenfunctions would form the ideal basis for performing the data projection necessary to preserve the concentration of the learned probability measure. Given the high dimensionality and small training datasets, a change of scale (smoothing) is introduced, similar to DMAPS, to calculate the eigenvectors of a matrix that approximates the FKP operator. This approximation is essential because computing the eigenfunctions of the FKP operator in high dimension is beyond the capabilities of classical computers. In this paper, we revisit this spectral problem and construct a formulation adapted to quantum computing, potentially enabling the direct calculation of eigenfunctions in high dimensions.

## 2.2. Defining the probability measure of random vector $H$

Let $\mathcal{D}_{\text{train}}(\eta) = \{\eta^1, \dots, \eta^{n_d}\}$ be the training dataset of $n_d > 1$ independent realizations $\eta^j \in \mathbb{R}^v$, with $v \ge 1$, of a second-order $\mathbb{R}^v$-valued random variable $H$, defined on a probability space $(\Theta, \mathcal{T}, \mathcal{P})$ and which results of the PCA of random vector $X$ (see Section 2.1. Let $\underline{\eta}_d \in \mathbb{R}^v$ and $[C_d] \in \mathbb{M}_{n_d}$ be the associated empirical estimates of the mean value and the covariance matrix constructed with the points of $\mathcal{D}_{\text{train}}(\eta)$,

$$ \underline{\eta}_d = \frac{1}{n_d} \sum_{j=1}^{n_d} \eta^j, \quad [C_d] = \frac{1}{n_d - 1} \sum_{j=1}^{n_d} (\eta^j - \underline{\eta}_d) \otimes (\eta^j - \underline{\eta}_d). \quad (2.1) $$

The properties of the PCA imply that

$$ \underline{\eta}_d = \mathbf{0}_v, \quad [C_d] = [I_v]. \quad (2.2) $$

Let $y \mapsto p_H(y)$ be the probability density function on $\mathbb{R}^v$, with respect do the Lebesgue measure $dy$, defined by

$$ p_H(y) = \frac{1}{n_d} \sum_{j=1}^{n_d} \frac{1}{(\sqrt{2\pi} \hat{s})^v} \exp\left(-\frac{1}{2\hat{s}^2} \|y - \underline{\eta}^j\|^2\right), \quad \forall y \in \mathbb{R}^v, \quad (2.3) $$

where $\hat{s}$ and $s$ are defined by

$$ s = \left( \frac{4}{n_d(2+v)} \right)^{1/(v+4)}, \quad \hat{s} = \frac{s}{\sqrt{s^2 + (n_d - 1)/n_d}}. \quad (2.4) $$
Eqs. (2.3) and (2.4) correspond to the Gaussian kernel-density estimation (KDE) constructed using the $n_d$ independent realizations of $D_{train}(\eta)$ involving the modification [29] of the usual formulation [30, 31, 32, 33], in which $s$ is the Silverman bandwidth. Let $\mathbf{H}$ be the second-order $\mathbb{R}^v$-valued random variable, defined on a probability space $(\Theta, \mathcal{T}, \mathcal{P})$, whose probability measure $P_H(dy) = p_H(\mathbf{y}) dy$ on $\mathbb{R}^v$ is defined by the probability density function $p_H$ given by Eq. (2.3). It can be seen that, for any fixed $n_d > 1$, we have exactly,

$$E\{\mathbf{H}\} = \int_{\mathbb{R}^v} \mathbf{y} P_H(dy) = \frac{1}{2\hat{s}^2} \underline{\eta}_d = \mathbf{0}_v, \qquad (2.5)$$

$$E\{\mathbf{H} \otimes \mathbf{H}\} = \int_{\mathbb{R}^v} \mathbf{y} \otimes \mathbf{y} P_H(dy) = \hat{s}^2 [I_v] + \frac{\hat{s}^2}{s^2} \frac{(n_d - 1)}{n_d} [C_d] = [I_v]. \qquad (2.6)$$

Eqs. (2.5) and (2.6) show that $\mathbf{H}$ is a normalized $\mathbb{R}^v$-valued random variable. The probability density function $p_H$ defined by Eq. (2.3) is rewritten, for all $\mathbf{y}$ in $\mathbb{R}^v$, as

$$p_H(\mathbf{y}) = c_v \xi(\mathbf{y}) \quad , \quad \xi(\mathbf{y}) = e^{-\Phi(\mathbf{y})}, \qquad (2.7)$$

in which $c_v = (\sqrt{2\pi} \hat{s})^{-v}$ and where $\Phi(\mathbf{y}) = -\log(\xi(\mathbf{y}))$ is such that

$$\Phi(\mathbf{y}) = -\log\left(\frac{1}{n_d} \sum_{j=1}^{n_d} \exp\left(-\frac{1}{2\hat{s}^2}\|\mathbf{y} - \frac{\hat{s}}{s}\mathbf{\eta}'\|^2\right)\right) \quad , \quad \forall \mathbf{y} \in \mathbb{R}^v. \qquad (2.8)$$

### 3. eigenvalue problem of the FKP operator and a representation adapted to quantum computing algorithm

This section summarizes essential results concerning the Fokker-Planck (FKP) operator and its associated eigenvalue problem, which are formally presented. Then, the FKP operator is transformed to obtain a Schrödinger-type formulation in which $\mathcal{V}(\mathbf{y})$ will be the potential on $\mathbb{R}^v$. Finally, we propose a methodology to construct a representation of the potential $\mathcal{V}(\mathbf{y})$ adapted to quantum computing algorithms, which is presented in Section 6. This representation will be based on a polynomial chaos expansion (PCE) in Gaussian space of $\mathcal{V}(\mathbf{y})$, which will be explicitly constructed in Section 4.

#### 3.1. Itô stochastic differential equation related to probability density function $p_H$

We introduce an Itô stochastic differential equation (ISDE) on $\mathbb{R}^v$, with initial condition, for which $P_H(dy) = p_H(\mathbf{y}) dy$ is the invariant measure. A classical candidate to such an ISDE is written as

$$d\mathbf{Y}(t) = \mathbf{b}(\mathbf{Y}(t)) dt + d\mathbf{W}(t) \quad , \quad t > 0, \qquad (3.1)$$

$$\mathbf{Y}(0) = \mathbf{x} \in \mathbb{R}^v, \quad \text{a.s.}, \qquad (3.2)$$

where the drift vector is the function $\mathbf{y} \mapsto \mathbf{b}(\mathbf{y})$ from $\mathbb{R}^v$ into $\mathbb{R}^v$ defined by

$$\mathbf{b}(\mathbf{y}) = -\frac{1}{2} \nabla \Phi(\mathbf{y}) \quad , \quad \forall \mathbf{y} \in \mathbb{R}^v. \qquad (3.3)$$

In Eq. (3.1), $\mathbf{W}(t) = (W_1(t), \dots, W_v(t))$ is the normalized Wiener stochastic process [34] on $\mathbb{R}^+ = [0, +\infty[$, with values in $\mathbb{R}^v$, which is a stochastic process with independent increments, such that $\mathbf{W}(0) = \mathbf{0}_v$, a.s., and for $0 \le \tau < t < +\infty$, the increment $\Delta\mathbf{W}_{\tau t} = \mathbf{W}(t) - \mathbf{W}(\tau)$ is a Gaussian $\mathbb{R}^v$-valued second-order random variable, centered and with a covariance matrix that is written as

$$[C_{\Delta\mathbf{W}_{\tau t}}] = E\{\Delta\mathbf{W}_{\tau t} \otimes \Delta\mathbf{W}_{\tau t}\} = (t - \tau)[I_v]. \qquad (3.4)$$

We will see that $\{\mathbf{Y}(t), t \in \mathbb{R}^+\}$ is a homogeneous diffusion stochastic process, which is asymptotically stationary for $t \to +\infty$. Assuming that the transition probability measure of $\mathbf{Y}(t)$ given $\mathbf{Y}(0) = \mathbf{x}$ admits a density with respect to $dy$, such that, for all $t > 0$, for all $\mathbf{x}$ and $\mathbf{y}$ in $\mathbb{R}^v$, and for any Borelian $\mathcal{B}$ in $\mathbb{R}^v$, we have

$$\mathcal{P}(\{\mathbf{Y}(t) \in \mathcal{B} | \mathbf{Y}(0) = \mathbf{x}\} = \int_{\mathcal{B}} \rho(\mathbf{y}, t | \mathbf{x}, 0) dy, \qquad (3.5)$$

$$\lim_{t \to 0^+} \rho(\mathbf{y}, t | \mathbf{x}, 0) dy = \delta_0(\mathbf{y} - \mathbf{x}), \qquad (3.6)$$

$$\int_{\mathbb{R}^v} \rho(\mathbf{y}, t | \mathbf{x}, 0) dy = 1. \qquad (3.7)$$
### 3.2. FKP equation associated with the ISDE, and its steady-state solution

For all $\mathbf{x}$ in $\mathbb{R}^v$, the transition probability density function $(y, t) \mapsto \rho(y, t | \mathbf{x}, 0)$ from $\mathbb{R}^v \times \mathbb{R}^+$ into $\mathbb{R}^+$ verifies the following Fokker-Planck (FKP) equation (see for instance [35, 36, 37]),

$$ \frac{\partial \rho}{\partial t} + L_{\text{FKP}}(\rho) = 0, \quad t > 0, \qquad (3.8) $$

with the initial condition for $t = 0$ defined by Eq. (3.6). The Fokker-Planck operator $L_{\text{FKP}}$ can be written, after a small algebraic manipulation and for any sufficiently differentiable function $y \mapsto v(y)$ from $\mathbb{R}^v$ into $\mathbb{R}$, as

$$ \{L_{\text{FKP}}(v)(y)\} = -\frac{1}{2}\text{div}\left\{p_H(y)\nabla\left(\frac{v(y)}{p_H(y)}\right)\right\}. \qquad (3.9) $$

The detailed balance (the probability current vanishes) is satisfied and the steady-state solution of Eq. (3.8) is the pdf $p_H$ defined by Eq. (2.3) [37, 12, 11]. We then have, for $v = p_H$,

$$ L_{\text{FKP}}(p_H) = 0. \qquad (3.10) $$

The invariant measure $P_H(dy) = p_H(y)dy$ is such that, for all $y$ in $\mathbb{R}^v$ and for all $t \ge 0$,

$$ p_H(y) = \int_{\mathbb{R}^v} \rho(y, t | \mathbf{x}, 0) p_H(\mathbf{x}) d\mathbf{x}. \qquad (3.11) $$

The ISDE defined by Eqs. (3.1) and (3.2) admits an asymptotic ($t \to +\infty$) stationary solution whose marginal probability density function of order one is $p_H$. Consequently, for all $\mathbf{x}$ and $y$ in $\mathbb{R}^v$, we have

$$ \lim_{t \to +\infty} \rho(y, t | \mathbf{x}, 0) = p_H(y). \qquad (3.12) $$

### 3.3. Formal formulation of the eigenvalue problem of the FKP operator

The eigenvalue problem, posed in an adapted functional space, is written as

$$ L_{\text{FKP}}(v) = \lambda v, \qquad (3.13) $$

for which the current must vanish at infinity, yielding the condition,

$$ \lim_{\|y\| \to +\infty} p_H(y) \| \nabla(p_H(y)^{-1}v(y)) \| = 0. \qquad (3.14) $$

Continuing the development within a formal framework, such as that used in [11], we introduce the change of function,

$$ v(y) = p_H(y)^{1/2} q(y), \quad y \in \mathbb{R}^v, \quad q: \mathbb{R}^v \to \mathbb{R}. \qquad (3.15) $$

Let $\hat{L}_{\text{FKP}}$ be the linear operator defined, for $q$ belonging to an admissible set of functions,

$$ \{\hat{L}_{\text{FKP}}(q)\}(y) = p_H(y)^{-1/2} L_{\text{FKP}}(p_H(y)^{1/2} q(y)), \quad y \in \mathbb{R}^v. \qquad (3.16) $$

Therefore, the eigenvalue problem defined by Eqs. (3.13) and (3.14) can be rewritten in $q$ as

$$ \hat{L}_{\text{FKP}}(q) = \lambda q, \qquad (3.17) $$

with the condition at infinity,

$$ \lim_{\|y\| \to +\infty} p_H(y) \| \nabla(p_H(y)^{-1/2} q(y)) \| = 0. \qquad (3.18) $$
### 3.4. Schrödinger-type formulation of the FKP operator $\hat{L}_{FKP}$

Using Eq. (2.7), which shows that $p_H(\mathbf{y})^{-1}\nabla p_H(\mathbf{y}) = -\nabla\Phi$, and using Eqs. (3.9) and (3.16), it can be seen that

$$ \{\hat{L}_{FKP}(q)\}(\mathbf{y}) = \mathcal{V}(\mathbf{y}) q(\mathbf{y}) - \frac{1}{2} \nabla^2 q(\mathbf{y}) , \quad \mathbf{y} \in \mathbb{R}^v, \qquad (3.19) $$

in which $\nabla^2$ is the Laplacian operator in $\mathbb{R}^v$ and where $\mathbf{y} \mapsto \mathcal{V}(\mathbf{y})$ is the function from $\mathbb{R}^v$ into $\mathbb{R}$, which is defined, for all $\mathbf{y}$ in $\mathbb{R}^v$, as

$$ \mathcal{V}(\mathbf{y}) = \frac{1}{8} \| \nabla \Phi(\mathbf{y}) \|^2 - \frac{1}{4} \nabla^2 \Phi(\mathbf{y}). \qquad (3.20) $$

### 3.5. Properties of operator $\hat{L}_{FKP}$ and hypothesis on its spectrum

Let $\delta q$ be a function from $\mathbb{R}^v$ into $\mathbb{R}$, belonging to the admissible set that allows the evaluation of the bracket

$$ \langle \hat{L}_{FKP}(q), \delta q \rangle = \int_{\mathbb{R}^v} \{\hat{L}_{FKP}(q)\}(\mathbf{y}) \delta q(\mathbf{y}) \, d\mathbf{y}. $$

Removing $\mathbf{y}$ and using Eq. (3.16) with Eq. (3.9) yields

$$ \langle \hat{L}_{FKP}(q), \delta q \rangle = -\frac{1}{2} \int_{\mathbb{R}^v} (p_H^{-1/2} \delta q) \operatorname{div}\{p_H \nabla (p_H^{-1/2} q)\} d\mathbf{y}. \qquad (3.21) $$

Using the condition at infinity, defined by Eq. (3.18), Eq. (3.21) can be rewritten as,

$$ \langle \hat{L}_{FKP}(q), \delta q \rangle = \frac{1}{2} \int_{\mathbb{R}^v} p_H \langle \nabla(p_H^{-1/2} q), \nabla(p_H^{-1/2} \delta q) \rangle_{\mathbb{R}^v} d\mathbf{y}. \qquad (3.22) $$

(a) Eq. (3.22) shows that $\hat{L}_{FKP}$ is a symmetric and positive operator.

(b) Eqs. (3.10) and (3.16) show that

$$ \hat{L}_{FKP}(q_0) = 0 \quad \text{for} \quad q_0 = p_H^{1/2}. \qquad (3.23) $$

Since $p_H(\mathbf{y}) d\mathbf{y}$ is a bounded positive measure (probability measure), the right-hand side of Eq. (3.22) shows that the null space of $\hat{L}_{FKP}$ is of dimension 1 and is constituted of the function $q_0 = p_H^{1/2}$. For $q = \delta q \neq q_0$, and $\|q_0\| \neq 0$, we have $\langle \hat{L}_{FKP}(q), q \rangle > 0$. Therefore, $\hat{L}_{FKP}$ is a positive operator (in the quotient space by the null space).

It is assumed that $p_H$ defined by Eq. (2.3), which is constructed with the $n_d$ points $\{\eta^j, j = 1, \dots, n_d\}$ of the training dataset, is such that the spectrum of $\hat{L}_{FKP}$ is countable. Due to (a) and (b), we then deduce that the eigenvalues of $\hat{L}_{FKP}$ (defined by Eqs. (3.17) and (3.18)) are positive except one that is zero. We will also assume that the multiplicity of each eigenvalue is finite and that the eigenfunctions are continuous on $\mathbb{R}^v$.

### 3.6. Eigenvalue problem for operator $\hat{L}_{FKP}$ and objective

Under the hypothesis introduced in Section 3.5, the eigenvalue problem

$$ \hat{L}_{FKP}(q_m) = \lambda_m q_m , \quad m \in \mathbb{N}, \qquad (3.24) $$

for operator $\hat{L}_{FKP}$, with the condition defined by Eq. (3.18), is such that

$$ 0 = \lambda_0 < \lambda_1 \leq \lambda_2 \leq \dots, \qquad (3.25) $$

the multiplicity of each eigenvalue being finite and the associated eigenfunctions $\{q_m, m \in \mathbb{N}\}$ belong to $L^2(\mathbb{R}^v)$. We will admit that the family $\{q_m, m \in \mathbb{N}\}$ is a Hilbert basis of $L^2(\mathbb{R}^v)$. We then have

$$ \langle q_m, q_n \rangle_{L^2} = \int_{\mathbb{R}^v} q_m(\mathbf{y}) q_n(\mathbf{y}) d\mathbf{y} = \delta_{mn}. \qquad (3.26) $$

The eigenfunction $q_0$ associated with $\lambda_0 = 0$, is such that (see Eq. (3.23)),

$$ q_0 = p_H^{1/2} , \quad \| q_0 \|_{L^2} = 1, \qquad (3.27) $$
and we have

$$ \sum_{m \in \mathbb{N}} q_m(\mathbf{y}) q_m(\mathbf{x}) d\mathbf{y} = \delta_0(\mathbf{y} - \mathbf{x}). \qquad (3.28) $$

From Eqs. (3.26) and (3.27), it can be deduced that

$$ \forall m \ge 1, \quad \int_{\mathbb{R}^v} p_H(\mathbf{y})^{1/2} q_m(\mathbf{y}) d\mathbf{y} = 0. \qquad (3.29) $$

Here, the objective is to compute with quantum computing algorithm, the values of the eigenfunctions $q_m$, for $m = 1, \dots, M$, associated with the eigenvalues $0 = \lambda_0 < \lambda_1 < \dots < \lambda_M$ with $M < n_d$ at the sampling points $\eta^1, \dots, \eta^{n_d}$, that is to say, for fixed $m, q_m(\eta^1), \dots, q_m(\eta^{n_d})$.

### 3.7. Proposed methodology to construct a representation of the potential $V$ adapted to quantum computing algorithm

As we explained in Section 1, within the formulation we propose for the quantum computing algorithm, it is necessary to construct a polynomial representation in $\mathbf{y}$ of the potential $V(\mathbf{y})$. Clearly, the expression of $V(\mathbf{y})$, as defined by Eq.(3.20), where $\Phi(\mathbf{y})$ is defined by Eq.(2.8), is not a polynomial. Therefore, we need to construct an approximation of $V(\mathbf{y})$ in the space of multivariate polynomials on $\mathbb{R}^v$, for instance, using monomials on $\mathbb{R}^v$. Such a construction is not straightforward for the following reasons:

(i) The problem is set in high dimensions, meaning that $v$ is large.

(ii) The sequence of polynomials chosen from the selected subset of all the polynomials on $\mathbb{R}^v$ must be complete in the adapted vector space, and the convergence must be fast with respect to the number of polynomials and their maximum degrees.

(iii) Eq. (3.20) shows that potential $V(\mathbf{y})$ is a linear combination of a positive-valued term $\|\nabla\Phi(\mathbf{y})\|^2$ and a real-valued term $\nabla^2\Phi(\mathbf{y})$, where $\Phi(\mathbf{y})$ is defined by Eq. (2.8). Although $V(\mathbf{y})$ is not, a priori, positive for all $\mathbf{y}$ in $\mathbb{R}^v$, there is an underlying algebraic structure in the expression of $V(\mathbf{y})$, and these induced algebraic properties must be preserved.

There is an important property of the problem that can and will be used to select a well-adapted family of polynomials. The probability density function $p_H$ corresponds to the normalized random variable $H$, which is therefore centered and whose covariance matrix is equal to the identity matrix. Thus, if the training dataset consisted of $n_d$ independent realizations of a normalized Gaussian random vector, then as $n_d \to +\infty$, $p_H$ would be Gaussian, and $p_H$ would coincide with the canonical normal density on $\mathbb{R}^v$, denoted $p_Y$. In this case, we would have $\Phi(\mathbf{y}) = \|\mathbf{y}\|^2/2$, leading to $\|\nabla\Phi(\mathbf{y})\|^2 = \|\mathbf{y}\|^2$ and $\nabla^2\Phi(\mathbf{y}) = v$, which shows that $V(\mathbf{y})$ would be written as $\|\mathbf{y}\|^2/8 - v/4$. Consequently, the choice of normalized Hermite polynomials, which constitute a Hilbert basis in the Hilbert space of all square integrable functions on $\mathbb{R}^v$ for the measure $p_Y(\mathbf{y})d\mathbf{y}$, is certainly a good choice with respect to (i) and (ii). Therefore, the choice will be the polynomial chaos expansion (PCE) in the Gaussian space. Nevertheless, to satisfy (iii), it is preferable to construct such a PCE for $f(\mathbf{y}) = \log(p_H(\mathbf{y}))$ and then, to deduce of it the PCE of $\|\nabla f(\mathbf{y})\|^2$ and of $\nabla^2 f(\mathbf{y})$, which allows for the PCE of $V(\mathbf{y})$ to be obtained, knowing that $\nabla\Phi(\mathbf{y}) = -\nabla f(\mathbf{y})$ and $\nabla^2\Phi(\mathbf{y}) = \nabla^2 f(\mathbf{y})$ and so

$$ V(\mathbf{y}) = \frac{1}{8} \| \nabla f(\mathbf{y}) \|^2 + \frac{1}{4} \nabla^2 f(\mathbf{y}). \qquad (3.30) $$

We will then use the Gaussian Sobolev space [18]. Finally, the following comments should be noted:

(a) Performing a PCE of $p_H$, instead of $f(\mathbf{y}) = \log(p_H(\mathbf{y}))$, would not be convenient because it is difficult to preserve the positivity of the PCE of $p_H$ for all $\mathbf{y}$ in $\mathbb{R}^v$ with a small number of polynomials of low degrees.

(b) Performing a direct construction of the PCE of $\nabla\Phi(\mathbf{y})$, and even more so for $\nabla^2\Phi(\mathbf{y})$, would also be difficult to satisfy (ii) because $\nabla\Phi(\mathbf{y}) = -p_H(\mathbf{y})^{-1}\nabla p_H(\mathbf{y})$, showing that $p_H(\mathbf{y})$ appears to the denominator.

### 3.8. Use of the eigenfunctions of the FKP operator in PLoM

The PLoM approach was originally developed using a reduced-order diffusion-maps basis (RODB) {$g_{DM}^1, \dots, g_{DM}^{m_{opt}}$} in $\mathbb{R}^{n_d}$, represented by the matrix $[g_{DM}] \in \mathbb{M}_{n_d, m_{opt}}$ with $m_{opt} < n_d$. Such an RODB is constructed using the isotropic kernel built from the training dataset $\mathcal{D}_{train}(\eta)$, which consists of the $n_d$ points $\eta^1, \dots, \eta^{n_d}$ in $\mathbb{R}^d$ (see [19, 20]). The optimal dimension $m_{opt}$ is determined as explained in [25] or possibly by using the algorithm presented in [38].
If the data in the training dataset are heterogeneous, we proposed in [1] to replace the RODB with a reduced-order basis {$g_{TB}^1, \dots, g_{TB}^{m_{opt}}$} (ROTB) in $\mathbb{R}^{n_d}$, represented by the matrix $[g_{TB}] \in \mathbb{M}_{n_d, m_{opt}}$. This ROTB is constructed using a transient anisotropic kernel based on the transition probability density function, which is the transient solution of the Fokker-Plank equation defined by Eq. (3.8) and estimated by solving the Itô stochastic differential equation defined by Eqs. (3.1) and (3.2). It should be noted that in this construction, we introduced a scaling change in the transient anisotropic kernel to connect the ROTB to the RODB in the vicinity of the initial time $t = 0$. This approach was proposed because the use of eigenfunctions from the eigenvalue problem of the $L_{FKP}$ operator, defined by Eq. (3.13), could not be explored due to the fact that, in high dimensions, such an eigenvalue problem cannot be solved using classical computers.

This is the reason why, in the present paper, we explore a methodology for calculating the eigenfunctions associated with the smallest eigenvalues of the $L_{FKP}$ operator using quantum computing, with the aim of developing an approach that could allow, in the near future, replacing the ROTB with the reduced-order basis {$g_{FKP}^1, \dots, g_{FKP}^{m_{opt}}$} in $\mathbb{R}^{n_d}$, represented by the matrix $[g_{FKP}] \in \mathbb{M}_{n_d, m_{opt}}$. For $m = 1, \dots, m_{opt}$, the vector $g_{FKP}^m$ in $\mathbb{R}^{n_d}$ is constructed as the value of the eigenfunction $y \mapsto v_m(y)$ associated with the eigenvalues $\lambda_m$ at the point $\eta^j$ in $\mathbb{R}^d$. From Eqs. (3.13), (3.15), and (3.24) with (3.25), we obtain the component $g_{FKP,j}^m$ of vector $g_{FKP}^m \in \mathbb{R}^{n_d}$, where $j \in \{1, \dots, n_d\}$,

$$g_{FKP,j}^m = v_m(\eta^j) = p_H(\eta^j)^{1/2} q_m(\eta^j). \qquad (3.31)$$

## 4. Polynomial chaos expansion of potential $V$ in the Gaussian space for quantum computing algorithm

In this section, we introduce the functional spaces necessary to construct the Gaussian polynomial chaos expansion. We then present and prove such an expansion of $V$ based on the strategy proposed in Section 3.7. Finally we introduce the convergence criteria of the truncated Chaos representation of $V$.

### 4.1. Properties of the function $f = \log(p_H)$ and its two first derivatives

We begin with Definition 1 of functional spaces and then proceed with Proposition 1, which provides the properties of the function $f = \log(p_H)$ and its two first derivatives.

**Definition 1 (Hilbert spaces $L^2(\Theta, \mathbb{R})$ and $H = L^2(\mathbb{R}^v; p_Y(y) dy)$).** (i) Let $Y = (Y_1, \dots, Y_v)$ be a $\mathbb{R}^v$-valued normalized Gaussian random variable defines on a probability space $(\Theta, \mathcal{T}, \mathcal{P})$, whose probability measure $p_Y(y) dy$ is defined by the canonical normal density $y \mapsto p_Y(y)$ on $\mathbb{R}^v$,

$$p_Y(y) = (2\pi)^{-v/2} \exp(-\|y\|^2/2). \qquad (4.1)$$

Let $L^2(\Theta, \mathbb{R})$ be the Hilbert space of all second-order real-valued random variables defined on the probability space $(\Theta, \mathcal{T}, \mathcal{P})$ (equivalence classes of almost-surely equal random variables), equipped with the inner product and the associated norm,

$$\langle U, V \rangle_{\Theta} = E\{UV\} = \int_{\Theta} U(\theta) V(\theta) d\mathcal{P}(\theta) \quad , \quad \|U\|_{\Theta} = (E\{U^2\})^{1/2}. \qquad (4.2)$$

(ii) Let $H = L^2(\mathbb{R}^v; p_Y(y) dy)$ be the Hilbert space of all the real-valued functions on $\mathbb{R}^v$, which are square integrable on $\mathbb{R}^v$ with respect to the measure $p_Y(y) dy$, equipped with the inner product and the associated norm,

$$\langle u, v \rangle_{H} = \int_{\mathbb{R}^v} u(y) v(y) p_Y(y) dy \quad , \quad \|u\|_{H}^{2} = \langle u, u \rangle_{H}. \qquad (4.3)$$

For all continuous functions $u$ and $v$ in $H$, $U = u(Y)$ and $V = v(Y)$ are real-valued random variables in $L^2(\Theta, \mathbb{R})$,

$$\langle u, v \rangle_{H} = E\{u(Y) v(Y)\} = E\{UV\} = \langle U, V \rangle_{\Theta} \quad , \quad \|u\|_{H}^{2} = E\{u(Y)^2\} = E\{U^2\} = \|U\|_{\Theta}^{2}. \qquad (4.4)$$

**Proposition 1 (Properties of function $f = \log(p_H)$).** Let $v \ge 1$ and $n_d \ge 1$ be fixed integers, and let $\mathcal{D}_{train}(\eta)$ be the training dataset of points $\eta^1, \dots, \eta^{n_d}$ in $\mathbb{R}^v$, defined in Section 2.2. Let $p_H$ and $p_Y$ be defined by Eqs. (2.7) and (4.1). (i) The real-valued function $f$ is continuous on $\mathbb{R}^v$ and belongs to $H$,

$$E\{f(Y)^2\} = \|f\|_{H}^{2} < +\infty. \qquad (4.5)$$
(ii) Let $\mathbf{g} = (g_1, \dots, g_\nu)$ be the gradient function defined by $\mathbf{g} = \nabla f = (\partial_1 f, \dots, \partial_\nu f)$ where $\partial_j = \partial/\partial y_j$ is a $\mathbb{R}^\nu$-valued continuous function on $\mathbb{R}^\nu$. For all $j = 1, \dots, \nu$, the function $g_j$ belongs to $\mathbb{H}$ and the function $\hat{g} = \|\mathbf{g}\|$ belongs to $\mathbb{H}$,

$$ E\{g_j(\mathbf{Y})^2\} = \|g_j\|_{\mathbb{H}}^2 < +\infty, \quad E\{\hat{g}(\mathbf{Y})^2\} = \|\hat{g}\|_{\mathbb{H}}^2 < +\infty, \quad E\{\hat{g}(\mathbf{Y})^4\} < +\infty \quad (4.6) $$

(ii) The Laplacian function $h = \nabla^2 f = \partial_1^2 f + \dots + \partial_\nu^2 f$ is a real-valued continuous function on $\mathbb{R}^\nu$ and the function $h$ belongs to $\mathbb{H}$,

$$ E\{|h(\mathbf{Y})|^2\} = \|h\|_{\mathbb{H}}^2 < +\infty. \quad (4.7) $$

PROOF. (Proposition 1). The continuity of functions $f$, $g$, and $h$ is simple to be prove.

(i) Let $p_H(\mathbf{y}) = c_\nu \xi(\mathbf{y})$ be defined by Eq. (2.7) where $\xi(\mathbf{y}) = \frac{1}{n_d} \sum_{j=1}^{n_d} \exp(-\frac{1}{2s^2} \|\mathbf{y} - \hat{s}_s \boldsymbol{\eta}^j\|^2)$. For all $\mathbf{y}$ in $\mathbb{R}^\nu$, we have $0 < \xi(\mathbf{y}) \le 1$ and $(\log p_H(\mathbf{y}))^2 = (\log c_\nu + \log \xi(\mathbf{y}))^2 \le 2(\log c_\nu)^2 + 2(\log \xi(\mathbf{y}))^2$. Consequently, $E\{f(\mathbf{Y})^2\} \le 2(\log c_\nu)^2 + 2E\{(\log \xi(\mathbf{Y}))^2\}$. Since $z \mapsto (\log z)^2$ is a convex function on $[0, 1]$, and since $\xi(\mathbf{y}) \in [0, 1]$ for all $\mathbf{y}$ in $\mathbb{R}^\nu$, using the Jensen inequality yields, $(\log \xi(\mathbf{y}))^2 \le \frac{1}{n_d} \sum_{j=1}^{n_d} \{\log(\exp(-\frac{1}{2s^2} \|\mathbf{y} - \hat{s}_s \boldsymbol{\eta}^j\|^2))\}^2 = \frac{1}{n_d} \sum_{j=1}^{n_d} \frac{1}{4s^4} \|\mathbf{y} - \hat{s}_s \boldsymbol{\eta}^j\|^4$. Therefore, we have $E\{(\log \xi(\mathbf{Y}))^2\} < +\infty$, which proves Eq. (4.5).

(ii) We have $g(\mathbf{y}) = p_H(\mathbf{y})^{-1} \nabla p_H(\mathbf{y}) = \xi(\mathbf{y})^{-1} \nabla \xi(\mathbf{y})$ and since $\nabla \xi(\mathbf{y}) = -\frac{1}{n_d} \frac{1}{s^2} \sum_{j=1}^{n_d} (\mathbf{y} - \hat{s}_s \boldsymbol{\eta}^j) \exp(-\frac{1}{2s^2} \|\mathbf{y} - \hat{s}_s \boldsymbol{\eta}^j\|^2)$, we deduce that $\|g(\mathbf{y})\| \le \frac{1}{s^2} \|\mathbf{y}\| + \frac{1}{n_d} \frac{1}{s^2 \xi(\mathbf{y})} \sum_{j=1}^{n_d} \hat{s}_s \|\boldsymbol{\eta}^j\| \exp(-\frac{1}{2s^2} \|\mathbf{y} - \hat{s}_s \boldsymbol{\eta}^j\|^2)$. Since $n_d$ is finite and since $E\{\|\mathbf{H}\|^2\} = \frac{1}{n_d} \sum_{j=1}^{n_d} \|\boldsymbol{\eta}^j\|^2 = \nu$ is also finite, we have $\sup_{j=1, \dots, n_d} \|\boldsymbol{\eta}^j\| = c_\eta < +\infty$. Therefore, $\|g(\mathbf{y})\| \le \frac{1}{s^2} \hat{s}_s c_\eta + \frac{1}{s^2} \|\mathbf{y}\|$ and $\|g(\mathbf{y})\|^4 \le 8(\frac{1}{s^2} \hat{s}_s c_\eta)^4 + 8(\frac{1}{s^2} \|\mathbf{y}\|)^4$. Since $E\{\|\mathbf{Y}\|\} < +\infty$, $E\{\|\mathbf{Y}\|^2\} < +\infty$, and $E\{\|\mathbf{Y}\|^4\} < +\infty$, we have proven the second and third parts of Eq. (4.6). Since, for all $j$, $E\{g_j(\mathbf{Y})^2\} \le E\{\hat{g}(\mathbf{Y})^2\} < +\infty$ yields the first part of Eq. (4.6).

(iii) It can be seen that $h(\mathbf{y}) = \nabla^2 f(\mathbf{y}) = -\|g(\mathbf{y})\|^2 + \xi(\mathbf{y})^{-1} \nabla^2 \xi(\mathbf{y})$, and $\xi(\mathbf{y})^{-1} \nabla^2 \xi(\mathbf{y}) = -\frac{\nu}{s^2} + \frac{1}{n_d} \frac{1}{s^4 \xi(\mathbf{y})} \sum_{j=1}^{n_d} \|\mathbf{y} - \hat{s}_s \boldsymbol{\eta}^j\|^2 \exp(-\frac{1}{2s^2} \|\mathbf{y} - \hat{s}_s \boldsymbol{\eta}^j\|^2)$, which allows us to write $|\xi(\mathbf{y})^{-1} \nabla^2 \xi(\mathbf{y})| \le \frac{\nu}{s^2} + \frac{1}{s^4} \|\mathbf{y}\|^2 + \frac{1}{n_d} \frac{1}{s^2} \frac{\hat{s}_s^2}{s^2} \frac{1}{\xi(\mathbf{y})} \sum_{j=1}^{n_d} \|\boldsymbol{\eta}^j\|^2 \exp(-\frac{1}{2s^2} \|\mathbf{y} - \hat{s}_s \boldsymbol{\eta}^j\|^2)$. Since $\sup_j \|\boldsymbol{\eta}^j\|^2 = \tilde{c}_\eta < +\infty$, we obtain $|\xi(\mathbf{y})^{-1} \nabla^2 \xi(\mathbf{y})| \le \frac{\nu}{s^2} + \frac{1}{s^4} \frac{\hat{s}_s^2}{s^2} \tilde{c}_\eta + \frac{1}{s^4} \|\mathbf{y}\|^2$. It can then be deduced that $|h(\mathbf{y})| \le \|g(\mathbf{y})\|^2 + \frac{\nu}{s^2} + \frac{1}{s^4} \frac{\hat{s}_s^2}{s^2} \tilde{c}_\eta + \frac{1}{s^4} \|\mathbf{y}\|^2$ and using (ii) yields $|h(\mathbf{y})| \le c_0 + c_1 \|\mathbf{y}\| + c_2 \|\mathbf{y}\|^2$ with $c_0 > 0$, $c_1 > 0$, and $c_2 > 0$. We then have proven Eq. (4.7) because $E\{\|\mathbf{Y}\|^r\} < +\infty$ for $r = 0, 1, \dots, 4$.

## 4.2. Polynomial chaos expansion in Gaussian space of potential $\mathcal{V}$

We begin with Definition 2 of the normalized Hermite polynomials on $\mathbb{R}^\nu$ as a Hilbert basis of $\mathbb{H}$. Then, Proposition 2 yields the polynomial chaos expansions in $\mathbb{H}$ on the normalized Hermite polynomials for functions $f$, $g$, and $h$. Finally, Proposition 3 provides the Polynomial chaos expansion in Gaussian space for the potential $\mathcal{V}$.

**Definition 2 (Normalized Hermite polynomials on $\mathbb{R}^\nu$ as a Hilbert basis of $\mathbb{H}$).** (i) Normalized Hermite polynomials on $\mathbb{R}$. Let $\{H_n(\mathbf{y}), n \in \mathbb{N}\}$ be the Hermite polynomials on $\mathbb{R}$ related to the canonical normal measure $p_Y(\mathbf{y}) d\mathbf{y}$ on $\mathbb{R}$ where $p_Y(\mathbf{y}) = (2\pi)^{-1/2} \exp(-y^2/2)$, whose the first three polynomials are $H_0(\mathbf{y}) = 1$, $H_1(\mathbf{y}) = \mathbf{y}$, and $H_2(\mathbf{y}) = y^2 - 1$. It is well known that $\{\psi_n(\mathbf{y}) = H_n(\mathbf{y})/\sqrt{n!}, n \in \mathbb{N}\}$ is a Hilbert basis of the Hilbert space of $L^2(\mathbb{R}; p_Y(\mathbf{y}) d\mathbf{y})$ (see for instance [39, 40, 41, 18]).

(ii) Normalized Hermite polynomials on $\mathbb{R}^\nu$. Let $\boldsymbol{\alpha} = (\alpha_1, \dots, \alpha_\nu)$ be the multi-index in $\mathbb{N}^\nu$. Therefore, the family of normalized Hermite polynomials $\{\psi_{\boldsymbol{\alpha}}, \boldsymbol{\alpha} \in \mathbb{N}^\nu\}$, defined by

$$ \psi_{\boldsymbol{\alpha}}(\mathbf{y}) = \frac{H_{\alpha_1}(y_1)}{\sqrt{\alpha_1!}} \times \dots \times \frac{H_{\alpha_\nu}(y_\nu)}{\sqrt{\alpha_\nu!}} , \quad \forall \mathbf{y} = (y_1, \dots, y_\nu) \in \mathbb{R}^\nu, \quad (4.8) $$

constitutes a Hilbert basis of $\mathbb{H}$, which implies that

$$ \langle \psi_{\boldsymbol{\alpha}}, \psi_{\boldsymbol{\beta}} \rangle_{\mathbb{H}} = E\{\psi_{\boldsymbol{\alpha}}(\mathbf{Y}) \psi_{\boldsymbol{\beta}}(\mathbf{Y})\} = \delta_{\boldsymbol{\alpha}\boldsymbol{\beta}} , \quad \forall \boldsymbol{\alpha} \in \mathbb{N}^{\nu}, \forall \boldsymbol{\beta} \in \mathbb{N}^{\nu}. \quad (4.9) $$

Let $|\boldsymbol{\alpha}| = \alpha_1 + \dots + \alpha_\nu$ be the degree of $\psi_{\boldsymbol{\alpha}}$. For $|\boldsymbol{\alpha}| = 0$ we have $\boldsymbol{\alpha}^{(0)} = (0, \dots, 0)$ and $\psi_{\boldsymbol{\alpha}^{(0)}}(\mathbf{y}) = 1$. It can then be deduced that, for all $\boldsymbol{\alpha}$ such that $|\boldsymbol{\alpha}| \ge 1$, $E\{\psi_{\boldsymbol{\alpha}}(\mathbf{Y})\} = 0$.
**Proposition 2 (Polynomial chaos expansions in $\mathbb{H}$ for functions $f$, $g$, and $h$).** Let $f$, $g = (g_1, \dots, g_v)$, and $h$ be the functions on $\mathbb{R}^v$ defined in Proposition 1.

(i) Function $f$ in $\mathbb{H}$ admits the polynomial chaos expansion $f = \sum_{\alpha \in \mathbb{N}^v} f_\alpha \psi_\alpha$, which is convergent in $\mathbb{H}$. The coefficients $f_\alpha = \langle f, \psi_\alpha \rangle_{\mathbb{H}}$ are real, $f_{\alpha^{(0)}} = E\{f(Y)\}$, and $\|f\|_{\mathbb{H}}^2 = \sum_{\alpha \in \mathbb{N}^v} f_\alpha^2 < +\infty$.

(ii) For all $j = 1, \dots, v$, function $g_j = \partial_j f$ in $\mathbb{H}$ admits the polynomial chaos expansion $g_j = \sum_{\alpha \in \mathbb{N}^v} g_{j,\alpha} \psi_\alpha$, which is convergent in $\mathbb{H}$, where the real coefficients $g_{j,\alpha}$ are defined by

$$g_{j,\alpha} = \sqrt{\alpha_j + 1} f_{\alpha_1, \dots, \alpha_{j+1}, \dots, \alpha_v}, \quad \forall \alpha = (\alpha_1, \dots, \alpha_j, \dots, \alpha_v) \in \mathbb{N}^v. \quad (4.10)$$

We have $g_{j,\alpha^{(0)}} = E\{g_j(Y)\}$ and $\|g_j\|_{\mathbb{H}}^2 = \sum_{\alpha \in \mathbb{N}^v} g_{j,\alpha}^2 < +\infty$.

(iii) Function $h = \sum_{j=1}^v \partial_j^2 f$ in $\mathbb{H}$ admits the following polynomial chaos expansion $h = \sum_{\alpha \in \mathbb{N}^v} h_\alpha \psi_\alpha$, which are convergent in $\mathbb{H}$, where the real coefficients $h_\alpha$ are defined by

$$h_{\alpha} = \sum_{j=1}^{v} \sqrt{(\alpha_j + 1)(\alpha_j + 2)} f_{\alpha_1, \dots, \alpha_{j+2}, \dots, \alpha_v}, \quad \forall \alpha = (\alpha_1, \dots, \alpha_j, \dots, \alpha_v) \in \mathbb{N}^v. \quad (4.11)$$

We have $h_{\alpha^{(0)}} = E\{h(Y)\}$ and $\|h\|_{\mathbb{H}}^2 = \sum_{\alpha \in \mathbb{N}^v} h_{\alpha}^2 < +\infty$.

**PROOF.** (Proposition 2). We use simultaneously Proposition 1 and Definition 1.

(i) Since $f$ is in $\mathbb{H}$ and since $\{\psi_\alpha, \alpha \in \mathbb{N}^v\}$ is a Hilbert basis of $\mathbb{H}$, this yields (i).

(ii) Since $g_j$ is in $\mathbb{H}$, for all $j \in \{1, \dots, v\}$, $f$ has a partial derivative with respect to $y_j$ in $\mathbb{H}$. Consequently, $g_j = \partial_j f = \partial_j \sum_{\alpha \in \mathbb{N}^v} f_\alpha \psi_\alpha = \sum_{\alpha \in \mathbb{N}^v} f_\alpha \partial_j \psi_\alpha$. It is known that $\partial_j H_{\alpha_j}(y_j) = \alpha_j H_{\alpha_j-1}(y_j)$. Therefore, for $\alpha_j \ge 1$, we have $\partial_j H_{\alpha_j}(y_j)/\sqrt{\alpha_j!} = \sqrt{\alpha_j} H_{\alpha_j-1}(y_j)/\sqrt{(\alpha_j-1)!}$ and for $\alpha_j = 0$, we have $\partial_j H_{\alpha_j}(y_j)/\sqrt{\alpha_j!} = 0$. Function $g_j$ can then be rewritten as

$$g_j = \sum_{\alpha_1 \ge 0} \dots \sum_{\alpha_j \ge 1} \dots \sum_{\alpha_v \ge 0} \sqrt{\alpha_j} f_{\alpha_1, \dots, \alpha_{j-1}, \dots, \alpha_v} \frac{H_{\alpha_1}}{\sqrt{\alpha_1!}} \times \dots \times \frac{H_{\alpha_{j-1}}}{\sqrt{(\alpha_{j-1})!}} \times \dots \times \frac{H_{\alpha_v}}{\sqrt{\alpha_v!}}$$

Performing the change of index $\alpha_j = \hat{\alpha}_j + 1$ and renaming $\hat{\alpha}_j$ to $\alpha_j$ allow (ii) to be proven.

(iii) Since $h = \nabla^2 f = \partial_1^2 f + \dots + \partial_v^2 f$ is in $\mathbb{H}$, the function $\partial_j^2 f$ exists in $\mathbb{H}$. The proof of (iii) is then similar to that of (ii) and is left to the reader.

**Proposition 3 (Polynomial chaos expansions in Gaussian space of potential $V$).** The potential function $V$, defined by Eq. (3.20) or even by Eq. (3.30), is in $\mathbb{H}$ and its polynomial chaos expansion in Gaussian space is written for all $y$ in $\mathbb{R}^v$ as

$$V(y) = \sum_{\alpha' \in \mathbb{N}^v} V_{\alpha'} \psi_{\alpha'}(y), \quad E\{\|V(Y)\|^2\} = \|V\|_{\mathbb{H}}^2 = \sum_{\alpha' \in \mathbb{N}^v} V_{\alpha'}^2 < +\infty. \quad (4.12)$$

The polynomial chaos expansion of $V(y)$ can also be rewritten as

$$V(y) = \frac{1}{8} \sum_{\alpha \in \mathbb{N}^v} \sum_{\beta \in \mathbb{N}^v} \mathfrak{G}_{\alpha\beta} \psi_{\alpha}(y) \psi_{\beta}(y) + \frac{1}{4} \sum_{\alpha \in \mathbb{N}^v} h_{\alpha} \psi_{\alpha}(y), \quad (4.13)$$

in which $\mathfrak{G}_{\alpha\beta} = \sum_{j=1}^v g_{j,\alpha} g_{j,\beta}$ and where $g_{j,\alpha}$ and $h_\alpha$ are defined by Eqs. (4.10) and (4.11). For any $\alpha'$ in $\mathbb{N}^v$, the real-valued coefficient $V_{\alpha'}$ can be calculated by

$$V_{\alpha'} = \frac{1}{8} \sum_{\alpha \in \mathbb{N}^v} \sum_{\beta \in \mathbb{N}^v} \mathfrak{G}_{\alpha\beta} E\{\psi_{\alpha}(Y) \psi_{\beta}(Y) \psi_{\alpha'}(Y)\} + \frac{1}{4} h_{\alpha'} \quad (4.14)$$

The mean value of the random variable $V(Y)$ is such that

$$E\{V(Y)\} = \frac{1}{8} \sum_{\alpha \in \mathbb{N}^v} \mathfrak{G}_{\alpha\alpha} + \frac{1}{4} h_{\alpha^0}, \quad \alpha^0 = (0, \dots, 0) \in \mathbb{N}^v. \quad (4.15)$$
PROOF. (Proposition 3). Eq. (3.30) and notations of Proposition 1 yield $\mathcal{V}(y) = \frac{1}{8} \|g(y)\|^2 + \frac{1}{4} h(y)$. Let us prove that $\mathcal{V}$ is in $\mathbb{H}$, that is to say, $\mathcal{V}(Y)$ is a second-order random variable. Since $\mathcal{V}$ is continuous on $\mathbb{R}^v$, $\mathcal{V}(Y)$ is a real-valued random variable. We have $\mathcal{V}(y)^2 \le \frac{1}{32} \|g(y)\|^4 + \frac{1}{8} h(y)^2$, which shows that $E\{\mathcal{V}(Y)^2\} \le \frac{1}{32} E\{\|g(Y)\|^4\} + \frac{1}{8} E\{h(Y)^2\}$. Using Eq. (4.6) yields $E\{\mathcal{V}(Y)^2\} = \|\mathcal{V}\|_{\mathbb{H}}^2 < +\infty$, which proves that $\mathcal{V}(Y)$ is of second order and that $\mathcal{V}$ is in $\mathbb{H}$. Consequently, we obtain the three parts of Eq. (4.12). Rewriting $\mathcal{V}(y)$ as $\mathcal{V}(y) = \frac{1}{8} \sum_{j=1}^v g_j(y)^2 + \frac{1}{4} h(y)$ and using Proposition 2-(ii) and (iii) yields Eq. (4.13). Combining the first part of Eq. (4.12) with Eq. (4.13), using Definition 2-(ii), and projecting onto $\psi_{\alpha'}$ yields Eq. (4.14). Finally, Eq. (4.15) is directly deduced from Eq. (4.13) by replacing $y$ by $Y$, taking the mathematical expectation, and using Definition 2-(ii).

**Remark 1 (Choice of the polynomial chaos expansion for $\mathcal{V}$).** The representation defined by Eq.(4.12), which requires the computation of the coefficients defined by Eq.(4.14), is given for informational purposes, as this approach is not suitable for the following two reasons. First, and most importantly, it results in the loss of the underlying algebraic structure (presence of a positive term) that we discussed in Section 3.7-(iii), which is essential to preserve when truncating the representation to a finite number of terms. Second, it necessitates additional computation of the coefficients $\mathcal{V}_{\alpha'}$ using Eq.(4.14), a computation that is not required with the representation defined by Eq.(4.13). In conclusion, we will use Eq. (4.13) going forward, as it is more efficient when truncated to a finite number of terms.

### 4.3. Truncated polynomial chaos expansion of potential $\mathcal{V}$ and convergence criterion

Following the discussion in Remark 1, we consider the representation defined by Eq. (4.13), which arises from the substitution in Eq. (3.30) of the representation $g_j = \sum_{\alpha \in \mathbb{N}^v} g_{j, \alpha} \psi_{\alpha}$ and $h = \sum_{\alpha \in \mathbb{N}^v} h_{\alpha} \psi_{\alpha}$ established in Proposition 2-(ii) and (iii). Here, $g_{j, \alpha}$ and $h_{\alpha}$, which are defined by Eqs. (4.10) and (4.11), are derived from the coefficients $f_{\alpha}$ of the polynomial chaos expansion $f = \sum_{\alpha \in \mathbb{N}^v} f_{\alpha} \psi_{\alpha}$ established in Proposition 2-(i). Consequently, we introduce the truncated polynomial chaos expansion of $f$ and we then deduce the truncated polynomial chaos expansion of $g_j$ and $h$.

**Definition 3 (Truncated polynomial chaos expansions of $f$, $g_j$, and $h$).** Let $\mu$ be a finite integer such that $\mu \ge 2$. Let $\kappa$ be an integer equal to 0, 1, or 2. We define the subset $A_{\kappa, \mu}$ of $\mathbb{N}^v$ such that

$$A_{\kappa, \mu} = \{ \alpha \in \mathbb{N}^v \mid 0 \le |\alpha| \le \mu - \kappa \}. \quad (4.16)$$

The number of multi-indices in the set $A_{\kappa, \mu}$ is $|A_{\kappa, \mu}| = \frac{(\nu+\mu-\kappa)!}{\nu!(\mu-\kappa)!}$. The truncated polynomial chaos expansions $g_j^{A_{\kappa, \mu}}$ of $g_j$ for $j = 1, \dots, \nu$, and $h^{A_{\kappa, \mu}}$ of $h$, based on the truncated polynomial chaos expansion $f^{A_{\kappa, \mu}}$ of $f$ limited to the maximum degree $\mu$, along with their norms in $\mathbb{H}$, are defined as follows,

$$f^{A_{0,\mu}} = \sum_{\alpha \in A_{0,\mu}} f_{\alpha} \psi_{\alpha} , \quad \|f^{A_{0,\mu}}\|_{\mathbb{H}}^2 = \sum_{\alpha \in A_{0,\mu}} f_{\alpha}^2 \quad (4.17)$$

$$g_j^{A_{1,\mu}} = \sum_{\alpha \in A_{1,\mu}} g_{j,\alpha} \psi_{\alpha} , \quad \|g_j^{A_{1,\mu}}\|_{\mathbb{H}}^2 = \sum_{\alpha \in A_{1,\mu}} g_{j,\alpha}^2 , \quad j = 1, \dots, \nu \quad (4.18)$$

$$h^{A_{2,\mu}} = \sum_{\alpha \in A_{2,\mu}} h_{\alpha} \psi_{\alpha} , \quad \|h^{A_{2,\mu}}\|_{\mathbb{H}}^2 = \sum_{\alpha \in A_{2,\mu}} h_{\alpha}^2 \quad (4.19)$$

The truncated polynomial chaos expansion $\mathcal{V}^{\mu}$ of $\mathcal{V}$, based on the truncated polynomial chaos expansion $f^{A_{0,\mu}}$ of $f$ limited to the maximum degree $\mu$ is defined by,

$$\mathcal{V}^{\mu}(y) = \frac{1}{8} \sum_{\alpha \in A_{1,\mu}} \sum_{\beta \in A_{1,\mu}} \mathfrak{G}_{\alpha\beta} \psi_{\alpha}(y) \psi_{\beta}(y) + \frac{1}{4} \sum_{\alpha \in A_{2,\mu}} h_{\alpha} \psi_{\alpha}(y), \quad (4.20)$$

**Proposition 4 (Convergence criterion of the sequences in $\mu$).** (i) The sequences of approximations $\{f^{A_{0,\mu}}\}_{\mu}$, $\{g_j^{A_{1,\mu}}\}_{\mu}$ for $j = 1, \dots, \nu$, and $\{h^{A_{2,\mu}}\}_{\mu}$ strongly converge to $f$, $g_j$ for $j = 1, \dots, \nu$, and $h$, respectively, as $\mu$ approaches $+\infty$, and we have

$$\lim_{\mu \to +\infty} \|f^{A_{0,\mu}}\|_{\mathbb{H}} = \|f\|_{\mathbb{H}} , \quad \lim_{\mu \to +\infty} \|g_j^{A_{1,\mu}}\|_{\mathbb{H}} = \|g_j\|_{\mathbb{H}} , \quad \lim_{\mu \to +\infty} \|h^{A_{2,\mu}}\|_{\mathbb{H}} = \|h\|_{\mathbb{H}}. \quad (4.21)$$

(ii) The sequence of approximation $\{\mathcal{V}^{\mu}\}_{\mu}$ converges to $\mathcal{V}$ as $\mu$ approaches $+\infty$, and we have

$$\lim_{\mu \to +\infty} \|\mathcal{V}^{\mu}\|_{\mathbb{H}} = \|\mathcal{V}\|_{\mathbb{H}}. \quad (4.22)$$
PROOF. (Proposition 4). (i) The strong convergence in $\mathbb{H}$, that is, convergence with respect to the norm in $\mathbb{H}$, is simply a restatement of Proposition 2. It is well known that if a sequence $x_\mu$ in a Hilbert space $\mathbb{H}$ strongly converges to and element $x$ of $\mathbb{H}$, then the norm $\|x_\mu\|_{\mathbb{H}}$ converges to $\|x\|_{\mathbb{H}}$.

(ii) Eq. (4.12) also shows that the sum $\sum_{\alpha'} \mathcal{V}_{\alpha'} \psi_{\alpha'}$ strongly converges to $\mathcal{V}$ in $\mathbb{H}$. However, we previously explained, we will not use the representation $\mathcal{V} = \sum_{\alpha' \in \mathbb{N}^\nu} \mathcal{V}_{\alpha'} \psi_{\alpha'}$ for $\mathcal{V}$, but instead of the representation defined by Eq. (4.13). Therefore, we need to prove that if Eq. (4.21) holds, then Eq. (4.22) follows. In the proof of Proposition 3, we showed that $\mathcal{V}(y)^2 \le \frac{1}{32} \|g(y)\|^4 + \frac{1}{8} h(y)^2$, which allows us to write $E\{\mathcal{V}(Y)^2\} \le \frac{1}{32} E\{\|g(Y)\|^4\} + \frac{1}{8} E\{h(Y)^2\}$. This expression shows that, to complete the proof of Eq. (4.22), we need only prove that there exists a finite positive constant $c$ such that $E\{\|g(Y)\|^4\} \le c (E\{\|g(Y)\|^2\})^2$. In the proof of Proposition 1, we obtained the inequality $\|g(y)\| \le a + b \|y\|$, where $a$ and $b$ are finite positive constants. We then have $(E\{\|g(Y)\|^2\})^2 \le 4a^4 + 4b^4(E\{\|Y\|^2\})^2$ and $E\{\|g(Y)\|^4\} \le 8a^4 + 8b^4 E\{\|Y\|^4\}$. Since $E\{\|Y\|^2\} = m_2 < +\infty$ and $E\{\|Y\|^4\} = m_4 < +\infty$, there exists a finite positive constant $c$ such that $E\{\|g(Y)\|^4\} \le c (E\{\|g(Y)\|^2\})^2$, which concludes the proof.

**Remark 2 (Practical criterion for checking the convergence).** The convergence criterion for each considered sequence cannot directly be based on the strong convergence definition or on the use of the Cauchy sequence in Hilbert space $\mathbb{H}$, which would be numerically expensive. This would require the numerical evaluation of a large number of mathematical expectations for a random function on $\mathbb{R}^\nu$. Therefore, we propose using a criterion based on Eq. (4.21) with Eqs. (4.17) to (4.19), without evaluating the norm of the limit. The criteria we propose involve the study of the functions,

$$ \mu \mapsto \|f^{A_{0,\mu}}\|_{\mathbb{H}} = (\sum_{\alpha \in A_{0,\mu}} f_{\alpha}^2)^{1/2}, \quad \mu \mapsto \|g_j^{A_{1,\mu}}\|_{\mathbb{H}} = (\sum_{\alpha \in A_{1,\mu}} g_j^2)^{1/2}, \quad j = 1, \dots, \nu, \quad \mu \mapsto \|h^{A_{2,\mu}}\|_{\mathbb{H}} = (\sum_{\alpha \in A_{2,\mu}} h_{\alpha}^2)^{1/2}. \quad (4.23) $$

**Remark 3 (About the quantification of error induced by the truncation).** It should be noted that an upper bound on the approximation error induced by the use of a finite representation can be calculated using classical results for approximations in Sobolev spaces. For example, suppose that $f$ is a function in the Sobolev space $\mathbb{H}^s(\mathbb{R}^\nu, p_Y(y) dy)$. Its truncated PCE of order $\mu$ is given by $f^\mu(y) = \sum_{|\alpha| \le \mu} f_\alpha \psi_\alpha(y)$. The truncation error satisfies $\|f - f^\mu\|_{\mathbb{H}}^2 = \sum_{|\alpha| > \mu} f_\alpha^2$, where $\mathbb{H} = L^2(\mathbb{R}^\nu, p_Y(y) dy)$. The decay rate of this error depends on the regularity of $f$. If $f \in \mathbb{H}^s(\mathbb{R}^\nu, p_Y(y) dy)$, then its expansion coefficients decay like $|f_\alpha| \le \tilde{c} (|\alpha|+1)^{-r}$, where $r$ is determined by the smoothness $s$. This type of decay is typical in spectral approximations; smoother functions (larger $s$) have faster decaying coefficients. Note that the relation $r = s + \nu/2$ is standard under natural assumptions of polynomial chaos expansions. Using the decay of coefficients, the error $\|f - f^\mu\|_{\mathbb{H}}^2 \le c \sum_{\alpha>\mu} (|\alpha|+1)^{-2r}$. By approximating this sum with an integral (in high dimensions, with an integration in spherical coordinates in $\mathbb{R}^\nu$), one obtains a bound of the form $\|f - f^\mu\|_{\mathbb{H}}^2 \le c \mu^{-p}$, where $p$ is a positive exponent that depends on both the decay rate $r$ and the space dimension $\nu$. For example, under typical conditions one finds $p = 2r - \nu$ (with the requirement $2r > \nu$ for convergence). The constant $c$ is problem-dependent and cannot be determined analytically, which is standard in such estimates.

## 4.4. Algebraic formulas for computation and convergence analyses

In this Section, we propose a numerical method for the effective calculation of the coefficients $f_\alpha$ of the polynomial chaos expansion. We summarize the calculation process for $g_{j,\alpha}$, $h_\alpha$, and $g_{\alpha\beta}$, and address the questions related to convergence control.

(i) **Computation of an approximation of $f_\alpha$ and convergence criterion for the approximation.** Let $A^{0,\mu}$ be the set of multi-indices $\alpha = (\alpha_1, \dots, \alpha_\nu) \in \mathbb{N}^\nu$, where the maximum degree $\mu$ is fixed. For a given $\alpha$ in $A^{0,\mu}$, Proposition 2-(i) defines the coefficient $f_\alpha$, which we rewrite in terms of a mathematical expectation as $f_\alpha = \langle f, \psi_\alpha \rangle_{\mathbb{H}} = E\{f(Y) \psi_\alpha(Y)\}$. This mathematical expectation is classically estimated using $N$ independent realizations $\{y^l \in \mathbb{R}^\nu, l = 1, \dots, N\}$ of the $\mathbb{R}^\nu$-valued normalized Gaussian random variable $Y$ (see Definition 1), which yields the following approximation $f_\alpha^N$ of $f_\alpha$,

$$ f_{\alpha}^{N} = \frac{1}{N} \sum_{l=1}^{N} f(y^l) \psi_{\alpha}(y^l) , \quad \alpha \in A^{0,\mu}, \quad (4.24) $$
in which $f(y^l) = \log(p_H(y^l))$, computed using Eq. (2.3). Let $\mathbf{Y}^1, \dots, \mathbf{Y}^N$ be $N$ independent copies of the random vector $\mathbf{Y}$. Let $F_{\alpha}^N$ be the real-valued random variable defined on probability space $(\Theta, \mathcal{T}, \mathcal{P})$ such that $F_{\alpha}^N = \frac{1}{N} \sum_{l=1}^N f(\mathbf{Y}^l) \psi_{\alpha}(\mathbf{Y}^l)$. Then (see, for instance, [42]), the sequence of real-valued random variables $\{F_{\alpha}^N\}_N$ on $(\Theta, \mathcal{T}, \mathcal{P})$ is convergent in probability to $f_{\alpha}$,

$$ \forall \varepsilon > 0, \quad \lim_{N \to +\infty} \mathcal{P}\{|F_{\alpha}^N - f_{\alpha}| \ge \varepsilon\} = 0. \qquad (4.25) $$

There is also the almost sure convergence, thanks to the strong law of large numbers

$$ \mathcal{P}\{\lim_{N \to +\infty} F_{\alpha}^N = f_{\alpha}\} = 1. \qquad (4.26) $$

It should be noted that the speed of convergence is proportional to $1/\sqrt{N^2} = 1/N$ and is independent of dimension $v$. The quantification of the approximation error could traditionally be estimated using the central limit theorem, which involves the variance of the estimator (see, for instance, [42, 33, 43]).

Concerning the convergence with respect to $\mu$, it is understood that the value of $N$, determined for a convergence tolerance based on the criterion defined by Eq.(4.25), is valid for all multi-indices $\alpha$ in $A^{0,\mu}$, and therefore $N$ does not depend on $\alpha$. However, since $A^{0,\mu}$ depends on $\mu$, $N$ will, *a priori*, depend on $\mu$. Under these conditions, the convergence analysis with respect to $\mu$ is carried out using the first equation in Eq.(4.23), for which the function $\mu \mapsto \|f^{A^{0,\mu}}\|_{\mathcal{H}}$ will depend on the value of $N$ that has been determined. This means that the value of $N$ is updated for each value of $\mu$ in the convergence analysis with respect to $\mu$.

(ii) *Computation of the corresponding approximations of $g_{j,\alpha}$, $\mathfrak{G}_{\alpha\beta}$, $h_{\alpha}$, and $\mathcal{V}$. It is assumed that $\mu$ and the associated value of $N$ are fixed to the values identified in (i) above, for which all the approximations $f_{\alpha}^N$ of the coefficient $f_{\alpha}$ are computed for all $\alpha$ in $A^{0,\mu}$. The corresponding approximations of the coefficients $g_{j,\alpha}$ for $j = 1, \dots, v$, $\mathfrak{G}_{\alpha\beta}$, and $h_{\alpha}$ defined by Eq. (4.10), in Proposition 3, and by (4.11), respectively, are thus computed by the formulas:*

$$ g_{j,\alpha}^N = \sqrt{\alpha_j + 1} f_{\alpha_1, \dots, \alpha_j+1, \dots, \alpha_v}^N, \quad \forall \alpha = (\alpha_1, \dots, \alpha_j, \dots, \alpha_v) \in A^{1,\mu}, \quad j = 1, \dots, v, \qquad (4.27) $$

where $A^{1,\mu}$ is defined by Eq. (4.16) for $\kappa = 1$,

$$ \mathfrak{G}_{\alpha\beta}^N = \sum_{j=1}^{v} g_{j,\alpha}^N g_{j,\beta}^N, \quad \forall \alpha \in A^{1,\mu}, \quad \forall \beta \in A^{1,\mu}, \qquad (4.28) $$

and

$$ h_{\alpha}^{N} = \sum_{j=1}^{v} \sqrt{(\alpha_j + 1)(\alpha_j + 2)} f_{\alpha_1, \dots, \alpha_{j+2}, \dots, \alpha_v}^{N}, \quad \forall \alpha = (\alpha_1, \dots, \alpha_j, \dots, \alpha_v) \in A^{2, \mu}. \qquad (4.29) $$

where $A^{2,\mu}$ is defined by Eq. (4.16) for $\kappa = 2$. The corresponding approximation of $\mathcal{V}^{\mu}$ defined by Eq. (4.20) are thus computed by

$$ \mathcal{V}^{\mu,N}(\mathbf{y}) = \frac{1}{8} \sum_{\alpha \in A_{1,\mu}} \sum_{\beta \in A_{1,\mu}} \mathfrak{G}_{\alpha\beta}^N \psi_{\alpha}(\mathbf{y}) \psi_{\beta}(\mathbf{y}) + \frac{1}{4} \sum_{\alpha \in A_{2,\mu}} h_{\alpha}^N \psi_{\alpha}(\mathbf{y}), \quad \forall \mathbf{y} \in \mathbb{R}^v. \qquad (4.30) $$

### 4.5. Adaptation of the truncated polynomial chaos expansion of potential $\mathcal{V}$ to a quantum computing formulation

According to the explanations given in Sections 1 and 3.7, we need to rewrite the right-hand side of Eq. (4.29) as a sum of monomials, which is straightforward.

Let $\psi_{\alpha}$ be the normalized Hermite polynomial of multi-index $\alpha$ defined by Eq. (4.8), where $\alpha$ belongs to $A_{1,\mu}$ or $A_{2,\mu}$ as defined by Eq. (4.16). For $j$ fixed in $\{1, \dots, v\}$, let $\alpha_j$ be the component $j$ of $\alpha = (\alpha_1, \dots, \alpha_v)$, and the Hermite polynomial $H_{\alpha_j}(y_j)$ for the component $y_j$ of $\mathbf{y} = (y_1, \dots, y_v)$ is written as $H_{\alpha_j}(y_j) = \sum_{m_j=0}^{\alpha_j} a_{m_j, \alpha_j} y_j^{m_j}$. The coefficients $a_{m_j, \alpha_j}$ are known. For instance, if $\alpha_j = 2$, then $H_{\alpha_j}(y_j) = y_j^2 - 1$, and consequently, $a_{0,2} = -1$, $a_{1,2} = 0$, and $a_{2,2} = 1$. However, the maximum value of $\alpha_j$ is $\mu - 1$ if $\alpha \in A_{1,\mu}$ and $\mu - 2$ if $\alpha \in A_{2,\mu}$. To avoid dependence on the upper bound
$\alpha_j$ in the sum that defines $H_{\alpha_j}(y_j)$, which allows us to simplify the algebraic representation, we complete, if necessary, by adding zero coefficients. For $\kappa = 1$ or 2, we therefore have,

$$ \forall \alpha \in A_{\kappa, \mu}, \quad H_{\alpha_j}^{(\kappa)}(y_j) = \sum_{m_j=0}^{\mu-\kappa} a_{m_j, \alpha_j}^{(\kappa)} y_j^{m_j}, \qquad (4.31) $$

with $a_{m_j, \alpha_j}^{(\kappa)} = a_{m_j, \alpha_j}$ if $m_j \le \alpha_j$ and, if $\alpha_j+1 > \mu-\kappa$, then $a_{m_j, \alpha_j}^{(\kappa)} = 0$ for all $m_j \in \{\alpha_j+1, \dots, \mu-\kappa\}$. Let $\mathbf{m} = (m_1, \dots, m_v)$ and $\mathbf{m}' = (m'_1, \dots, m'_v)$, the vectors of integers belonging to the set $\mathcal{M}_\kappa = \{0, 1, \dots, \mu-\kappa\}^v$. Using Eqs. (4.8) and (4.31), Eq. (4.30) can be rewritten as

$$ \mathcal{V}^{\mu,N}(\mathbf{y}) = \frac{1}{8} \sum_{\mathbf{m} \in \mathcal{M}_1} \sum_{\mathbf{m}' \in \mathcal{M}_1} \hat{\mathcal{G}}_{\mathbf{m}\mathbf{m}'}^N \mathbf{y}^\mathbf{m} \mathbf{y}^{\mathbf{m}'} + \frac{1}{4} \sum_{\mathbf{m} \in \mathcal{M}_2} \hat{h}_m^N \mathbf{y}^\mathbf{m}, \quad \forall \mathbf{y} \in \mathbb{R}^v. \qquad (4.32) $$

In Eq. (4.31), $\mathbf{y}^\mathbf{m} = y_1^{m_1} \times \dots \times y_v^{m_v}$, and for all $\mathbf{m}$ and $\mathbf{m}'$ in $\mathcal{M}_1$, we have

$$ \hat{\mathcal{G}}_{\mathbf{m}\mathbf{m}'}^N = \sum_{\alpha \in A_{1,\mu}} \sum_{\beta \in A_{1,\mu}} \mathfrak{G}_{\alpha\beta}^N \frac{a_{\mathbf{m},\alpha}^{(1)}}{\sqrt{\alpha!}} \frac{a_{\mathbf{m}',\alpha}^{(1)}}{\sqrt{\beta!}}, \qquad (4.33) $$

where $\alpha! = \alpha_1! \times \dots \times \alpha_v!$ and $a_{\mathbf{m},\alpha}^{(1)} = a_{m_1,\alpha_1}^{(1)} \times \dots \times a_{m_v,\alpha_v}^{(1)}$. For all $\mathbf{m}$ in $\mathcal{M}_2$, we have

$$ \hat{h}_m^N = \sum_{\alpha \in A_{2,\mu}} \hat{h}_\alpha^N \frac{a_{\mathbf{m},\alpha}^{(2)}}{\sqrt{\alpha!}}. \qquad (4.34) $$

## 5. Numerical validation of the software code for polynomial-chaos-expansion coefficient computation

A software code has been developed based on the algebraic formulas for computation and convergence analyses, which we have presented and summarized in Section 4.4. However, such a validation is not easy if the probability measure $p_H(\mathbf{y})$, $d\mathbf{y}$ is arbitrary, because in that case, we do not have a reference to compare with. Given that our purpose is to validate the software code, we can choose any probability measure, particularly one for which we have a reference. We opted $\mathbf{H}$ to be a $\mathbb{R}^v$-valued normalized Gaussian random variable. Certainly, there are choices of probability distributions for which the truncation order would be very high and so not accessible to quantum computers with the proposal function. Nevertheless, the formulation is general for any non-Gaussian probability distribution, and the convergence of the approximation is mathematically proven in the paper. Thus, a numerical simulation of a non-Gaussian probability distribution would only bring numerical complexity at this stage. The choice of a normalized random variable is justified by the fact that Eq. (2.2) must be satisfied. We will begin this section by providing the reference algebraic formulas for this Gaussian case, and then we will continue by presenting the numerical results. However, we must first specify the parameterization of the numerical analysis that we are going to present.

(i) Since $\mathbf{H}$ is chosen as Gaussian for this validation, we can skip the convergence analysis with respect to $\mu$ for the following reason. For the considered validation case, the reference probability density function is $p_H^{\text{ref}}(\mathbf{y}) = (2\pi)^{-v/2} \exp(-\|\mathbf{y}\|^2/2)$ that is independent of $n_d$, while its GKDE approximation is given by Eq. (2.3) and depends on $n_d$. This case defined by $p_H^{\text{ref}}$ will be named the Gaussian reference case. We then have

$$ f(\mathbf{y}) = \log(p_H^{\text{ref}}(\mathbf{y})) = -\frac{\nu}{2} \log(2\pi) - \frac{1}{2} \|\mathbf{y}\|^2, \qquad (5.1) $$

which shows that $\mu = 2$ could be sufficient. Under this condition, we have fixed the value $\mu = 4$.

(ii) The value of $\mu$ being fixed, in order to minimize the number of figures presented for the convergence analyses with respect to $N$ we will limit the calculated quantities to the functions defined by Eq. (4.23). Since $\mu$ is fixed and since these functions depend on $N$, we adapt the definition and notation of these functions as follows,

$$ N \mapsto \mathcal{E}_f(N) = \{\sum_{\alpha \in A_{0,\mu}} (f_{\alpha}^N)^2\}^{1/2}, \qquad (5.2) $$

$$ N \mapsto \mathcal{E}_g(N) = \{\sum_{\alpha \in A_{1,\mu}} \sum_{j=1}^{\nu} (g_{j,\alpha}^N)^2\}^{1/2}, \qquad (5.3) $$

$$ N \mapsto \mathcal{E}_h(N) = \{\sum_{\alpha \in A_{2,\mu}} (h_{\alpha}^N)^2\}^{1/2}. \qquad (5.4) $$
This choice of functions defined by Eqs. (5.2) to (5.4) allows us to control the convergence of all the coefficients of the polynomial chaos expansions. For this Gaussian reference case, we will also use the GKDE representation. In that case, these functions also depend on $n_d$, and we will consider these functions as a family of functions indexed by the parameter $n_d$.

### 5.1. Algebraic formulas for the Gaussian reference case

Using Eq. (4.30), it can be seen that $g(y) = \nabla f(y) = -y$ and $h(y) = \nabla^2 f(y) = -v$. Using the classical formulas: $E\{Y_j\} = 0$, $E\{Y_j^2\} = 1$, and $E\{Y_j^4\} = 3$, it can easily proven the following results,

$$E_f^{\text{ref}} = \|f\|_{H} = \{\sum_{\alpha \in \mathbb{N}^{\nu}} f_{\alpha}^{2}\}^{1/2} = \{v/2 + (1+2\log(2\pi) + (\log(2\pi))^{2})v^{2}/4\}^{1/2}, \quad (5.5)$$

$$E_g^{\text{ref}} = \{\sum_{j=1}^{\nu} \|g_j\|_{H}^{2}\}^{1/2} = \{\sum_{j=1}^{\nu} \sum_{\alpha \in \mathbb{N}^{\nu}} g_{j,\alpha}^{2}\}^{1/2} = \sqrt{v}, \quad (5.6)$$

$$E_h^{\text{ref}} = \|h\|_{H} = \{\sum_{\alpha \in \mathbb{N}^{\nu}} h_{\alpha}^{2}\}^{1/2} = v. \quad (5.7)$$

### 5.2. Numerical validation of the software code

The convergence analysis with respect to $N$ is done for two values of $v$. Fig. 1 concerns the results for $v = 1$, while Fig. 2 concerns the results for $v = 10$.

(i) Gaussian reference case. For the Gaussian reference case, for which the probability density function of $H$ is $p_H^{\text{ref}}$ defined at the beginning of Section 5, the use of formulas defined by Eqs. (5.5) to (5.7) yields $(E_f^{\text{ref}}, E_g^{\text{ref}}, E_h^{\text{ref}}) = (1.585, 1.000, 1.000)$ for $v = 1$ and $(14.364, 3.162, 10.000)$ for $v = 10$. For analyzing the convergence of the polynomial chaos expansions of this Gaussian reference case, the functions $N \mapsto E_f(N), E_g(N)$, and $E_h(N)$ are computed using Eqs. (5.2), (5.3), and (5.4), respectively. The graphs of these functions are displayed in Figs. 1a, 1b, and 1c for $v = 1$, and in Figs. 1d, 1e, and 1f for $v = 10$. For $N = 10^7$, we have $(E_f(10^7), E_g(10^7), E_h(10^7)) = (1.585, 0.999, 0.999)$ for $v = 1$ and $(14.364, 3.171, 9.971)$ for $v = 10$. These values are very close to the reference values that show a complete convergence and a full validation of the software code devoted to the computation of the polynomial-chaos-expansion coefficients based on the theory presented in Section 4.

(ii) Gaussian case with GKDE representation. With respect to software code validation, this case does not provide anything more than the Gaussian reference case presented in (i), because the algorithms for calculating the coefficients of polynomial chaos expansions are independent of the choice of the probability density of $H$. What this analysis shows for this Gaussian case is the influence of the number $n_d$ of points in the training dataset on the estimation of $p_H$ by the GKDE method (given by Eq. (2.3)), compared to the exact expression $p_H^{\text{ref}}$ (which corresponds to $p_H$ as $n_d \to +\infty$). Therefore, for the Gaussian case with the GKDE estimation of $p_H$ defined by Eq. (2.3), we have considered four values of $n_d$, which are 100, 1000, 5000, and 10000. For the convergence analysis of the polynomial chaos expansions, the functions $N \mapsto E_f(N), E_g(N)$, and $E_h(N)$ are computed using Eqs. (5.2), (5.3), and (5.4), respectively. The graphs of these functions are displayed in Figs. 2a, 2b, and 2c for $v = 1$, and in Figs. 2d, 2e, and 2f for $v = 10$. For $N = 10^7$, we have $(E_f(10^7), E_g(10^7), E_h(10^7)) = (1.594, 1.026, 1.027)$ for $v = 1$ and $(15.249, 4.502, 14.148)$ for $v = 10$. These values are close to the reference values for $v = 1$ and slightly further away for $v = 10$. This does not call into question either the theoretical formulation proposed or the validation of the software code. This difference is simply due to the convergence of $p_H$ towards $p_H^{\text{ref}}$ as a function of $n_d$ for a fixed value of $v$. We know that the larger $v$ is, the larger $n_d$ must be. However, in the context of the probabilistic learning on manifolds approach, the eigenfunctions associated with the first eigenvalues of the Fokker-Planck operator that we are looking for are only used to construct a vector basis for the projection of the data, and replacing the exact operator with an approximate operator, constructed with a small training dataset, does not introduce any difficulties (see, for example, [1]).

## 6. Quantum computing formulation

The eigenproblem given in Eq. (3.17) is fully equivalent to a time-independent Schrödinger equation, where $\hat{L}_{FKP}$ acts as the Hamiltonian and $q$ represents the corresponding eigenstate (or wavefunction). Solving this eigenproblem is commonly approached in the quantum mechanics community via projection (using the variational principle) onto a
<figure><img src="image_2.png" /><figcaption>Figure 1: Convergence analysis for the Gaussian reference case, for $v = 1$ (a,b,c) and $v = 10$ (d,e,f). Graph of functions $N \mapsto \mathcal{E}_f^{\text{ref}}(N)$ (a,d), $N \mapsto \mathcal{E}_g^{\text{ref}}(N)$ (b,e), and $N \mapsto \mathcal{E}_h^{\text{ref}}(N)$ (c,f). The horizontal axes are $\log_{10}(N)$.</figcaption></figure>

tensor-product basis spanning the full space $\mathbb{R}^v$ [44, 45, 46, 47]. This approach is feasible only for low dimensionality $v$, as its computational cost scales exponentially with $v$. For large dimensional problems, alternative approximations using specific tensor decompositions are sometimes employed [48, 49, 50], though these schemes do not completely alleviate the exponential scaling. Recent advances in quantum computing offer potential for encoding the wavefunctions $q$ and their operations, even in high-dimensional problems (see for instance that quantum computers now have more than 1000 qubits [51]). This encoding becomes possible by mapping the system, as defined by the operator $\hat{L}_{\text{FKP}}$, onto a quantum computer, thereby representing wavefunctions and operations as qubits. While much of the research community is focused on electronic structure problems [52], new efforts are emerging to adapt other quantum problems for quantum computation [13, 14, 53]. For completeness, we depict here the general approach as follows: (i) choice of a basis, (ii) second quantization, and (iii) mapping onto qubits.

## 6.1. Choice of the basis

The first step is to adopt a finite basis representation (spectral method), for which we can easily rewrite Eq. (3.17) in a second quantized form. In our case, all dimensions are distinguishable, such that operators acting on different coordinates $y_j$ commute, and we must use different basis functions for each dimension. In line with section 3.7, we employ a basis of Hermite polynomials as a working basis. A common practice in the quantum mechanical community is to absorb the measure $p_Y(y) dy = \otimes_{j=1}^v p_Y(y_j) dy_j$ into the basis, where $p_Y$ is defined in Definition 2. We also switch to the physicists’ convention for the Hermite polynomials during the process of building the basis. Furthermore, we perform an affine transformation on each coordinate, $y_j \to \sqrt{2\omega_j}(y_j - y_{j,0})$, to adapt the basis to improve convergence, where $y_{j,0}$ allows for a shift of the coordinate $y_j$, if needed. For coordinate $y_j$, the resulting basis $\{\varphi_n^{(j)}(y_j), n \in \mathbb{N}\}$ is such that

$$ \varphi_n^{(j)}(y_j) = \sqrt{2^n} \psi_n(\sqrt{2\omega_j}(y_j - y_{j,0})) \{p_Y(\sqrt{2\omega_j}(y_j - y_{j,0}))\}^{1/2}, $$

for which the measure is $dy_j$ and where $\psi_n$ is defined in Definition 2. These basis functions are the normalized eigenfunctions of the operator $\omega_j (y_j - y_{j,0})^2 - (1/\omega_j) \partial^2/\partial y_j^2$. By projecting Eq. (3.17) onto this basis, the eigenproblem
<figure><img src="image_3.png" /><figcaption>Figure 2: Convergence analysis for the Gaussian GKDE case, for $v = 1$ (a,b,c) and $v = 10$ (d,e,f). For $n_d = 100$ (dashed line), 1000 (thin line), 5000 (thick line), and 10 000 (line with diamonds), graph of functions $N \mapsto \mathcal{E}_f(N)$ (a,d), $N \mapsto \mathcal{E}_g(N)$ (b,e), and $N \mapsto \mathcal{E}_h(N)$ (c,f). The horizontal axes are $\log_{10}(N)$.</figcaption></figure>

becomes a matrix diagonalization, where the eigenstate $q_n$ is expanded as

$$q_n(\mathbf{y}) = \sum_{\alpha \in M} c_{n,\alpha} \varphi_{\alpha}(\mathbf{y}), \qquad (6.1)$$

where $M = \{\alpha \in \mathbb{N}^v | \alpha_j \le m_j, j = 1, 2, \dots, v\}$, $c_{n,\alpha} \in \mathbb{R}$, and $\varphi_{\alpha}(\mathbf{y}) = \prod_{j=1}^{v} \varphi_{\alpha_j}^{(j)}(y_j)$, where $\{m_j\}_{j=1}^v$ are the maximum index that basis functions of the degree of freedom $y_j$ can have. Therefore, there are $1 + m_j$ basis functions per degree of freedom, and the total number of one-dimensional basis functions, including all degrees of freedom, is $m_T = \sum_{j=1}^{v}(1 + m_j)$. It should be noted that $m_T$ is not the total number of basis functions in the full tensor product space, whose explicit use is infeasible in quantum computing. Instead, quantum algorithms exploit superposition and entanglement to efficiently process high-dimensional Hilbert spaces, a fundamental aspect of quantum advantage. As a result, $m_T$ increases linearly with the number of dimensions, while the total number of basis functions in the full tensor product space increases exponentially with dimensionality.

## 6.2. Second quantization

On our way to express the eigenproblem on a quantum computer, we use the second quantization formalism [54, 55], which is particularly suited for encoding the qubits. In this formalism, all the basis functions can be generated from $\varphi_0^{(j)}(\mathbf{y})$ by employing the so-called creation operator

$$a_j^+ = \frac{1}{\sqrt{2}} (\sqrt{\omega_j} y_j - \frac{1}{\sqrt{\omega_j}} \frac{\partial}{\partial y_j}),$$

with $a_j^+ \varphi_n^{(j)} = \sqrt{n+1} \varphi_{n+1}^{(j)}$ such that $\varphi_{\alpha}(\mathbf{y}) = (\prod_{j=1}^{v} (a_j^+)^{\alpha_j}) \varphi_{\alpha^{(0)}}(\mathbf{y})$ with $\alpha^{(0)} = (0, 0, \dots, 0)$. By use of the creation operator, one can rewrite the solution $q_n$ in the so-called Fock space, where the basis functions are described by occupation numbers that we denote by

$$\varphi_{\alpha} \equiv |\alpha_1 \alpha_2 \dots \alpha_v\rangle = \left(\prod_{j=1}^{v} (a_j^+)^{\alpha_j}\right) |00\dots0\rangle.$$
Similarly, we define the annihilation operator

$$a_j^- = \frac{1}{\sqrt{2}}(\sqrt{\omega_j} y_j + \frac{1}{\sqrt{\omega_j}} \frac{\partial}{\partial y_j}),$$

with $a_j^- \varphi_n^{(j)} = \sqrt{n} \varphi_{n-1}^{(j)}$. Note that $a_j^+$ and $a_j^-$ are adjoint one to the other. Knowing the effect of $a_j^+$ and $a_j^-$ on the basis, we can rewrite all operators involved in Eq. (4.32) using the identities

$$y_j = \frac{1}{\sqrt{2\omega_j}}(a_j^+ + a_j^-) \quad , \quad \frac{\partial}{\partial y_j} = (a_j^- - a_j^+) \sqrt{\omega_j/2}.$$

Then, using the commutator $a_j^- a_k^+ - a_k^- a_j^+ = \delta_{jk}$, one can express $\nabla^2$ and all monomials in the so-called normal ordering (with creation operators on the left of the annihilation operators) as follows,

$$\nabla^2 = \sum_{j=1}^{\nu} \frac{\omega_j}{2} ( (a_j^+)^2 + (a_j^-)^2 - 2a_j^+ a_j^- - 1 ), \qquad (6.2)$$

$$y^{\alpha} = \prod_{j=1}^{\nu} \frac{\alpha_j!}{\sqrt{(2\omega_j)^{\alpha_j}}} \sum_{s=0}^{\lfloor \alpha_j/2 \rfloor} \sum_{p=0}^{\alpha_j - 2s} \frac{(2s-1)!!}{(\alpha_j - 2s - p)!(2s)!p!} (a_j^+)^{\alpha_j - 2s - p} (a_j^-)^p, \qquad (6.3)$$

where $n!! = 1 \times \dots \times (n-4) \times (n-2) \times n$ is the double factorial, and $\lfloor a \rfloor$ returns the integer part of $a$. For example, a quadratic term is given by $y_j^2 = ((a_j^+)^2 + (a_j^-)^2 + 2a_j^+ a_j^- + 1)/(2\omega_j)$. Note that Eq. (6.3) results from the expansion of the operator $y_j^{\alpha_j} = \frac{1}{\sqrt{(2\omega_j)^{\alpha_j}}} (a_j^+ + a_j^-)^{\alpha_j}$. The double sum and the double factorial reflect the combinatorics arising from the commutation relation $a_j^- a_j^+ - a_j^+ a_j^- = 1$. With these expressions at hand, the approximation $\hat{L}_{\text{FKP}}^{\mu,N}$ of the Fokker-Planck operator $\hat{L}_{\text{FKP}}$, as defined by Eq. (3.19), can be expressed in a general form for a truncation at order $\mu$ with $N$ realizations of the Gaussian random variable $Y$,

$$\hat{L}_{\text{FKP}}^{\mu,N} = \sum_{\substack{\alpha \in \mathbb{N}^{\nu} \\ |\alpha| \le 2(\mu-1)}} \sum_{\substack{k \in \mathbb{N}^{\nu} \\ a^{(0)} \le k \le \alpha}} G_{\alpha,k}^{N} \prod_{j=1}^{\nu} (a_j^+)^{\alpha_j} (a_j^-)^{k_j - \alpha_j}, \qquad (6.4)$$

where the coefficients $G_{\alpha,k}^N$ are obtained by substituting the expression from Eq. (6.3) into Eq. (4.32) and combining it with Eq. (6.2) into Eq. (3.19). Employing the form given in Eq. (6.4) is convenient, as the action of the creation and annihilation operators can easily be computed using the fact that

$$(a_j^+)^{\beta} (a_j^-)^{\alpha} |\varphi_n^{(j)}\rangle = \frac{\sqrt{(n+\beta-\alpha)! n!}}{(n-\alpha)!} |\varphi_{n+\beta-\alpha}^{(j)}\rangle.$$

### 6.3. Map onto a system of qubits

Qubits are 2-level quantum systems, expressed as $\{c_0|0\rangle + c_1|1\rangle, c = (c_0, c_1) \in \mathbb{C}^2\}$, that compose the quantum computers. To exploit quantum computing, we need to map the system of interest to the system of qubits. Rewriting the Fokker-Planck operator in a basis of qubits is not unique. Two mappings exist: the “direct” map [56], also named the “uniary” map [14], which we employ here due to its simplicity, and the “compact” map [57], also named “binary” map [14], which reduces the number of required qubits. An alternative approach involves using quantum computers based on bosonic modes (harmonic oscillators) instead of qubits, or hybrid architectures combining oscillators and qubits, where additional mappings and other quantum circuits construction stategies can also be explored [58, 59]. Here, we utilize the direct map [56, 13, 53]: the occupation of $|\varphi_n^{(j)}\rangle$ is represented by the occupation of a corresponding qubit state $|s_{j,k}\rangle$, such that a single qubit describes the occupation of a single basis function for a given degree of freedom, $y_j$. In the direct map, the total number of qubit is given by $m_T$, the total number of basis functions across all degrees of freedom as defined in Section 6.1. Suppose that the different degrees of freedom require an average of $\bar{m}$ basis functions, the number of qubits, here $m_T$, scales as $\nu \times \bar{m}$. The qubit state corresponding to $\varphi_n^{(j)}$, associated with
the degree of freedom $y_j$, is rewritten, assuming that excatly one qubit is in the state $|1\rangle$ at position $n$, while all others are in the state $|0\rangle$:

$$ |\varphi_n^{(j)}\rangle = \underbrace{|(s_{j,0} = 0)\rangle \otimes \dots \otimes |(s_{j,n-1} = 0)\rangle}_{\text{these qubits in state } |0\rangle} \otimes |(s_{j,n} = 1)\rangle \otimes \underbrace{|(s_{j,n+1} = 0)\rangle \otimes \dots \otimes |(s_{j,m_j} = 0)\rangle}_{\text{remaining qubits in state } |0\rangle} . \quad (6.5) $$

By employing this notation, products of creation and annihilation operators on a $n^{th}$ basis function $|\varphi_n^{(j)}\rangle$ take the form

$$ [(a_j^+)^\beta (a_j^-)^\alpha]_n \equiv \frac{1}{4} \frac{\sqrt{(n+\beta-\alpha)!n!}}{(n-\alpha)!} (\hat{X}_{j,n} + i\hat{Y}_{j,n}) \otimes (\hat{X}_{j,n+\beta-\alpha} - i\hat{Y}_{j,n+\beta-\alpha}), \quad (6.6) $$

where $i$ is the imaginary number, $i^2 = -1$. The operators $\hat{X}_{j,n} = (\begin{smallmatrix} 0 & 1 \\ 1 & 0 \end{smallmatrix})$ and $\hat{Y}_{j,n} = (\begin{smallmatrix} 0 & -i \\ i & 0 \end{smallmatrix})$ are Pauli matrices associated with the $k^{th}$ qubit used in the description of the $j^{th}$ degree of freedom $y_j$. The hat indicates the matrix representation in the basis where individual qubits are represented as vectors $|0\rangle \equiv (\begin{smallmatrix} 1 \\ 0 \end{smallmatrix})$ and $|1\rangle \equiv (\begin{smallmatrix} 0 \\ 1 \end{smallmatrix})$. Then, we can build the operator $(a_j^+)^\beta (a_j^-)^\alpha$ as

$$ (a_j^+)^\beta (a_j^-)^\alpha \equiv \frac{1}{4} \sum_{n=\alpha}^{(2m_j+\alpha-\beta+|\alpha-\beta|)/2} \frac{\sqrt{(n+\beta-\alpha)!n!}}{(n-\alpha)!} (\hat{X}_{j,n} + i\hat{Y}_{j,n}) \otimes (\hat{X}_{j,n+\beta-\alpha} - i\hat{Y}_{j,n+\beta-\alpha}). \quad (6.7) $$

Equation (6.7) sums over all valid $n$ to form the full operator in qubit space, ensuring that each $|n\rangle$ can transition to $|n + \beta - \alpha\rangle$ with the usual harmonic-oscillator amplitude.

## 7. Quantum operations and measurements: quantum circuit

The construction of quantum states and the measurement of their observables are fundamental steps in the development of more advanced quantum algorithms. In this section, we provide an overview of these foundational concepts, particularly for readers who may be less familiar with quantum computing.

Manipulating qubits is achieved through quantum operations that are unitary transformations. Hence, the eigen-solution that we seek can be expressed as $|q_n\rangle = U_q|0^{m_T}\rangle$, where $|0^{m_T}\rangle$ represents the vacuum state of $m_T$ qubits (i.e., all qubits initialized in their ground state). Because individual qubits are 2-level systems, we can employ single-qubit transformations from $U(2)$, such as $\exp(-i\theta_x \hat{X})$, $\exp(-i\theta_y \hat{Y})$, and $\exp(-i\theta_z \hat{Z})$, where $\hat{X}$, $\hat{Y}$, and $\hat{Z}$ are the Pauli matrices (noting that $\hat{Z} = (\begin{smallmatrix} 1 & 0 \\ 0 & -1 \end{smallmatrix})$). Transformations acting on pairs, triplets, or more qubits can be built as tensor products of unitaries from $U(2)$, and more complex multi-qubit transformations generally involve entanglement and interactions between qubits. For example, the Pauli group (composed of tensor products of Pauli operators and the identity) forms a complete basis for constructing any multi-qubit unitary transformation [60]. In view of Eq. (6.5), constructing a specific state is equivalent to creating a superposition of collective qubit states by applying suitable unitary transformations to single qubits, pairs of qubits, or larger subsets of qubits. Such transformations are implemented on quantum computers via quantum gates, which are the fundamental building blocks of quantum circuits. Examples of commonly used quantum gates [61, 62] include the Hadamard gate $H$, the phase gate $S$, and the Controlled-NOT gate for two qubits, such that

$$ H = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}, \quad S = \begin{pmatrix} 1 & 0 \\ 0 & i \end{pmatrix}, \quad \text{CNOT} = \begin{pmatrix} \hat{1} & \hat{0} \\ \hat{0} & \hat{X} \end{pmatrix}, $$

where $\hat{1}$ is the 2x2 identity and $\hat{X}$ is the Pauli-x matrix. To construct eigenstates, one needs a set of universal quantum gates that can reconstruct any unitary transformation in the qubit Hilbert space can be approximated to arbitrary precision [61, 60]. By successively applying these gates, one can construct quantum circuits that generate specific quantum states.

The target solutions of the FKP equation span only a subset of the qubits' Hilbert space. This subset is restricted because certain qubit states are not valid representations of the basis states defined in Section 6.1. For instance, qubit states of the form $\dots \otimes |(s_{j,k} = 1)\rangle \otimes \dots \otimes |(s_{j,l} = 1)\rangle \otimes \dots$ do not correspond to any valid basis state. This is due to the
<figure><img src="image_4.png" /><figcaption>Figure 3: General form of a circuit: the left dashed box indicates the trial state construction for $q_n$ in terms of qubits, while the right dashed box indicates the measurement of an observable $A_β$, where $U_β^m$ is a unitary transformation that diagonalizes $A_β$.</figcaption></figure>

constraint that, for each degree of freedom $y_j$, only one qubit can take the value $s_{j,k} = 1$ among the $m_j$ qubits assigned to describe $y_j$. Mathematically, this constraint is expressed as

$$ \sum_{k=0}^{m_j} (1 - \langle \hat{Z}_{j,k} \rangle) / 2 = \sum_{k=0}^{m_j} s_{j,k} = 1, \quad \forall j \in \{1, \dots, v\}. $$

To address this issue, two possible strategies can be employed: i) adding constraints to the optimization process to explicitly enforce the condition $\sum_{k=0}^{m_j} s_{j,k} = 1$ for all $j$, or ii) restricting the set of unitary transformations so that they inherently respect the constraint. This can be achieved by constructing transformations exclusively from the creation and annihilation operators, ensuring that $\hat{X}$ and $\hat{Y}$ are always employed in pairs of the form $\hat{X}_{j,k} \hat{Y}_{j,l} - \hat{Y}_{j,k} \hat{X}_{j,l}$.

Extracting information from quantum states is accomplished through quantum measurements of physical observables of interest. This measurement process is inherently probabilistic and results in the irreversible collapse of the quantum state. For example, assuming that an observable $\hat{A}$ has an eigendecomposition $\hat{A} = \sum_{k=1}^{2^M} |k\rangle a_k \langle k|$, where $M$ is the number of qubits, measuring $\hat{A}$ on a quantum state $|\chi\rangle$ returns an eigenvalue $a_k$ with probability $|\langle k|\chi\rangle|^2$. After the measurement, the quantum state collapses to the corresponding eigenstate $|k\rangle$. As a result, many measurements (realizations) are required to statistically estimate the expectation value $\langle \hat{A} \rangle = \langle \chi | \hat{A} | \chi \rangle$. For example, in quantum chemical applications using the Variational Quantum Eigensolver (VQE) framework, the number of measurements typically required is on the order of $10^9$ [63]. Only observables that are diagonal in the computational (qubit) basis can be directly measured, i.e., observables built as tensor products of $\hat{1}$ and $\hat{Z}$ for each qubit. Hence, it is necessary to decompose the observable as $\hat{A} = \sum_{\beta} A_{\beta}$, where each fragment $A_{\beta}$ can be efficiently transformed into a diagonal form via a suitable unitary transformation $U_{\beta}^m$. This decomposition implies that the state of interest must also be transformed according to the same unitary operation before measuring each fragment. Reconstructing the expectation value $\langle \hat{A} \rangle$ requires combining all measurement outcomes on a classical computer (see, for instance, [64, 65, 66] for efficient measurement strategies). An alternative approach to measurement involves using the Hadamard test [60], which will be discussed in Section 8.2 and serves as an example of a measurement technique.

Considering the two key steps—construction and measurement of a quantum state—a quantum processing unit operates in two corresponding phases: (i) construction and (ii) measurement, as illustrated in Fig. 3. The measurement process must be repeated a sufficient number of times to ensure convergence to the expectation value $\langle \hat{A}_k \rangle$ for a given fragment. Since each measurement irreversibly collapses the quantum state, the entire process, including the state construction, must be repeated for every realization. Then, a classical computer combines the results from all fragments to reconstruct the total expectation value $\langle \hat{A} \rangle$.

## 8. Proposal for eigenstates construction and overlap extraction

As explained in Section 3.8, we have to compute $g_{FKP,j}^m = v_m(\eta^j)$, for all $m \in \{1, \dots, m_{opt}\}$ and for all $j \in \{1, \dots, n_d\}$, that is to say, using Eq. (3.31), to compute $q_m(\eta^j)$. This is achieved by first generating eigenstates as solution of Eq. (3.24), and second by calculating the following overlap $q_n(\eta^j) = \langle \eta^j | q_n \rangle$.
## 8.1. Eigenstates construction

Some of the most extensively studied algorithms for solving eigenproblems are the Variational Quantum Eigen-solver (VQE) [16, 63] and the Quantum Phase Estimation (QPE) [17, 60].

(i) QPE algorithm. The QPE algorithm typically involves the following steps: 1) Preparing an initial quantum state with a significant overlap with the target eigenstates. 2) Decomposing the propagator $\exp(-i\hat{L}_{FKP}t)$ into a sequence of quantum gates for various times $t$. 3) Computing and storing the overlaps between the resulting states and the initial state, also referred to as the autocorrelation function, for the chosen time intervals. 4) Performing a Fourier transform on the autocorrelation function to extract the eigenspectrum and measuring the resulting state to identify eigenstates with significant overlap with the initial state.

(ii) Rationale for presenting the QPE or VQE algorithm. One objective of this section is to illustrate for readers who are not specialists in scientific quantum computing how the essential coupling between a quantum computer and a classical computer operates. Given the complexity and non-trivial nature of the QPE algorithm, we choose to focus on the VQE approach, which can be introduced using simpler concepts. Furthermore, QPE requires substantially longer quantum circuits than VQE [67, 63], making the latter more practical for quantum computers in the Noisy Intermediate-Scale Quantum (NISQ) era, where quantum hardware remains limited [68]. We therefore focus on the VQE approach and provide here a brief overview.

(iii) VQE algorithm. The VQE approach is a hybrid quantum-classical algorithm [69] designed to approximate the lowest, “ground state”, eigensolutions of an eigen equation (see Fig. 4). It is based on the Rayleigh-Ritz variational principle, which states that any normalized trial wavefunction $\tilde{q}(y)$ satisfies the inequality $R(\tilde{q}) \ge \lambda_0 = 0$, where

$$R(\tilde{q}) = \langle \hat{L}_{FKP} \tilde{q}, \tilde{q} \rangle_{L^2} = \int_{\mathbb{R}^V} \tilde{q}(y) (\hat{L}_{FKP} \tilde{q}(y)) dy, \quad \|\tilde{q}\|_{L^2} = \int_{\mathbb{R}^V} \tilde{q}(y)^2 dy = 1,$$

is the Rayleigh quotient of the symmetric and positive operator $\hat{L}_{FKP}$ (see Sections 3.5 and 3.6). Equality occurs when $\tilde{q}(y) = q_0(y)$ is the ground-state wavefunction. Algorithm VQE combines the capabilities of quantum and classical computers to iteratively refine a parameterized trial wavefunction in order to minimize $R(\tilde{q})$. The quantum computer is used to construct the trial wavefunction and estimate the Rayleigh quotient, while the classical computer optimizes the parameters of the wavefunction.

(iii-1) Initial guess, assumed to be a good approximation to the exact solution. The first step in VQE is to define an initial guess $\tilde{q}^{init}(y)$, which serves as an approximation to the exact eigenstate. A good initial guess reduces the number of iterations required for convergence. The construction of $\tilde{q}^{init}(y)$ is performed on a classical computer, using a mean-field approximation (e.g., a single Hartree product solution), or by solving a zeroth-order Hamiltonian [70, 55]. The initial guess is then encoded on the quantum computer using a unitary transformation $U^{init}$, such that $U^{init}|0^{mT}\rangle = |\tilde{q}^{init}\rangle$.

(iii-2) Parameterization of the trial wavefunction. The second step is parameterizing the trial wavefunction using a unitary ansatz $U^{VQE}(\tau)$, where $\tau = \{\tau_\alpha\}_{\alpha=1}^{N_{VQE}}$ are the $N_{VQE}$ parameters to be optimized. Ideally, the ansatz would explore the entire Hilbert space, for example employing operations from the Pauli group. However, the exponential scaling of the number of elements with the number of qubits makes this computationally infeasible, as the optimization is performed on a classical computer. Instead, a structured, parameterized ansatz, named the VQE ansatz, is employed. Common choices include Unitary Coupled Cluster [16], Hardware-Efficient Ansatz [71], and Quantum Coupled Cluster [72]. The ansatz is typically expressed as a product of parameterized unitary operations, $U^{VQE}(\tau) = \prod_{\alpha=1}^{N_{VQE}} U_{\alpha}^{VQE}(\tau_{\alpha})$, and the trial state is now expressed as $|\tilde{q}\rangle = (\prod_{\alpha=1}^{N_{VQE}} U_{\alpha}^{VQE}(\tau_{\alpha})) U^{init}|0^{mT}\rangle$.

(iii-3) Optimization employing the classical computer. Once the ansatz is defined, the parameters $\tau$ are optimized to minimize the Rayleigh quotient $R(\tilde{q})$, which corresponds to the expectation value of $\hat{L}_{FKP}$. The optimization process is carried out on a classical computer using gradient-free methods, such as the Nelder-Mead simplex algorithm [73],
<figure><img src="image_5.png" /><figcaption>Figure 4: Schematic view of a Variational Quantum Eigensolver (VQE) hybrid quantum classical circuit. The slash ‘/’ indicates a multi-qubit line. The quantum state is constructed on the quantum circuit in two steps: i) the initial guess using $U^{\text{init}}$, and ii) the optimized unitary transformation $U^{\text{VQE}}$, which depend on a set of parameters $\tau$. The classical computer drives the measurements, with unitaries $\{U_{\alpha}^m\}$ to reconstruct the expectation value of the Rayleigh quotient and its gradients with respect to the parameters $\tau$. The classical computer then operates a feedback on the quantum computer to optimize the trial state by changing $U^{\text{VQE}}$.</figcaption></figure>

or gradient-based methods like gradient descent. Gradients can be computed either analytically or through finite-difference approximations (see, for instance, [74, 75, 76]). These methods require additional quantum measurements to compute extra points or derivatives, which increases the computational overhead.

(iii-4) *Computing other eigenstates.* Generating other eigenstates $q_k(y)$ with $k > 0$, the “excited states”, is possible within the VQE framework. Assuming that $q_0(y)$ has already been obtained, the algorithm can be modified to target excited states by either adding a penalty or imposing a constraint on the symmetry of the excited state, when applicable [77], or by enforcing the orthogonality requirement $\langle \tilde{q} | q_0 \rangle = 0$ [78].

The number of generators for the target unitary transformations is roughly $2^{m_T}(2^{m_T-1}-1)$, implying an exploration of a space scaling as $O((2^{2\tilde{m}})^v)$, and thus exponential complexity with the number of dimensions $v$. Nevertheless, the Solovay-Kitaev theorem ensures that any unitary can be approximated to accuracy $\epsilon$ using a number of gates scaling as $O(\log^c(1/\epsilon))$ where $1 \le c < 4$ [79]. This allows circuit depth - the maximum number of sequential gates - to grow sub-exponentially, making low-depth approximations feasible. Most VQE ansatzes exploit this to avoid exponential scaling while maintaining reasonable accuracy. Circuit depth determines both computation time and cumulative noise, and depends on the number of qubits and their entanglement structure. In practice, depths range from a few dozen to dozens of thousands, depending on the ansatz and problem complexity [63, 80, 81]. In our case, this depth will depend on the number $v$ of dimensions, the truncation of the basis $\sim \tilde{m}$, and the order of $\hat{L}_{\text{FKP}}$'s expansion, described here by $\mu$. We can obtain a quick estimate of the circuit depth by fixing $v = 7$, $\tilde{m} = 8$, and $\mu = 4$, with the Unitary Coupled Cluster ansatz. This estimate is comparable to the calculation for the $Cr_2$ molecule using the 1-UpCCGSD method, for which the depth is reported to be approximately $\sim 27,000$ in Ref. [63]. The current goal is to present the general VQE approach without focusing on optimization, and more quantitative estimates based on realistic simulations will be provided in future work.

## 8.2. Overlap extraction

Assuming that eigenstates are obtained, it becomes possible to extract overlaps of the form $\langle \eta | q_n \rangle$. To achieve this, we first express the overlap as $\langle \eta^\parallel | q_n \rangle$, where $|\eta^\parallel\rangle = \sum_\alpha \varphi_\alpha(\eta)|\varphi_\alpha\rangle$. Note that $\varphi_\alpha$ is the tensor-product basis; consequently, this expression cannot be computed directly. This is why we rewrite it as a manageable single product. Then, the objectives are twofold: i) construct $|q_n\rangle$ and $|\eta^\parallel\rangle$ on the quantum computer, and ii) extract the overlap as a measurement process.

The construction of $|q_n\rangle$ depicted in Sec. 8.1. Regarding the construction of $|\eta^\parallel\rangle$, we rewrite this state as a direct product:

$$|\eta^\parallel\rangle = N(\eta) \prod_{k=1}^{v} \left( \sum_{\alpha=0}^{m_k} C_{\alpha}^{(k)}(\eta_k) |\varphi_{\alpha}^{(k)}\rangle \right), \quad (8.1)$$
<figure><img src="image_6.png" /><figcaption>Figure 5: Circuit for the construction of one factor in (8.1), i.e. $\sum_{\alpha=1}^{m_k} \varphi_{\alpha}^{(k)}(\eta_k) |\varphi_{\alpha}^{(k)}\rangle$. The lines connecting circles at the nodes indicates an anti-control (control on the zero value) on the corresponding qubits.</figcaption></figure>

where $N(\eta)$ and $C_{\alpha}^{(k)}(\eta_k)$ are written as

$$N(\eta) = \prod_{k=1}^{\nu} \left( \sum_{\alpha=0}^{m_k} |\varphi_{\alpha}^{(k)}(\eta_k)|^2 \right) , \quad C_{\alpha}^{(k)}(\eta_k) = \frac{\varphi_{\alpha}^{(k)}(\eta_k)}{\sum_{\alpha=0}^{m_k} |\varphi_{\alpha}^{(k)}(\eta_k)|^2} . \qquad (8.2)$$

Note that substituting Eq. (8.2) into (8.1) effectively yields $|\eta^{\parallel}\rangle = \prod_{k=1}^{\nu} (\sum_{\alpha=0}^{m_k} \varphi_{\alpha}^{(k)}(\eta_k) |\varphi_{\alpha}^{(k)}\rangle)$. Each factor in the product on the right-hand-side of Eq. (8.1) is constructed on the quantum computer by using a product of unitaries. We give here a simple algorithm for their constructions. Starting with all qubits associated with $\eta_k$ in their ground state, the vacuum state, we apply successive rotations to distribute the coefficients $\{C_{\alpha}^{(k)}\}_{\alpha=0}^{m_k}$ on the qubit states that encode $\{|\varphi_{\alpha}^{(k)}\rangle\}_{\alpha=0}^{m_k}$. These rotations are defined by

$$\hat{R}_{k,\alpha} = \exp(-i\theta_{k,\alpha} \hat{Y}_{k,\alpha}) , \quad \theta_{k,\alpha} = \operatorname{atan}\left(C_{\alpha}^{(k)}\left[1 - \sum_{\beta=1}^{\alpha} (C_{\beta}^{(k)})^2\right]^{-1/2}\right) . \qquad (8.3)$$

As expected, the first rotation applied on the vacuum gives a superposition of two terms: $[1-(C_0^{(k)})^2]^{1/2}|0^{m_k}\rangle+C_0^{(k)}|\varphi_0^{(k)}\rangle$, the second term being the one of interest. If we apply the second rotation on this last obtained superposition, it would destroy $C_0^{(k)}|\varphi_0^{(k)}\rangle$. Instead, we must apply the second rotation only on the term containing $|0^{m_k}\rangle$. This is achieved by using an “anti-controlled” rotation,

$$\bar{C}\hat{R}_{k,1;0} = (\hat{I}_{k,0} + \hat{Z}_{k,0}) \otimes \hat{R}_{k,1}/2 + (\hat{I}_{k,0} - \hat{Z}_{k,0})/2,$$

which applies the rotation $\hat{R}_{k,1;0}$ only if $s_{k,0} = 0$. Indeed, with the controlled rotation, we now obtain

$$\bar{C}\hat{R}_{k,1;0} \left[ \sqrt{1 - (C_0^{(k)})^2} |0^{m_k}\rangle + C_0^{(k)} |\varphi_0^{(k)}\rangle \right] = \sqrt{1 - (C_0^{(k)})^2} \hat{R}_{k,1} |0^{m_k}\rangle + C_0^{(k)} |\varphi_0^{(k)}\rangle \\ = \sqrt{1 - (C_0^{(k)})^2 - (C_1^{(k)})^2} |0^{m_k}\rangle + C_1^{(k)} |\varphi_1^{(k)}\rangle + C_0^{(k)} |\varphi_0^{(k)}\rangle . \qquad (8.4)$$

To generalize this scheme, we need to define anti-controlled rotations of the form $\bar{C}\hat{R}_{k,\alpha;0,1,\dots,\alpha-1}$, where the anti-control test is applied on all preceeding qubits, $\{|s_{k,\beta}\rangle\}_{\beta=0}^{\alpha-1}$. These anti-controlled rotations are then used in place of the rotations $\hat{R}_{k,\alpha}$. The circuits that achieve this construction of a factor is depicted in Fig. 5. Equipped with this circuit, the unitary transformation, $U_{\eta}$, that builds $|\eta^{\parallel}\rangle/N(\eta)$ is given by

$$U_{\eta} = \prod_{k=1}^{\nu} \prod_{\alpha=0}^{m_k} \bar{C} \hat{R}_{k,\alpha;0,1,\dots,\alpha-1}.$$

To extract the overlap $\langle\eta^{\parallel}|q_n\rangle$, we propose using the Hadamard test [60]. In combination with the qubits used to construct $|q_n\rangle$, the Hadamard test employs an additional “ancilla” qubit. Using the Hadamard gate and controlled-unitaries, this ancilla qubit is used to construct a superposition of $|q_n\rangle$ and $|\eta^{\parallel}\rangle$ in the entangled form (see Fig. 6),

$$U_{\eta}^{\dagger}((|\eta^{\parallel}\rangle/N(\eta) + |q_n\rangle) \otimes |0\rangle + (|\eta^{\parallel}\rangle/N(\eta) - |q_n\rangle) \otimes |1\rangle)/2.$$
<figure><img src="image_7.png" /><figcaption>Figure 6: Circuit for Hadamard test employed in measuring the overlap $\langle 0^{m_T} | U_\eta^\dagger U_{\text{VQE}}(\tau) U_{\text{init}} | 0^{m_T} \rangle = \langle \eta^\|$|q_n\rangle$. The slash '/' indicates a multi-qubit line. The lines connecting black point at the node indicates a control (control on the value equal to 1) on the corresponding qubits.</figcaption></figure>

By measuring $\hat{Z}$ for the ancilla qubit, the result will be -1 with probability $p_− = (1 + \Re\langle\eta^\|$|q_n\rangle/N(η))/2$, or +1 with probability $p_+ = (1 - \Re\langle\eta^\|$|q_n\rangle/N(η))/2$, and where $\Re$ denotes the real part. After reconstructing the probabilities $p_±$, the overlap is extracted as $\langle\eta^\|$|q_n\rangle = (p_− − p_+) \times N(η)$ (accounting for the fact that $\langle\eta^\|$|q_n\rangle ∈ ℝ$). The minimal number of qubits required is $m_T + 1$, accounting for the additional ancilla qubit used in the Hadamard test. With this minimal setup, number of operations used for the anti-controlled rotations $\tilde{C}\hat{R}_{k,α;0,1,\dots,α−1}$ is quadratic with the number of controlled qubits [60, 61]. Consequently, the number of layers needed to construct the state $|\eta^\|$\rangle/N(η) scales cubically with $m_T$. An alternative approach to reduce the number of layers is to introduce additional working qubits. In particular, if $m_T$ working qubits are added, the number of layers can be reduced to scale linearly with $m_T$ [60, 61]. Moreover, by the Solovay-Kitaev theorem, the depth of the rotations $\hat{R}_{k,α}$ scales as $O(\log^c(1/ε))$ [60], for a target error $ε$. As a result, the overall circuit depth for the overlap extraction scales linearly with the number of qubits, or equivalently, with the number $v$ of dimensions. However, the number of qubits is now increased to $2m_T + 1$.

## 9. Conclusions

In this paper, we have developed a methodology and demonstrated a quantum computing framework to address a key challenge in probabilistic learning on manifolds (PLoM). This challenge involves solving the spectral problem of the high-dimensional FKP operator, a task that classical computational methods cannot efficiently handle. Notably, this problem may also be relevant to frameworks and applications beyond those that motivated the developments presented in this work.

In addition to the methodological aspects, we introduced a polynomial chaos expansion in the Gaussian Sobolev space for the potential in the Schrödinger operator associated with the FKP operator. This approach preserves the algebraic properties of the potential and enables second quantization through the introduction of creation and annihilation operators. By reformulating the eigenvalue problem in Fock space, we derived explicit formulas for both the Laplacian and the potential, facilitating the implementation of the FKP operator as a Hamiltonian in quantum circuits.

To provide a comprehensive perspective, we also included a brief overview of foundational concepts related to the construction of quantum states and the measurement of their observables, particularly for readers less familiar with quantum computing. To illustrate the essential coupling between a quantum computer and a classical computer, we chose to present the Variational Quantum Eigensolver (VQE) approach instead of the Quantum Phase Estimation (QPE) algorithm, given the latter's greater complexity and non-trivial nature. Furthermore, QPE requires longer quantum circuits, with a larger the number of qubits than what is currently available in the NISQ era [51] and the ability for error correction [67, 82].

The next step in this work will be the implementation of the proposed method on a quantum computer.

## References

[1] C. Soize, R. Ghanem, Transient anisotropic kernel for probabilistic learning on manifolds, Computer Methods in Applied Mechanics and Engineering 432 (2024) 117453. doi:10.1016/j.cma.2024.117453.

[2] B. Liu, M. Ortiz, F. Cirak, Towards quantum computational mechanics, Computer Methods in Applied Mechanics and Engineering 432 (2024) 117403. doi:10.1016/j.cma.2024.117403.

[3] A. Wulff, B. Chen, M. Steinberg, Y. Tang, M. Möller, S. Feld, Quantum computing and tensor networks for laminate design: a novel approach to stacking sequence retrieval, Computer Methods in Applied Mechanics and Engineering 432 (2024) 117380. doi:10.1016/j.cma.2024.117380.

[4] Y. Xu, J. Yang, Z. Kuang, Q. Huang, W. Huang, H. Hu, Quantum computing enhanced distance-minimizing data-driven computational mechanics, Computer Methods in Applied Mechanics and Engineering 419 (2024) 116675. doi:10.1016/j.cma.2023.116675.
[5] H. Ma, M. Govoni, G. Galli, Quantum simulations of materials on near-term quantum computers, npj Computational Materials 6 (1) (2020) 85. doi:10.1038/s41524-020-00353-z.

[6] C. Micheletti, P. Hauke, P. Faccioli, Polymer physics by quantum computing, Physical Review Letters 127 (8) (2021) 080501. doi:10.1103/PhysRevLett.127.080501.

[7] K. Mizuta, M. Fujii, S. Fujii, K. Ichikawa, Y. Imamura, Y. Okuno, Y. O. Nakagawa, Deep variational quantum eigensolver for excited states and its application to quantum chemistry calculation of periodic materials, Physical Review Research 3 (4) (2021) 043121. doi:10.1103/PhysRevResearch.3.043121.

[8] O. M. Raisuddin, S. De, Feqa: Finite element computations on quantum annealers, Computer Methods in Applied Mechanics and Engineering 395 (2022) 115014. doi:10.1016/j.cma.2022.115014.

[9] Y. Wang, J. E. Kim, K. Suresh, Opportunities and challenges of quantum computing for engineering optimization, Journal of Computing and Information Science in Engineering 23 (6) (2023) 060817. doi:10.1115/1.4062969.

[10] Z. Ye, X. Qian, W. Pan, Quantum topology optimization via quantum annealing, IEEE Transactions on Quantum Engineering 4 (2023) 1–15. doi:10.1109/TQE.2023.3266410.

[11] H. Risken, The Fokker-Planck Equation, Second Edition, Springer Verlag, Berlin, Heidelberg, 1989.

[12] C. W. Gardiner, Handbook of Stochastic Methods, Second Edition, Springer Verlag, Berlin, Heidelberg, 1985.

[13] S. McArdle, A. Mayorov, X. Shan, S. Benjamin, X. Yuan, Digital quantum simulation of molecular vibrations, Chem. Sci. 10 (22) (2019) 5725–5735. doi:10.1039/C9SC01313J.

[14] N. P. D. Sawaya, J. Huh, Quantum Algorithm for Calculating Molecular Vibronic Spectra, J. Phys. Chem. Lett. 10 (13) (2019) 3586–3591. doi:10.1021/acs.jpclett.9b01117.

[15] Ó. Amaro, L. I. I. Gamiz, M. Vranic, Variational Quantum Simulation of the Fokker-Planck Equation applied to Quantum Radiation Reaction, arXiv (Nov. 2024). arXiv:2411.17517, doi:10.48550/arXiv.2411.17517.

[16] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, J. L. O'Brien, A variational eigenvalue solver on a photonic quantum processor, Nature Communications 5 (1) (2014) 4213. doi:10.1038/ncomms5213.

[17] A. Yu. Kitaev, Quantum measurements and the Abelian Stabilizer Problem, arXiv (Nov. 1995). doi:10.48550/arXiv.quant-ph/9511026. arXiv:quant-ph/9511026.

[18] P. Malliavin, Gaussian Sobolev spaces and stochastic calculus of variations, in: Integration and Probability, Vol. 157 of Graduate Texts in Mathematics, Springer, New York, NY, 1995, pp. 229–252. doi:10.1007/978-1-4612-4202-4-5.

[19] C. Soize, R. Ghanem, Data-driven probability concentration and sampling on manifold, Journal of Computational Physics 321 (2016) 242–258. doi:10.1016/j.jcp.2016.05.044.

[20] C. Soize, R. Ghanem, Probabilistic learning on manifolds, Foundations of Data Science 2 (3) (2020) 279–307. doi:10.3934/fods.2020013.

[21] C. Soize, R. Ghanem, Physics-constrained non-Gaussian probabilistic learning on manifolds, International Journal for Numerical Methods in Engineering 121 (1) (2020) 110–145. doi:10.1002/nme.6202.

[22] C. Soize, R. Ghanem, Probabilistic learning on manifolds constrained by nonlinear partial differential equations for small datasets, Computer Methods in Applied Mechanics and Engineering 380 (2021) 113777. doi:10.1016/j.cma.2021.113777.

[23] C. Soize, Probabilistic learning inference of boundary value problem with uncertainties based on Kullback-Leibler divergence under implicit constraints, Computer Methods in Applied Mechanics and Engineering 395 (2022) 115078. doi:10.1016/j.cma.2022.115078.

[24] C. Soize, R. Ghanem, Probabilistic-learning-based stochastic surrogate model from small incomplete datasets for nonlinear dynamical systems, Computer Methods in Applied Mechanics and Engineering 418 (2024) 116498. doi:10.1016/j.cma.2023.116498.

[25] C. Soize, R. Ghanem, Probabilistic learning on manifolds (PLoM) with partition, International Journal for Numerical Methods in Engineering 123 (1) (2022) 268–290. doi:10.1002/nme.6856.

[26] C. Soize, Probabilistic learning constrained by realizations using a weak formulation of fourier transform of probability measures, Computational Statistics 38 (4) (2023) 1879–1925. doi:10.1007/s00180-022-01300-w.

[27] R. Coifman, S. Lafon, A. Lee, M. Maggioni, B. Nadler, F. Warner, S. Zucker, Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps, PNAS 102 (21) (2005) 7426–7431. doi:10.1073/pnas.0500334102.

[28] R. Coifman, S. Lafon, Diffusion maps, Applied and Computational Harmonic Analysis 21 (1) (2006) 5–30. doi:10.1016/j.acha.2006.04.006.

[29] C. Soize, Polynomial chaos expansion of a multimodal random vector, SIAM-ASA Journal on Uncertainty Quantification 3 (1) (2015) 34–60. doi:10.1137/140968495.

[30] A. Bowman, A. Azzalini, Applied Smoothing Techniques for Data Analysis: The Kernel Approach With S-Plus Illustrations, Vol. 18, Oxford University Press, Oxford: Clarendon Press, New York, 1997. doi:10.1007/s001800000033.

[31] T. Duong, A. Cowling, I. Koch, M. Wand, Feature significance for multivariate kernel density estimation, Computational Statistics & Data Analysis 52 (9) (2008) 4225–4242. doi:10.1016/j.csda.2008.02.035.

[32] J. E. Gentle, Computational statistics, Springer, New York, 2009. doi:10.1007/978-0-387-98144-4.

[33] G. Givens, J. Hoeting, Computational Statistics, 2nd Edition, John Wiley and Sons, Hoboken, New Jersey, 2013.

[34] J. L. Doob, Stochastic processes, John Wiley & Sons, New York, 1953.

[35] I. I. Guikhman, A. Skorokhod, Introduction à la Théorie des Processus Aléatoires, Edition Mir, 1980.

[36] A. Friedman, Stochastic Differential Equations and Applications, Dover Publications, Inc., Mineola, New York, 2006.

[37] C. Soize, The Fokker-Planck Equation for Stochastic Dynamical Systems and its Explicit Steady State Solutions, Vol. Series on Advances in Mathematics for Applied Sciences: Vol 17, World Scientific, Singapore, 1994. doi:10.1142/2347.

[38] C. Soize, R. Ghanem, C. Safta, X. Huan, Z. P. Vane, J. C. Oefelein, G. Lacaze, H. N. Najm, Q. Tang, X. Chen, Entropy-based closure for probabilistic learning on manifolds, Journal of Computational Physics 388 (2019) 528–533. doi:10.1016/j.jcp.2018.12.029.

[39] N. Wiener, The homogeneous chaos, American Journal of Mathematics 60 (1) (1938) 897–936.

[40] R. H. Cameron, W. T. Martin, The orthogonal development of non-linear functionals in series of Fourier-Hermite functionals, Annals of Mathematics 48 (2) (1947) 385–392. doi:10.2307/1969178.

[41] P. Krée, C. Soize, Mathematics of Random Phenomena, Reidel Pub. Co, 1986, (first published by Bordas in 1983 and also published by Springer Science & Business Media in 2012).
[42] R. J. Serfling, Approximation theorems of mathematical statistics, Vol. 162, John Wiley & Sons, 1980.

[43] C. Soize, Uncertainty Quantification. An Accelerated Course with Advanced Applications in Computational Engineering, Springer, New York, 2017. doi:10.1007/978-3-319-54339-0.

[44] D. Gottlieb, S. A. Orszag, Numerical Analysis of Spectral Methods, Society for Industrial and Applied Mathematics, 1977. arXiv:https://epubs.siam.org/doi/pdf/10.1137/1.9781611970425, doi:10.1137/1.9781611970425. URL https://epubs.siam.org/doi/abs/10.1137/1.9781611970425

[45] J. C. Light, T. Carrington, Discrete-Variable Representations and their Utilization, in: Advances in Chemical Physics, John Wiley & Sons, Inc., Hoboken, NJ, USA, 2000, pp. 263–310. doi:10.1002/9780470141731.ch4.

[46] D. J. Tannor, Introduction to Quantum Mechanics: A Time-Dependent Perspective, Harvard, 2007. URL https://ui.adsabs.harvard.edu/abs/2007iqm..book....T/abstract

[47] F. Gatti, B. Lasorne, H.-D. Meyer, A. Nauts, Applications of quantum dynamics in chemistry, Vol. 98, Springer, 2017.

[48] H.-D. Meyer, F. Gatti, G. A. Worth, Multidimensional quantum dynamics: MCTDH theory and applications, John Wiley & Sons, 2009.

[49] U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Ann. Phys. 326 (1) (2011) 96–192. doi:10.1016/j.aop.2010.09.012.

[50] M. Bachmayr, Low-rank tensor methods for partial differential equations, Acta Numer. 32 (2023) 1–121.

[51] D. Castelvecchi, Ibm releases first-ever 1,000-qubit quantum chip, Nature 624 (2023) 238. doi:10.1038/d41586-023-03854-1.

[52] B. Fauseweh, Quantum many-body simulations on digital quantum computers: State-of-the-art and future challenges, Nat. Commun. 15 (2123) (2024) 1–13. doi:10.1038/s41467-024-46402-9.

[53] Y. Wang, L. M. Sager-Smith, D. A. Mazziotti, Quantum simulation of bosons with the contracted quantum eigensolver, New J. Phys. 25 (10) (2023) 103005. doi:10.1088/1367-2630/acf9c3.

[54] P. R. Surján, Second Quantized Approach to Quantum Chemistry, Springer, Berlin, Germany, 1989.

[55] J. M. Bowman, Vibrational Dynamics of Molecules, World Scientific Publishing Company, Singapore, 2021. doi:10.1142/12305.

[56] R. Somma, G. Ortiz, J. E. Gubernatis, E. Knill, R. Laflamme, Simulating physical phenomena by quantum networks, Phys. Rev. A 65 (4) (2002) 042323. doi:10.1103/PhysRevA.65.042323.

[57] L. Veis, J. Višņák, H. Nishizawa, H. Nakai, J. Pittner, Quantum chemistry beyond Born–Oppenheimer approximation on a quantum computer: A simulated phase estimation study, Int. J. Quantum Chem. 116 (18) (2016) 1328–1336. doi:10.1002/qua.25176.

[58] E. Crane, K. C. Smith, T. Tomesh, A. Eickbusch, J. M. Martyn, S. Kühn, L. Funcke, M. A. DeMarco, I. L. Chuang, N. Wiebe, A. Schuckert, S. M. Girvin, Hybrid Oscillator-Qubit Quantum Processors: Simulating Fermions, Bosons, and Gauge Fields, arXiv (Sep. 2024). arXiv:2409.03747, doi:10.48550/arXiv.2409.03747.

[59] S. Malpathak, S. D. Kallullathil, A. F. Izmaylov, Simulating Vibrational Dynamics on Bosonic Quantum Devices, arXiv (Nov. 2024). arXiv:2411.17950, doi:10.48550/arXiv.2411.17950.

[60] M. A. Nielsen, I. L. Chuang, Quantum computation and quantum information, Cambridge university press, 2010.

[61] A. Barenco, C. H. Bennett, R. Cleve, D. P. DiVincenzo, N. Margolus, P. Shor, T. Sleator, J. A. Smolin, H. Weinfurter, Elementary gates for quantum computation, Phys. Rev. A 52 (5) (1995) 3457–3467. doi:10.1103/PhysRevA.52.3457.

[62] I. Djordjevic, Chapter 3 - Quantum Circuits and Quantum Information Processing Fundamentals, in: Quantum Information Processing and Quantum Error Correction, Academic Press, Cambridge, MA, USA, 2012, pp. 91–117. doi:10.1016/B978-0-12-385491-9.00003-4.

[63] J. Tilly, H. Chen, S. Cao, D. Picozzi, K. Setia, Y. Li, E. Grant, L. Wossnig, I. Rungger, G. H. Booth, J. Tennyson, The Variational Quantum Eigensolver: A review of methods and best practices, Phys. Rep. 986 (2022) 1–128. doi:10.1016/j.physrep.2022.08.003.

[64] V. Verteletskyi, T.-C. Yen, A. F. Izmaylov, Measurement optimization in the variational quantum eigensolver using a minimum clique cover, The Journal of chemical physics 152 (12) (2020). doi:10.1063/1.5141458.

[65] T.-C. Yen, A. F. Izmaylov, Cartan Subalgebra Approach to Efficient Measurements of Quantum Observables, PRX Quantum 2 (4) (2021) 040320. doi:10.1103/PRXQuantum.2.040320.

[66] T.-C. Yen, A. Ganeshram, A. F. Izmaylov, Deterministic improvements of quantum measurements with grouping of compatible operators, non-local transformations, and covariance estimates, npj Quantum Inf. 9 (14) (2023) 1–7. doi:10.1038/s41534-023-00683-y.

[67] M. Reiher, N. Wiebe, K. M. Svore, D. Wecker, M. Troyer, Elucidating reaction mechanisms on quantum computers, Proc. Natl. Acad. Sci. U.S.A. 114 (29) (2017) 7555–7560. doi:10.1073/pnas.1619152114.

[68] J. Preskill, Quantum Computing in the NISQ era and beyond, Quantum 2 (18) 79. arXiv:1801.00862v3, doi:10.22331/q-2018-08-06-79.

[69] J. R. McClean, J. Romero, R. Babbush, A. Aspuru-Guzik, The theory of variational hybrid quantum-classical algorithms, New Journal of Physics 18 (2) (2016) 023023. doi:10.1088/1367-2630/18/2/023023.

[70] E. B. Wilson, Jr., J. C. Decius, P. C. Cross, Molecular Vibrations: The Theory of Infrared and Raman Vibrational Spectra (Dover Books on Chemistry), Dover Publications, Mineola, NY, USA, 1980. URL https://www.amazon.com/Molecular-Vibrations-Infrared-Vibrational-Chemistry/dp/048663941X

[71] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, J. M. Gambetta, Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets, Nature 549 (2017) 242–246. doi:10.1038/nature23879.

[72] I. G. Ryabinkin, T.-C. Yen, S. N. Genin, A. F. Izmaylov, Qubit coupled cluster method: a systematic approach to quantum chemistry on a quantum computer, Journal of chemical theory and computation 14 (12) (2018) 6317–6326. doi:10.1021/acs.jctc.8b00932.

[73] I. Miháliková, M. Friák, M. Pivoluska, M. Plesch, M. Saip, M. Šob, Best-Practice Aspects of Quantum-Computer Calculations: A Case Study of the Hydrogen Molecule, Molecules 27 (3) (2022) 597. doi:10.3390/molecules27030597.

[74] J. Romero, R. Babbush, J. R. McClean, C. Hempel, P. J. Love, A. Aspuru-Guzik, Strategies for quantum computing molecular energies using the unitary coupled cluster ansatz, Quantum Sci. Technol. 4 (1) (2018) 014008. doi:10.1088/2058-9565/aad3e4.

[75] L. Wang, J. Zhu, J. C. Spall, Mixed simultaneous perturbation stochastic approximation for gradient-free optimization with noisy measurements, in: 2018 Annual American Control Conference (ACC), IEEE, 2018, pp. 27–29. doi:10.23919/ACC.2018.8430974.

[76] X. Bonet-Monroig, H. Wang, D. Vermetten, B. Senjean, C. Moussa, T. Bäck, V. Dunjko, T. E. O’Brien, Performance comparison of optimization methods on variational quantum algorithms, Phys. Rev. A 107 (3) (2023) 032407. doi:10.1103/PhysRevA.107.032407.

[77] I. G. Ryabinkin, S. N. Genin, A. F. Izmaylov, Constrained Variational Quantum Eigensolver: Quantum Computer Search Engine in the Fock
Space, J. Chem. Theory Comput. 15 (1) (2019) 249–255. doi:10.1021/acs.jctc.8b00943.

[78] O. Higgott, D. Wang, S. Brierley, Variational Quantum Computation of Excited States, Quantum 3 (2019) 156. arXiv:1805.08138v5, doi:10.22331/q-2019-07-01-156.

[79] C. M. Dawson, M. A. Nielsen, The Solovay-Kitaev algorithm, Quantum Inf. Comput. 6 (1) (2006) 81–95. doi:10.5555/2011679.2011685.

[80] D. A. Fedorov, B. Peng, N. Govind, Y. Alexeev, VQE method: a short survey and recent developments, Mater. Theory 6 (1) (2022) 1–21. doi:10.1186/s41313-021-00032-6.

[81] Y. Zhang, L. Cincio, C. F. Negre, P. Czarnik, P. J. Coles, P. M. Anisimov, S. M. Mniszewski, S. Tretiak, P. A. Dub, Variational quantum eigensolver with reduced circuit complexity, npj Quantum Information 8 (1) (2022) 96. doi:10.1038/s41534-022-00599-z.

[82] L. Postler, S. Heußen, I. Pogorelov, M. Rispler, T. Feldker, M. Meth, C. D. Marciniak, R. Stricker, M. Ringbauer, R. Blatt, P. Schindler, M. Müller, T. Monz, Demonstration of fault-tolerant universal quantum gate operations, Nature 605 (2022) 675–680. doi:10.1038/s41586-022-04721-1.