Quantum computing for classical problems: Variational Quantum Eigensolver for activated processes

The theory of stochastic processes impacts both physical and social sciences. At the molecular scale, stochastic dynamics is ubiquitous because of thermal fluctuations. The Fokker-Plank-Smoluchowski equation models the time evolution of the probability density of selected degrees of freedom in the diffusive regime and it is therefore a workhorse of physical chemistry. In this paper we report the development and implementation of a Variational Quantum Eigensolver procedure to solve the Fokker-Planck-Smoluchowski eigenvalue problem. We show that such an algorithm, typically adopted to address quantum chemistry problems, can be applied effectively to classical systems paving the way to new applications of quantum computers. We compute the conformational transition rate in a linear chain of rotors experiencing nearest-neighbour interaction. We provide a method to encode on the quantum computer the probability distribution for a given conformation of the chain and assess its scalability in terms of operations. Performance analysis on noisy quantum emulators and quantum devices (IBMQ Santiago) is provided for a small chain showing results in good agreement with the classical benchmark without further addition of any error mitigation technique.


Introduction
Since the dawn of quantum mechanics the study of quantum phenomena has been one of the key focus of the scientific community. For many years quantum systems have been investigated with the sole purpose of deepening our understanding on their microscopic behaviour. In recent years, however, scientific and technological advancement opened the way to the direct manipulation and control of quantum systems with the aim of solving practical problems. This shift in paradigm, sometime referred as second quantum revolution [1], changed the way of looking at the quantum world moving it from an object of study to the status of a powerful new tool. Quantum computation is surely one of the key examples of these new quantum technologies and unearthing computational tasks for which the use of quantum resources can offer a significant increase in efficiency with respect to the best classical algorithmic counterpart is one of the key challenges of the field. Many efforts of the scientific community have been devoted towards such an alternative paradigm of computation leading to the development of several algorithms that heralded its disruptive potential [2]. Particularly, with regard to relevant breakthroughs in physics and chemistry, quantum simulation is considered among the first near-term applications of this field [3]. In the present era of Noisy Intermediate Quantum Devices (NISQ) [4], hybrid quantum algorithms represent a valuable tool to boost classical computational power using small quantum processors. In this framework, some of the more performing algorithms are the so called Variational Hybrid Algorithms (VHA) where the quantum device is tasked with the preparation of a parametrized trial state and the computation of a cost function which, in turn, is minimized by an external classical optimization routine [4]. This method requires the quantum device to run circuits of reasonably shortdepth that cope with the coherence time and error rate of modern NISQ devices [5].
The literature about VHA is constantly growing both in terms of different implementations and in terms of their applications. These range from the study of the electronic structure of small molecules [6,7] and molecular vibrations [8] to the simulation of condensed matter and high-energy physics [9,10]. In contrast, less attention has been dedicated to the application of quantum computing in the context of stochastic simulations (i.e., simulations of systems characterized by variables with a random character, and thus by classical probability distributions). We believe that leveraging standard methods of stochastic calculus, while harnessing the potential of quantum computers, may lead to the development of new quantum tools in the field of complex system simulations. This can possibly pave the way to important advances in chemistry, biochemistry and all the other connected fields in which stochastic analysis plays a relevant role such as quantitative finance, epidemiology, and computational fluid dynamics [11][12][13]. Here we propose a new VHA to solve the Fokker-Planck-Smoluchowski (FPS) eigenvalue problem based on the isomorphism existing between the FPS operator and the quantum Hamiltonian [14,15]. As an example, in this paper we will focus on the approximation of the kinetic rate constant for the conformational transition in linear chain molecules that is connected with the first non-vanishing eigenvalue of the FPS operator.
Despite of its simplicity, this model has been used as theoretical starting point for the conformational analysis in polymers [16][17][18] whose study has deep implications in many different areas [19].
The present contribution is organized as follows: Section 2 outlines the theoretical basis of our approach. We give a general overview on how the FPS eigenvalue problem can be treated on a quantum computer and then we formalize the stochastic description of the conformational dynamics of a chain of rotors characterized by a bistable intramolecular potential energy surface.
This system is the case-study of the present work and it models the molecular dynamics of simple polymers. An important aspect of the Fokker-Planck dynamics in bistable systems is that the first excited state can become almost degenerate with the ground state especially for large activation barriers. The corresponding small, but non-zero, eigenvalue determines the rate of passage of particles from one potential well to the other, so defining a theoretical framework for a microscopic derivation of the rate coefficients of elementary two-state reactions as initiated by Kramers [20,21]. Although an accurate calculation of the first non-zero eigenvalue assumes a central relevance in the study of kinetic processes, it is computationally challenging especially for large systems characterized by many coupled coordinates. In this context, the application of new quantum computational tools can be convenient to possibly tackle previously inaccessible stochastic problems.
We propose to harness the exponential storage capacity of a quantum computer to handle this problem and we detail the implementation of the stochastic problem in the quantum architecture in Section 3. The eigenvalue problem associated to the Fokker-Planck operator is encoded into the quantum register with a simple binary mapping of the basis set in the computational basis. By exploiting the symmetry of the FPS operator, we discuss a VQE-like algorithm to obtain the first non-zero eigenvalue of the stochastic operator, corresponding to the kinetic rate of the isomerization reaction. Section 4 presents the results obtained for the simulation of small chains characterized by different activation barriers. We discuss several aspects which play a critical role on the performance of the algorithm, such as the choice of the variational ansatz, the size of the basis set which is reflected on the number of qubits to be used and the effect of the noise on the accuracy of the result. In the concluding section, we critically discuss the potential quantum advantage of the proposed algorithm. The key advantage stems from the linear scaling in terms of memory resources, indeed the binary mapping implies a linear relation between the number of rotors in the polymeric chain and the number of qubits used in the implementation. This is an exponential gain compared to classical memory resources which grow with the dimension of the vector space representing the state of the system, dimension increasing exponentially with the number of relevant degrees of freedom. On the other hand, the present implementation suffers of an exponential scaling of the number of expectation values which needs to be measured to solve the problem and a somehow limited accuracy of the results for low barriers. We therefore point out future developments which are needed to bring the stochastic molecular dynamics amongst the applications which will take great advantage by the advent of quantum computers.

Theory
The stochastic description of molecular systems is a broad research field that embraces a wide range of theoretical tools to allow the interpretation of the fluctuating dynamics of microscopic systems in a plethora of physical conditions [22]. We will focus our attention on the simple case of a stationary and Markovian diffusive process. This kind of behavior is typical of the so called overdamped regime [23], in which a time-scale separation is observed due to the configurational variables' dynamics being slower than the evolution of the conjugated momenta. By looking at the system within a time-scale comparable with the fluctuations of the configurational variables q(t), the momentum variables rapidly loose correlation, so that one can assume they are relaxed in the equilibrium distribution.
In order to describe the time evolution of the probability distribution ρ(q, t), that depends on some configurational variables q, the Fokker-Planck-Smoluchowski (FPS) equation can be invoked [23]: where the FPS operatorΓ is defined as follows: Here D(q) is the diffusion tensor [23] and ρ eq (q) represents the Boltzmann equilibrium probability distribution: where U (q) is the mean-field potential and β = (k B T ) −1 .
The FPS operator is real and positive semi-definite such that each eigenvalue λ k obeys to the constraint λ k ≥ 0 ∀k. From now on we will consider the eigenvalues λ k sorted in ascending order such that λ k ≤ λ k+1 and we will indicate with ψ k (q) the corresponding eigenfunctions.
The lowest eigenfunction ψ 0 (q), corresponding to the null eigenvalue λ 0 = 0, is provided by the equilibrium distribution ρ eq (q). If the time evolution of a generic non-equilibrium distribution ρ(q, t) = k c k (t)ψ k (q) is considered, the following relation holds: From this equation we can easily conclude that, since all the non-zero eigenvalues play the role of exponential decay rates, a generic probability distribution must relax over time toward the equilibrium state ψ 0 (q).
The goal of this work is to investigate a possible strategy to solve, with the support of a quantum computer, the eigenvalue problem associated to the FPS operator from Eq. 2. This however cannot be done directly since the FPS operator is non Hermitian and, as such, cannot be converted into a qubit Hamiltonian to be used in the variational procedure. This issue has a simple workaround that involves the introduction of a symmetrized, self-adjoint operatorΓ: which leads to a "Schrödinger-like" representation of the FPS operator [15] withΓ in place of the Hamiltonian operator.
This new symmetrized operator conserves the same eigenvalues λ k of the original FPS operator while its eigenfunctionsψ k (q) are connected with the eigenfunctions ψ k (q) ofΓ by the relation: For more details about the symmetrized FPS operator we refer to [23,24].
Here, inspired by the Variational Quantum Eigensolver (VQE) [25], we propose an hybrid algorithm to obtain the first non-vanishing eigenvalue (i.e., λ 1 ) ofΓ. The procedure exploits the quantum computer to prepare, by means of a parametrized circuit, a trial wave-function to approximate the state correspondent to theψ 1 (q) target eigenfunction. Subsequently an external classical optimization routine updates the parameters entering the variational circuit on the basis of the measured expectation value Γ . The procedure is iterated until convergence. In order to measure Γ the symmetrized FPS operator will be represented as a linear combination of Pauli strings.Γ where the explicit formulation for the weights γ j and the Pauli stringsP j ∈ {σ x , σ y , σ z , I} ⊗N depends upon the adopted mapping (see Sec. 3).
Notably, the direct application of the VQE to the FPS eigenvalue problem would lead to the trivial solution λ 0 = 0, as already mentioned. In the following (Sec. 2.1), we will report the implementation of a particular problem for which the symmetry allows to easily target the desired eigenvalue. Nevertheless, various techniques can be applied to overcome this issue in analogy with the computation of excited states of a molecular Hamiltonian [26][27][28].
In the following section we will present, on the basis of this theoretical introduction, a linear chain molecule as model system to implement this procedure. We will focus on the calculation of the first excited state that, as will be discussed, assumes a central relevance in the study of kinetic processes.

Stochastic dynamics and conformational transition rate
The FPS equation describes the time evolution of the probability distribution over all the configurational domain allowing for a continuous description of the system dynamics. This degree of detail is surely very informative, but is also not always the best way of describing a reactive process. Often, either for simplicity or due to the lack of information about the initial distribution, less detailed descriptions are adopted, in which, rather than describing the molecular configurations, one usually thinks to the process as the inter-conversion between reactant and product structures [15]. In doing so we implicitly moved our description from a continuous FPS problem to a discrete Master equation with rates identifying the kinetic constants for the interconversions between different states.
In this paper, we will consider a symmetrical double minimum system in which two stable molecular configurations are kept apart by a potential barrier. This is the prototypical example of a simple isomerization process in which the kinetic constants k of the direct and reverse processes are the same. If the Master equation description is adopted, it is trivial to show how the population relaxation rate toward equilibrium is described by the exponential decay exp (−2kt). If the barrier is large enough, the FPS operator eigenvalue spectrum shows a large gap between the first non-zero eigenvalue λ 1 and the higher ones. In this picture, the states corresponding to λ n λ 1 have the character of fast (intra-minimum) relaxing modes, while the longer-lived state, corresponding to λ 1 , represents the kinetic relaxation mode connected with the population transfer between sites ( Figure 1). The exponential form dictating the relaxation dynamic can be obtained from Eq. 4 as exp (−λ 1 t) from which the simple relation, λ 1 = 2k, can be established between the kinetic constant k and the first non-vanishing FPS operator eigenvalue λ 1 .

The rotor chain as a model for molecular conformational dynamics
Now that the theoretical framework has been outlined let us introduce a model system to be used as a test subject. Let us consider a chain of N + 1 rotors free to rotate around a common axis, see Fig.2. Let ϕ k (for k = 0 , 1 , . . . , N ) be the angular coordinate representing the orientation of the k-th rotor with respect to a laboratory axis orthgonal to the chain axis and θ k := ϕ k − ϕ k−1 the relative orientation angles, hereafter called as dihedral angles, between the k-th rotor and the previous one in the chain. Let us assume that the rotors interact on the basis of a nearest-neighbour periodic potential U k (θ k ) = U k (θ k + 2π), depending on the corresponding dihedral angle θ k , such that the overall mean-field potential for the chain can be specified as: in which θ = (θ 1 , ..., θ N ) ∈ R N represents the vector encoding the internal configuration of the chain. In order to have a clear picture of the kinetic process, let us consider a chain in which N − 1 mono-stable dihedrals are coupled to a single bi-stable potential. In what follows the mono-stable potential will be specified as: while the bi-stable one by the form: where we have adopted the symbols ∆ r and ∆ nr to indicate the barrier height encountered when moving, respectively, along the reactive and non-reactive dihedral coordinates.
Furthermore, let us assume the absence of hydrodynamic interactions [19] amongst the rotors which, therefore, are considered as characterized by independent and constant diffusion coefficients (D k for k = 0 , 1 , . . . , N ). Under these assumptions the FPS operator for the process can be written as: in which ϕ = (ϕ 0 , ..., ϕ N ) ∈ R N +1 represents the vector encoding the orientation of each rotor. The FPS operator defined in Eq. 11 contains a degenerate degree of freedom. Because of the pairwise decomposition of the mean field potential (see Eq. 8), any homogeneous rotaion of the rotors brings the system into an equivalent conformational state. Such a degeneracy can be exploited by introducing a coordinate representation based on the dihedral angles θ and the overall chain orientation Φ defined as: If the coordinate transformation ϕ → [Φ, θ] is applied, the FPS operator can be decomposed as the sum of a global orientation operatorΓ Φ , N single dihedral termŝ Γ k and N − 1 two dihedrals interaction termsΓ k,k+1 : with these operators specified as: In order to study the conformational dynamics, the term Γ Φ , depending on the global orientation angle Φ, can be neglected since it does not affect the internal degrees of freedom.
At this point, in order to solve the eigenvalue problem associated with the FPS operator and to obtain the integrals required for the VQE algorithm, a convenient basis set must be selected. For the purpose of this paper a set {φ n (θ)} of composite basis functions has been constructed according to: where Ξ n k (θ k ) are the eigenfunctions of the single dihedral operatorΓ k and n = (n 1 , ..., n N ) represents a vector containing the order of the eigenfunctions for each dihedral angle. All the isolated dihedral eigenfunctions Ξ n k (θ k ) have been obtained by expanding the proper FPS operator over a Fourier basis set. Given the parameterized forms of Eq. 9 and Eq. 10 for the contributions to the mean field potential, one derives quite easily that the FPS operatorΓ is invariant with respect to the change of sign of all the dihedral angles. In mathematical terms this condition implies the commutation relation Γ ,Ŝ = 0, whereŜ is the symmetry operator for the change of sign of the dihedral angles. In essence the FPS operator acts on a symmetry parted Hilbert space H = H + ⊕ H − and its action does not mix functions belonging to the even H + and to the odd H − sub-spaces. So, the variational theorem applies to the lowest eigenvalue of both sub-spaces. Therefore both the first eigenstateψ 0 (q) = ρ eq (q) − 1 2 and the second eigenstateψ 1 , associated with the aforementioned kinetic mode, can be obtained from the variational procedure applied to the even and odd sub-spaces, respectively. In this framework, the application of a VQE procedure to get λ 1 requires the representation of FPS operator on the basis of odd basis functions only.
Given p n k ∈ (−1, 1) as the parity index of each dihedral eigenfuntion Ξ n k (θ k ), the condition: must be verified in order to consider only functions belonging to the odd H − sub-space. The results from the classical computations of the stochastic problem will be discussed, together with the ones coming from the VQE procedure, in Sec. 4.

Implementation
The first step in the implementation is represented by the translation of the targetψ 1 (q) eigenfunction from a linear combination of N -dihedral states (Eq. 17) to a linear combination of quantum computer bitstring states. Such a procedure will provide also an expression for the symmetrized FPS operator as a combination of Pauli strings. In this paper, we opted for a binary mapping in which a progressive numerical label j is assigned to every basis function φ n (θ). Each label is then directly mapped through its binary representation bin(j) in a state |j ≡ |bin(j) of the qubit register. Given a set of J basis functions a total of Q = log 2 (J) qubits is required.
We start the mapping process by considering the FPS operator representation in the previously introduced symmetry adapted basis: Let us indicate with δ q (j) the q-th digit of the binary representation of j such that |j ≡ |δ 1 (j), ..., δ Q (j) , we can rephrase the operator representation, Eq. 19, in a qubit-wise form: By substituing in the last equation each single qubit outer product as combination of Pauli matrices, we obtain the following decomposition in a linear combination of Pauli strings like in Eq. 7: As we can see each matrix element of the FPS operator is translated into a set of 2 Q strings, each of which is composed by Q Pauli matrices. Each set of strings is specific for the given binary bitstring representation of all the possible r and c states.
Here we discuss the potential quantum advantage of the approach. Considering that the number of basis states J to encode shows a roughly exponential increase with the system size (e.g. with the length of the rotor chain), the adoption of this mapping protocol allows to overcome the exponential growth of a classical storage space by requiring a number of qubits which scales linearly with the system size. To be more precise, let us indicate with m a positive integer related to the number of single-dihedral basis states used to set up the product states in Eq. 17 and let N be the number of dihedral angles in the chain. From these assumptions we can estimate the order of J to be J ∼ O(m N ) and, consequently, only O(N log 2 (m)) qubits are needed to represent the system on the quantum register. In order to fully characterize the efficiency associated with the procedure, a more in depth analysis of the mapping is required. To do so, the measurement process has to be taken into account. As mentioned in Sec. 2, the expectation value measurement is related, according to Eq. 7, to the expectation values of the Pauli strings which, in turn, are affected both by the structure of the operator and the mapping adopted. Since we have already discussed the features of the mapping now we should turn our attention to the estimation of the non-vanishing elements of the FPS operator.
Due to the nearest-neighbour nature of the interaction theΓ operator is sparse. Indeed, it has a block structure where only states differing for single-dihedral function of adjacent rotors are connected. Nevertheless, even with this great reduction, the number of Pauli strings to be measured is exponentially growing: the number of elements connecting states that differ for two single-dihedral functions is O(m 2 JN ) to which we add the O(mJN ) elements differing for a single-dihedral state and the J diagonal elements. It may be worth to notice that this estimate does not take into account any symmetry argument that could reduce the count and does not consider additional strategies that can be employed to reduce even more the number of required measurements. Some consideration about this point will be discussed in Sec. 5.
It is important to notice that this is not the only possible choice for the mapping and different strategies could be more or less convenient depending on the system at hand; ref. [29] details several possibilities available within the range set by the binary mapping (more dense) and the unary encoding (less dense, where d-level systems are mapped into d bitstring states with all the qubits in the state |0 but one distinguishing among different states, see for instance refs. [30,31]). As a general rule, a more or less dense mapping results in an accordingly higher or smaller cost in terms of length and number of Pauli strings generated by each mapped matrix element. These different factors must be balanced according to the system and its size.
That said, the binary mapping suffices for our scope which is to give a proof of concept of the application by considering a rather small system. The design of a more efficient mapping is left for future development and discussed in Sec. 5.

Computational details
Here we provide the details about the numerical results that are shown in Sec. 4. The matrix elements of the FPS operator have been obtained with a homemade C++ code [32] which has also been used as benchmark for the implementation of the hybrid algorithm in a Python code [33].
In order to evaluate the accuracy and the performance of our method we have considered a chain of three rotors. The two resulting dihedral angles have been parametrized as follows: the first one experiences a bi-stable potential, as dictated by Eq. 10, while the second one is described by a mono-stable potential as from Eq. 9. The barrier height ∆ r for the bistable dihedral has been varied from 0.5 k B T to 3.0 k B T while the ∆ nr parameter for the single-minimum dihedral has been kept fixed at 1 k B T . The diffusion coefficients for all three rotors have been set to the same relative value of 1. The term "relative" is here employed since, due to the structure of Eq. 11, any scale factor can be applied to the FPS operator to scale the result in term of the desired diffusion coefficient.
The algorithm has been tested on Q = 2, 3, 4 qubits quantum registers using a composite basis set filling completely all the 2 Q register states. In order to generate the proper composite basis set for the system an adequate number of single dihedral basis functions had to be selected according to Eq. 17 and Eq. 18. Indicating the number n 1 of isolated basis function for the double minimum dihedral and with n 2 the one selected for the single minimum one, the following (n 1 , n 2 ) basis set have been selected for each value of Q: For the 2 qubit system a (4, 2) basis set has been employed, for the 3 qubit we opted for a (4, 4) configuration and a (8, 4) set has been used for the 4 qubit calculations. Note that, due to symmetry restrains, only half of the composite basis functions generated by each (n 1 , n 2 ) set are used in the actual variational procedure.
A RyRz heuristic variational circuit has been employed in order to prepare the trial wave-function required in the VQE procedure [34]. Two main elements form such a circuit: a layer of parameterized R y and R z rotations, used to apply an arbitrary single qubit rotation to each qubit in the quantum register, and an entangler block to introduce correlation in the trial wavefunction. A variational circuit of depth d = 1 can be obtained enclosing an entangler block between two rotation layers. Higher depth circuits can be obtained by concatenating couples of entangler blocks and rotation layers to the depth d = 1 circuit. A Q-qubits RyRz circuit of depth d is defined by 2Q(d + 1) variational parameters. Two different entangler blocks configurations have been considered in this paper: one creating entanglement between adjacent qubits, hereafter addressed with the name "linear", and one creating entanglement between all couples of qubits in the register, called "full" from now on. In Figure 3 a schematic representation of these parameterized circuits is shown.
In order to start the VQE procedure the register is set to the all zeros state |0 that has been selected in order to correspond to the basis set function φ 0 (θ 1 , θ 2 ) = Ξ 1 (θ 1 )Ξ 0 (θ 2 ) given by the product of the first excited state Ξ 1 (θ 1 ) for the isolated double minimum system and the ground-state Ξ 0 (θ 2 ) of the single minimum one. The rotation angles of the variational form are set randomly at every run.
Four different optimization algorithm, namely the downhill simplex method, Sequential least square programming (SLSQP), Constrained optimization by linear approximation (COBYLA) and Simultaneous perturbation stochastic approximation (SPSA), have been applied with a variable degree of success based on the type of simulation. The IBM's Qiskit module [35] (version 0.22) has been employed for the simulation of the quantum circuits.  The implementation has been verified using the built-in Qiskit Statevector Simulator, an ideal simulator that return the statevector of the quantum register. Then, we characterized the statistics induced by a finite number of measurements of the quantum state in a noiseless quantum simulation varying the number of runs (that are the times a circuit is executed and measured). Finally the behavior of the algorithm has been tested on a noisy simulator. For this last task, the noise model of the IBMQ Santiago has been employed and a basic readout error mitigation protocol (based on Ref. [36]) has been applied to each circuit simulation as implemented in Qiskit.
In order to perform the simulations with the IBMQ Santiago noise model, the built-in Qiskit transpilation routines have been employed to map the circuit on the quantum computer topology. During the process the gates have been translated into U 1 , U 2 , U 3 and CNOT basis gates. Given the 5 inline qubits topology of the IBMQ Santiago quantum computer the transpilation of a linear entangler RyRz variational form on Q qubits appears to be very efficient resulting in a gate count of 4Q−3 (3Q−2 U n and Q − 1 CNOT) without requiring any qubit swap. The transpilation of a full entangler RyRz circuit represent instead a more complex task requiring qubit swaps due to the presence of CNOT entagling operations between non-adjacent physical qubits.
The expectation value measurement has been carried out either by measuring independently all the Pauli strings of the mapped operator or by retrieving from Qiskit the built-in method which groups all the terms that are simultaneously diagonalizable [37]. Both methods performed equally well and returned comparable results.

Results
In this section we discuss the results obtained with numerical simulations of the algorithm in several platforms: the Qiskit Statevector Simulator, the Qiskit Quantum Assembly Simulator (QASM) with and without a noise model and the IBMQ Santiago device [35].
In Sec. 4.1 we assess the variational circuit performance by discussing the results obtained with a noise-free optimization. The importance of this question is highlighted by the large amount of research on the ansatze efficiency when this kind of algorithms are applied to quantum chemistry problems [25,38,39]. Finally, in Sec. 4.2 we show the results of noisy simulations of the algorithm as a function of the barrier height and number of qubits adopted.

Variational network assessment
The first question we want to address is whether a RyRz variational form can produce a good trial wavefunction to approximate the classical computational result and if the optimization procedure is able to find its way towards such a configuration.
In order to give a reasonably general picture a few representative results have been reported in Tab. 1. All the data in Tab. 1 have been obtained from 60 independent VQE optimizations carried out on the Qiskit statevector simulator. For each evaluation the parameters of the RyRz variational form have been set randomly and optimized with the SPSA algorithm for a maximum of 600 iterations. The choice of the optimization algorithm has been dictated by a performance analysis carried out on different optimizers reported in appendix A. The expectation value has been evaluated at each step with an independent measurement for each of the Pauli strings.
We have chosen to report the minimum value obtained from the sampling as a measure of the heuristic network capability in generating a trial state reasonably close to the optimal configuration; while we have adopted the average converged result to quantify the performance of the optimizer in extracting the desired state from the overall accessible Hilbert space. We are aware that this is a non trivial analysis given that many factors contribute to the shape of the optimization landscape and the exploration of the underlying Hilbert space. Nevertheless, we will see as our considerations are in agreement with more structured studies in literature applied to parametrized quantum circuits [40,41].
As it can be seen from Tab. 1 increasing the number of qubits results in a steep rise of the average value error accompanied by a somewhat milder increase in the error on the minimum value that, even for 4 qubits, remains under the 5% limit. These data indicate how a RyRz variational form of depth d = 1 can be used to produce reasonable approximations of the target state with an error, on the expectation value, lower than 1.1% for 2 and 3 qubits systems. This finding is in agreement with the general analysis carried out by Sim. et al. [41] where parametrized quantum circuits similar to the ones adopted in this study have shown a good capability of uniformly represent the Hilbert space in which they are defined.
Conversely, the higher error on the distribution average is a clear sign of the challenging optimization process. It is known that such a complexity is not only arising from the dimensionality of the problem but it is strictly related to the phenomenon of quantum tipicality [42,43] meaning that the distribution of certain functions over the quantum states of a given Hilbert space is extremely peaked around a typical value. Specifically, for random quantum circuits such as the one considered in our case, the average value of the gradient objective function tends to zero and as the Hilbert dimension increases the more states will correspond to flat optimization landscape regions [44,45].
In order to clarify how the optimizer convergence impacts the overall VQE outcome let us observe how a state |φ obtained from a N -qubit RyRz variational form can be exactly reproduced, in the form of |φ ⊗ |0 , by a (N + 1)-qubit RyRz variational ansatz by properly setting the rotation parameters. This implies that: given two variational systems, composed respectively by N and Table 1: Results obtained by 60 independent VQE optimizations using the statevector simulations. Q indicates the number of qubits and ∆ r is the barrier height for the double mimimum dihedral. We report min(λ 1 ) as the minimum value obtained in the simulations, λ 1 as the average value, while λ ref 1 represents the theoretical target value for the considered basis set obtained with classical calculation. ε min (%) and ε avg (%) are respectively the percentage error on the minimum and on the average. N > N qubits, and considering the proper ordering of the basis-set, the accuracy of the best estimate accessible to the larger N -qubits system must be higher or at least equivalent to the accuracy of the estimate produced by the smaller N -qubits one. In order to test this idea a hierarchical procedure has been devised in which the optimized parameters of a N -qubit variational form have been used as the starting guess in a (N + 1)-qubit VQE procedure. A total of 500 VQE have been carried out using COBYLA as the optimizer. Differently from what observed in a regular VQE procedure, as the one reported in Tab. 1, a monotonic decrease in the minimum value is encountered in the hierarchical procedure; a minimum value of 1.516 is found for a 2 qubits system, a value of 1.484 is encountered for a 3 qubits system while a final value of 1.475 is encountered in the case of 4 qubits. A similar monotonic decrease is also observed in the case of the distribution average. We therefore verified that the accuracy of the obtained result increases by enlarging the computational space in the case the optimization procedure is guided by an educated guess. This observation is fundamental in understanding the nature of the results discussed above proving that the algorithmic bottleneck for the accuracy is the classical optimization procedure rather than the quantum portion. In conclusion, the trends emerging from these calculations are supporting the following considerations: the pure state distribution generated by the two entanglers allows one to reach the configuration corresponding to the global minimum of the optimization landscape but, on the other hand, we have seen that the landscape itself is affected by the dimension of the explorable Hilbert space. This drawback could be overcome developing physicallyinspired ansatze for the Smoluchowski operator similarly to what has been done in literature for the molecular Hamiltonian [38].
The results presented in this section will be useful in the following discussion: knowing the limit imposed by an ideal implementation will allow us to independently verify how various types of noise affect the overall stability of the method.

Statistics of finite measurements and quantum noise effects
At this point we continue our analysis by discussing the effects of different noise sources on the circuits presented above with the aim of answering the question: how does the introduced noise impact the expectation value measurement for a given circuit configuration? This is especially important in the case of nearly converged configurations in which the error introduced by noise can significantly affect the ability of the optimizer to operate. In order to investigate this point, a converged set of RyRz parameters has been used to study the effects of finite measurements and quantum noise on the expectation value. The optimized variational parameters have been obtained using the IBM Q Santiago quantum computer noise-model. To accumulate statistics, we have measured Γ with 1000 independent runs. Each Pauli string measurement was accomplished accumulating 20000 shots. Figure 4 reports the results of this investigation as two distinct distributions for the expectation value in the case of a sole sampling noise or with an addition of a noise model tailored on the IBM Q Santiago device. Figure 4: Expectation value distribution obtained from 1000 circuit (each one executed 20000 times) on the same converged RyRz parameter set for a 2 qubit system with ∆ r = 0.5k B T for the bi-stable dihedral. Blue distribution relates to the sampling under the simulated IBM Q Santiago hardware noise; orange distribution is due to finite measurements.
Looking at the distributions we can easily verify how the introduction of a noise model in the simulation causes both an increased distribution amplitude, with the loss of the peaked Gaussian profile observed for the noiseless simulator, and a shift of the average value. This last point is quite interesting since it highlights the profound difference existing between the two effects. The sampling effect stems from an inherent rounding error due to the limited number of measurement performed. As such, for a set of repeated measures it will appear as a Gaussian distribution centered around the exact expectation value for the state prepared by the variational form. Moreover, the variance of the sampling distribution tends to zero as the number of measurements increases. On the other hand, the simulated noise of a NISQ device alters the generation of the trial state producing slightly different wavefunctions. Since the trial state, obtained via the optimization procedure, is located near a minimum, most of the error in the state preparation will result in a higher expectation value estimate. The best value obtained in the presence of noise is, indeed, 0.7% higher than the best value obtained in an ideal calculation.
We now consider how the noise affects the quantum calculation of both different models, by changing the potential energy barrier of the conformational transition, and of different implementations, by changing the number of basis functions describing the target eigenfunction (i.e. using a different number of qubits).
Performance dependency on the potential barrier height For this evaluation we have performed 250 VQE routines starting with a random guess for the RyRz parameters. The variational parameters have been optimized using the SPSA algorithm under the effect of the Santiago quantum computer noise model. Every circuit has been transpiled before the simulation. In Fig. 5 we report the results where two qubits encode the composite dihedral states describing a system of three rotors in which ∆ r for the double minimum dihedral has been varied form 0.5 k B T to 3.0 k B T while ∆ nr for the single minimum dihedral was kept fixed at 1 k B T .
First of all we can notice how, as expected for activated processes, the kinetic constant decreases as we increase the potential energy barrier. Looking at the data we can observe the good agreement between the classically estimated target value (solid line) and the average results obtained with IBMQ Santiago noise model (downward triangles). In the same graph, we report (stars) the values of λ 1 obtained computing the expectation value on the IBMQ Santiago quantum computer starting from a set of variational parameters optimized using a VQE procedure carried out on the Qiskit Statevector Simulator. During the expectation value measurement on the IBMQ quantum hardware each circuit has been repeated, in order to build the required statistics, for the maximum allowed number of times (8192 runs). The obtained data, despite being computed from parameters obtained in a completely noiseless procedure, provide a natural metric to look at the results in term of the characteristic noise of a NISQ era quantum computer. These information, together with the ones coming from the noisy simulations, can be seen as a general proof of concept showing the feasibility to run the entire calculation on a NISQ device.
Bearing in mind that the quantum noise affects only the evaluation of the final expectation value we may notice that adding a further noise source causes sensible deviations from the noise-free average. To comment this results we leverage the intuition built from the statistical analysis in the previous section: the points collected on the IBM Q Santiago can be seen as samples of a third distribution generated by the noise of the actual quantum hardware. Under this light we may notice that in almost all cases we obtain absolute errors comparable with the maximum and mean deviations of the two distributions reported in Fig. 4.
This result on the one hand supports the effectiveness of the noise model in capturing the performance of real hardware, and on the other hand allows us to justify the greater deviation from the ideal average at the last point. Finally, we can also notice that the absolute width of the distribution remains substantially unchanged in the range of barriers explored. As a consequence, the relative error decreases moving from higher to smaller barriers passing from 13.6% at 3.0 k B T to 3.72% at 0.5 k B T . In the same figure (lower panel) the noise-free statevector results are reported as reference. Here we can appreciate with more the detail that, even though the distribution width is almost constant, the optimization is more sensitive with respect to the initial guess distribution for small barrier values.
Performance dependency on the basis set size Now that the accuracy of the algorithm has been discussed in terms of the physics of the simulated system, a final question still remains open: how does the performance of the protocol scales with the number of qubits employed? To answer to such a question we consider the same system of a three rotor chain with reactive barrier height of ∆ r = 0.5k B T . The probability density distribution has been described adopting an increasing number of basis elements from 4 up to 16 odd composite functions (see Eq. 17). This translates to an increasing number of qubits to encode the overall probability density distribution. The purpose of this analysis is not to assess the accuracy of the basis set size (which determines only a slight change in the kinetic rate constant) but to observe the effect of a noisy device upon performing the same calculation with an increasing number of quantum resources. The results obtained in such a simulation, together with the reference values, are reported in Fig. 6.
As we can see the distribution width rapidly increases with the number of qubits and so does also the average expectation value computed with both statevector and noisy simulators. This points out the increased complexity of the optimization procedure as previously mentioned in Sec. 4.1. For a 4 qubits system an error of the 37.6% is recovered for the noisy simulation while an error of the 23.3% is encountered in the case of a Statevector simulation. These data confirm the role of the noise in creating less accurate predictions but shows how a significant part of the error itself should be ascribed to a difficult convergence of the classical optimization as anticipated in section 4.1 . This massive performance degradation has also the effect of hiding the tiny margin for accuracy improvement induced by the basis-set extension. This is a critical aspect to consider in this early VQE implementations in which a good trade-off between algorithm stability and basis-set extension must be found. Further analyses of possible alternative variational networks could circumvent this issue.

Conclusions
In this paper we have presented a novel application for quantum computers. Particularly, we have addressed the solution of the Fokker-Planck-Smoluchowski eigenvalue problem exploiting its isomorphism with the Schrodinger equation. The proposed algorithm consists of three main ingredients: (i) a strategy to map the classical probability distribution onto the quantum register, (ii) the choice of a unitary ansatz to generate trial probability distributions (in the same spirit of the original VQEs for quantum systems) and (iii) the implementation of a classical optimization routine updating the unitary control parameters to move towards the exact configuration. Further investigation should follow along all the lines mentioned above.
Indeed, in order to devise a quantum algorithm useful for everyday applications and able to achieve quantum advantage with respect to a classical solution, a polynomial scaling of the quantum resources with the system size is desirable. Our analysis in Sec. 3 has shown how this has already partially been achieved considering the information storage point of view. On the other hand, a dense mapping such as the binary approach adopted insofar requires a number of measurements that scales exponentially with the number of rotors in the chain. This aspect can be mitigated by implementing different strategies that group sets of Pauli strings reducing the overall measurement count [37,[46][47][48]. For instance, considering the J diagonal elements, a single circuit execution is sufficient since all the strings are diagonal in the computational basis and can be measured simultaneously in a single post-rotation. Despite this, looking at systems of interest from an application point of view a step forward is still needed.
In a future work we will seek to overcome this issue by developing a mapping on the single-dihedral states which should reduce dramatically the number of terms in the FPS operator representation. In quantum chemistry parlance, this would correspond to a mapping on the spin-orbitals in place of a mapping on the Slater determinant basis.
Other matters to be addressed concern the unitary ansatz generating the trial probability distribution and the classical optimizer. These aspects have been considered in Sec. 4. There, we have shown that adopting a RyRz heuristic variational network allows to achieve good results for a small basis set. We believe that the results we have obtained in presence of a simulated noise can be further improved by the incorporation of a mitigation routine such as the zero noise error extrapolation [49]. However, we decided to not include also this feature in the implementation as it goes beyond the purpose of the present study of presenting a new application of hybrid algorithms in the field of stochastic processes.
At the same time, we have also noticed that the performance scales poorly with respect to the number of qubits and suffers of some problems deriving from the uniform and agnostic coverage of the underlying Hilbert space. In this direction we have suggested as a working strategy to lower the barrier of the optimization task by applying an hierarchical procedure in which the solution of the variational optimization with a smaller number of basis function is used as initial guess for the estimation of the kinetic rate constant using an higher number of qubits. Such an approach is in line with the development of algorithms in the NISQ era where small quantum processors are put to test by very shallow circuits; alternatively, we should look at further refinements of this implementation by developing physically-inspired ansatze as done in other contexts. Also this aspect will be investigated in a future work.
In essence, we achieved the proof of concept of the application of a quantum computer to solve stochastic dynamics. With some technical developments as discussed in the previous paragraphs, the proposed procedure brings stochastic dynamics amongst the promising applications of future quantum devices.
In what follows the optimizer selection and the effect that the entangler type has on the convergence behavior of the algorithm will be presented. Addressing the discussion of this problem is not trivial since the convergence profile of each VQE routine depends upon the selected set of initial parameters. In what follows we will proceed by analyzing, for different choices of optimizer and entangler, the convergence profile of a single VQE. This clearly represents only a single point in the multidimensional initial parameter space and consequently only comparative considerations can be done between different methods. Such a simple data representation, visible in Figure 7, can anyway already give a useful qualitative picture. For the present tests, a 3 qubit system is considered; both dihedral barriers (single and double) have been set to 1 k B T and a common set of randomly selected initial parameters has been kept fixed for all simulations. The built-in Qiskit SPSA optimizer together with the Nelder-Mead, SLSQP and COBYLA optimizer from the SciPy [50] Python module have been employed. The figure contains three plots that represent the profile of the VQE optimization in the case of various type of noise. The top panel contains data obtained using the noiseless statevector simulator. The center panel represents a noiseless simulator in which the number of measures on the same circuit has been set to 20000. The lower panel shows the results for the transpiled circuits simulated including the noise model of the IBMQ Santiago quantum computer. Also for this last case the number of runs has been set to 20000.
Analyzing Figure 7 many interesting comments can be made. First of all looking at the raw performance of the various optimizers one can easily observe how the Nelder-Mead algorithm consistently shows a slower convergence toward the target value if compared with the other optimization procedures. The SLSQP algorithm shows very good performance in the Statevector simulations converging in few iterations to the target value; however its sensitivity to noise prevent, as far as we can tell, its application to noisy simulations. Both COBYLA and SPSA shows consistently stable and reasonably good performances with the latter showing a more consistent behaviour and greater accuracy in simulations with higher noise. For this reason SPSA has been used as the default optimizer throughout the paper. Looking at the central panel one can easily see how little difference is observed between different entanglers in term of converged values with the full entangler resulting mildly more accurate. As it can be seen in the lower panel, the situation is rapidly reversed by the introduction of a realistic quantum computer noise model. The more complex structure of the full entangler joined with the less efficient transpilation, results in a slower converging and significantly less accurate result, clearly pointing to the linear entangler as the natural candidate for a realistic implementation.