Variational Quantum Simulation of Partial Differential Equations: Applications in Colloidal Transport

We assess the use of variational quantum imaginary time evolution for solving partial differential equations. Our results demonstrate that real-amplitude ansaetze with full circular entangling layers lead to higher-fidelity solutions compared to those with partial or linear entangling layers. To efficiently encode impulse functions, we propose a graphical mapping technique for quantum states that often requires only a single bit-flip of a parametric gate. As a proof of concept, we simulate colloidal deposition on a planar wall by solving the Smoluchowski equation including the Derjaguin-Landau-Verwey-Overbeek (DLVO) potential energy. We find that over-parameterization is necessary to satisfy certain boundary conditions and that higher-order time-stepping can effectively reduce norm errors. Together, our work highlights the potential of variational quantum simulation for solving partial differential equations using near-term quantum devices.

Although linear differential equations can be solved by the quantum linear solver algorithm (Berry et al., 2017;Harrow et al., 2009), the required resources are out of reach of the current noisy intermediate-scale quantum devices (Lau et al., 2022;Bharti et al., 2022;Preskill, 2018). In fact, practical near-term quantum algorithms are limited to those designed for short circuit depths, such as variational quantum algorithms (VQAs) (Cerezo et al., 2021), which use parameterized ansätze to optimize cost functions via variational updating.
VQAs can largely be classified into two categories, namely, optimization and simulation , each offering unique approaches to solving PDEs. Variational quantum optimization (VQO) aims to optimize a static target cost function through parameter tuning, an example of which is the popular variational quantum eigensolver (VQE) (Peruzzo et al., 2014) for minimizing energy states in the field of quantum chemistry. This led to the development of the variational quantum linear equation solver (Bravo-Prieto et al., 2019;Huang et al., 2021;Xu et al., 2021) for systems of linear equations and the variational quantum Poisson solver (Liu et al., 2021;Sato et al., 2021). Evolution of the Poisson equation allows parabolic PDEs to be solved through implicit time-stepping (Leong et al., 2022), which requires quantum information to be updated and encoded at each time-step.
On the other hand, variational quantum simulation (VQS) aims to simulate a dynamical quantum process, such as the Schrödinger time evolution (Li and Benjamin, 2017). This allows certain PDEs to be solved efficiently using imaginary quantum time evolution Endo et al., 2020;Yuan et al., 2019), including the Black-Scholes equation for option pricing (Miyamoto and Kubo, 2021;Radha, 2021;Stamatopoulos et al., 2020) and stochastic differential equations for stochastic processes . Recent work on the Feynman-Kac formulation (Alghassi et al., 2022) generalizes quantum simulation of parabolic PDEs, paving the way for potential near-term applications.
In this study, we explore applications of VQS (Miyamoto and Kubo, 2021;Alghassi et al., 2022) in solving PDEs, including the Smoluchowski equation for colloidal physics, with an emphasis on potential and non-homogeneous terms oft-neglected in quantum simulations. We select for high-fidelity real-amplitude ansätze, assess time complexity and propose an efficient encoding scheme for idealized pulse functions, as a proof of concept towards practical implementation of quantum simulation. HFF 2. Variational quantum simulation 2.1 Evolution equation Consider a 1-dimensional (1D) evolution equation expressed in the Feynman-Kac formulation (Alghassi et al., 2022): is a function of space x and time t, and a, b and c are the coefficients to the second-, firstand zeroth-order derivative terms in x, respectively (bold symbols denote vectors in space is a non-homogeneous source term and u 0 is the initial condition. Following , we rewrite equation (1) in Dirac notation [1]: where H(t): ¼ a@ xx þ b@ x þ c is the Hamiltonian operator, possibly non-Hermitian, and F (t) is a linear operator satisfying F (t)j0i ¼ jf (t) where u 0 (t) is a normalization parameter. To minimize the distance kju t ð Þi À jũ u t ð Þ ð Þik, we apply the McLachlan's variational principle : where jjvjj :¼ ffiffiffiffiffiffiffiffiffi ffi hvjvi p denotes the Euclidean norm and d denotes infinitesimal variation. This yields a system of ordinary differential equations: where _ u t ð Þ :¼ @ t h t ð Þ. The left-hand side matrix: and the right-hand side vector: Variational quantum simulation can be evaluated parametrically on quantum circuits (McArdle et al., 2019). See Appendix 1 for details.
With A and C specified, parameters h are evolved in time using the forward Euler method as: up to N t time-steps in each Dt. Higher-order methods, such as Runge-Kutta, are also available. Because the matrix A may be ill-conditioned, successful inversion may depend on methods such as the Moore-Penrose inverse or Tikhonov regularization (McArdle et al., 2019). We find that least-squares minimization with a 10 À6 cutoff is sufficient for stable solutions (Fontanela et al., 2021).

Decomposition of Hamiltonian
The Hamiltonian operator H introduced in equation (2) can be simplified through elimination of the skew-Hermitian term b@ x using substitution methods (Fontanela et al., 2021), such as u(t) ¼ e g v(t), where g is a function of a and b. If g(a, b) were constant in time, then the Hamiltonian operator reduces to Alghassi et al. (2022): where I is the identity operator. The potential vector is: where the last term can be neglected if {a; b} is independent of x. The Hamiltonian operator can be discretized in the space interval Dx, and decomposed into a linear combination of terms as: where 1 ¼ (1, 1, . . ., 1, 1) is the all-ones vector and I 0 ¼ j0ih0j. For the Neumann boundary condition, all six terms {H 1 , . . ., H 6 } are required. For the periodic boundary condition, only the first four terms {H 1 , . . ., H 4 } are required, and for the Dirichlet boundary condition, the first five terms {H 1 , . . ., H 5 } are required. Note that as an observable in the first term, the potential vector u does not increase the quantum complexity; measurement of the existing H 1 suffices to evaluate the expectation value of the potential. The operator S denotes the n-qubit cyclic shift operator (Sato et al., 2021): which can be implemented as a product of k-qubit Toffoli gates, for k in the range 1, . . ., n (see e.g. Sato et al., 2021, Figure 2).

Ansatz selection
For optimal algorithmic performance, a good choice of ansatz is crucial (Tilly et al., 2022;You et al., 2021). For PDEs that admit only real solutions, it is preferable to use a realamplitude ansatz formed by n l repeating blocks, each one consisting of a parameterized layer with one R Y rotation gate on each qubit, followed by an entangling layer with CNOT gates between consecutive qubits (Alghassi et al., 2022). Here, we consider two options for customization: the first between linear and circular entanglement and the second with or without an unentangled parameterized layer as the final block n l , as shown in Figure 1. We note that the circular entanglement with a final unentangled layer [ Figure 1(c)] is a popular choice of an ansatz for VQS Alghassi et al., 2022). For benchmarking, we perform numerical experiments on the various ansätze ( Figure 1) to solve a simple 1D heat or diffusion equation, expressed in Dirac notation as:

Variational quantum simulation
The initial trial state is set as a reverse step function (Sato et al., 2021): which can be implemented in practice by setting the final parameterized layer l as Figure 2(a) and 2(b) show the time evolution of the step function for four-qubit real-amplitude ansätze with four layers using time-step Dt ¼ 10 À4 up to T ¼ 10 À2 . We measure the fidelity of the VQS solution obtained from each ansatz compared to the classical solution and define the trace error as: Similarly, we define the norm error as: HFF where u 0 (t) is the normalization parameter as previously defined. Figure 2(c) and (d) shows the mean trace and norm errors depending on the number of ansatz layers n l using time-step Dt ¼ 10 À4 up to T ¼ 0.01, for periodic and Dirichlet boundary condition, the latter shown as closed symbols for jf(t)i ¼ 0.
The circular, fully entangled ansatz [ Figure 1(d)], here termed full circular ansatz for short, was found to outperform other ansätze, requiring fewer parameters for the same solution fidelity. For four qubits, the full circular ansatz is the only one able to produce a solution overlap with only two or three layers, which is less than the minimum required for convergence n l < 2 n /n. For five qubits, it delivered reduced solution and norm errors compared to other ansätze, independent of the number of layers. In this benchmark, the additional term introduced by the Dirichlet boundary condition does not diminish the superior performance of the full circular ansatz.

Initialization
An initial quantum state ju(0)i can be prepared through classical optimization and accepting converged solutions whose norms fall below a specified threshold (Fontanela et al., 2021;Alghassi et al., 2022) or direct encoding using quantum generative adversarial networks (Zoufal et al., 2019). In most cases, quantum encoding is cost-prohibitive, and subexponential encoding can be achieved only under limiting conditions Mitsuda et al., 2022).
The Dirac delta function is a popular initial probability distribution found in Fokker-Planck equations Alghassi et al., 2022). To encode the state jxi in the computational basis jxi ¼ Þi for an input state j0i. It turns out that for a full circular ansatz [ Figure 1(d)], encoding jxi does not necessarily require costly optimization. To access a given state jxi, one can search for a parameterized layer n lk such that a p bit-flip rotation on an R y u n l Àk 1;n ½ gate (or gates) yields an input state jx 0 i, which transforms to the output state after k circular entangling layers, i.e. C k n jx 0 i ¼ jxi, where the matrix C n represents a single circular entangling layer (Appendix 2). Figure 3 shows that for a four-qubit full circular ansatz, all 2 4 -1 ¼ 15 jxi states can be encoded by a single p bit-flip rotation of an R y gate within four parameterized layers.

Time complexity
To assess the time complexity of the VQS algorithm, we estimate the number of quantum circuits required per time-step as: where n p is the number of ansatz parameters (n p ¼ n ln for a real-amplitude ansatz) and n h is the number of terms in the Hamiltonian [equation (11)]. Likewise, the number of circuits required per time-step for a VQE implementation (Leong et al., 2022) can be estimated as: where n it is the number of iterations taken by the classical optimizer. Hence, the VQS algorithm is comparable with VQE in terms of circuit counts, if the number of ansatz Variational quantum simulation parameters is roughly double the expected number of iterations required for VQE, i.e. n p % 2n it ) n h .
For each circuit, the time complexity scales as (Sato et al., 2021): where d ansatz $ O(n l ) is the depth of the ansatz, d shift $ O(n 2 ) is the depth of the shift operator, and the denominator O(« 2 ) reflects the number of shots required for estimated expectation values up to a mean squared error of « 2 . Another consideration is the depth for amplitude encoding d enc , which can range from O(n 2 ) to O(2 n ). For VQS, encoding is performed once during initialization, but unlike VQE, repeated encoding is not necessary for time-stepping (Leong et al., 2022).
To solve an evolution PDE [e.g. equation (1)], a classical algorithm iterates a matrix of size 2 n Â 2 n , compared to a n p Â n p matrix for VQS, suggesting comparable performance at n l % 2 n /n.

Colloidal transport
With the VQS framework in place, one can explore applications in solving PDEs, such as heat, Black-Scholes and Fokker-Planck equations listed in Alghassi et al. (2022). In this  (9)], an aspect oft-neglected in quantum simulations.

Smoluchowski equation
Consider a spherical colloidal particle of radius a near a planar wall (Torres-Díaz et al., 2019). The generalized Smoluchowski equation (Smoluchowski, 1916) describes the probability p (h, t) of locating the particle at h, the distance of the particle centre from the wall at time t, as: where D is the diffusivity matrix. U(h) is the Derjaguin-Landau-Verwey-Overbeek (DLVO) sphere-wall interaction energy (Bhattacharjee et al., 1998), which is the sum of the electric double-layer and van der Waal's interaction energies, expressed as: where U ¼ U(h) and H ¼ (ha)/a is the dimensionless separation distance between the particle and wall. The electric double-layer coefficient, normalized inverse Debye length and van der Waal's coefficient are, respectively: where « is the relative permittivity of the medium, « 0 is the permittivity of free space, z v is the ionic valence, e is the electron charge, k B is the Boltzmann's constant and T is the temperature. z p and z w are the zeta potentials on the colloidal particle and wall, respectively. I is the ionic strength and A H is the Hamaker constant. With that, the first and second derivatives of the interaction energy in separation distance are: Rescaling time t ¼ tD/a 2 , we rewrite equation (20) in dimensionless form, which gives the evolution of the probability p(t) ¼ p(H, t) as: Variational quantum simulation which follows from equation (1), where a ¼ 1, b ¼ U 0 , c ¼ U 00 and f ¼ 0. Suitable boundary conditions are the absorbing condition on the surface at p(0) ¼ 0 (Dirichlet) and the no-flux condition in the far field @p(1)/@H ! 0 (Neumann) (Torres-Díaz et al., 2019).
Substituting p(t) ¼ q(t)e -U/2 (Fontanela et al., 2021), we express in Dirac notation: with the Hamiltonian operator: and the potential term: which can be evaluated classically [equation (23)] and implemented as a quantum observable u T°I . With Dirichlet-Neumann boundary conditions enforced, the Hamiltonian operator is decomposed as: where: e 2 n ¼ 0; 0; . . . ; 0; 1 ð Þ : 3.1.1 Potential-free case (w ¼ 0). Consider first the potential-free case where colloid-wall interactions are absent (w ¼ 0 where the initial impulse state is centred at jr 0 i ¼ j2 n-1 i. Using a full circular ansatz with six to eight layers, we evolve the initial pulse on a 2 n ¼ 16 grid using time-step Dt ¼ 10 À4 for early times up to 10 À2 [ Figure 4(a)] and late times up to T ¼ 10 À1 [Figure 4(b)]. The former is characterized by the spreading of the probability density due to diffusion, and the latter by the constraints imposed by the asymmetric boundary conditions, which reduces solution fidelity.
To assess the costs of over-parameterization, we calculate the mean trace error [ Figure 4(c)] and norm error [Figure 4(d)] depending on the total number of circuits required N q for the VQS with a run-time of T ¼ 10 À1 . Figure 4(c) shows that the mean trace error is insensitive to number of ansatz layers up to six and time-steps up to 5 Â 10 À4 ; it is reduced only with further HFF increase in the number of ansatz layers n l > 6, leading to optimal scaling of « tr $N À2 q . Closed symbols show results using a higher-order Runge-Kutta time-stepping in place of first-order Euler time-stepping [equation (8)]. We see that the cost scaling of the mean trace error is relatively unaffected by higher-order time-stepping, due to the additional circuit count required for four Runge-Kutta iterations per time-step. The peak trace errors follow a weaker scaling of max « tr ð Þ $N À0:8 q . Using Euler time-stepping, the mean norm error scales as « norm $N À0:8 q , regardless of n l [ Figure Figure 4(d)]. This cost scaling improves significantly up to « norm $N À3:4 q using Runge-Kutta time-stepping for circuits with n l > 6. Note, however, that this improvement does not extend to the peak norm errors, whose cost scaling remain as max « norm ð Þ $N À0:8 q , regardless of time-stepping scheme. probability densitŷ r H ð Þ profiles without interaction (U ¼ 0) on a 2 n ¼ 16 grid using time-step Dt ¼ 10 À4 plotted in increments of (a) 2 Â 10 À3 and (b) 2 Â 10 À2 up to  A, Z, k). In the presence of colloid-wall interactions, the DLVO potential term w depends minimally on three parameters, specifically A, Z and k [equation (22)]. Following the potential-free case (n ¼ 4, H [ [0, 1], Dt ¼ 10 À4 ), we perform VQS including w using eight ansatz layers in time t [ [0, T].
In the absence of the electric double layer (Z ¼ 0), the DLVO potential w(A) depends on only the van der Waal's interaction energy, assumed here to be attractive. Figure 5(a) shows that the DLVO potential w(H) profiles scaled by the square of the interval DH 2 for A [ {0.05, 0.1, 0.2, 0.5} is only short-ranged in H, so the quantum solution jri is insensitive to A. Recall, however, the earlier substitution p(t) ¼ q(t)e -U/2 , such that the actual solution p depends on the longer-ranged interaction energy U(H) [equation (21)], as shown in Figure 5(a, inset). Indeed, space-time plots show that the colloidal probability density p(H, t) for A ¼ 0.05 up to T ¼ 0.1 is depleted near wall [ Figure 5(c)] compared to the potential-free (w ¼ 0) case [ Figure 5(b)]. Increasing A further increases the depletion range [ Figure 5(d)].
Otherwise, the DLVO potential w(A, K, k) includes the electric double layer interaction energy, assumed here to be repulsive. For A ¼ 0.5, Figure 6(a) shows that the DLVO potential shows short-ranged dependence on Z and k. However, p depends on the longer-ranged interaction energy U(H) that can be either attractive or repulsive as shown in Figure 6(a, inset). A space-time plot of the colloidal probability density p(H, t) for {Z, k} ¼ {10, 10} shows long- 3.1.3 Trace and norm errors. Here, we characterize the effect of DLVO potential on the solution fidelity in time using the trace error « trace (t) [equation (15)] and the norm error « norm (t) [equation (16)]. Figure 7 shows that « trace (t) peaks and decreases during the early diffusion phase [ Figure 4(a)], then peaks and decreases again as the normalized probability densityr approaches a steady-state profile constrained by the imposed asymmetric boundary conditions [ Figure 4(b)]. Parametric analyses suggest that the electric double layer coefficient Z has the strongest effect on « trace (t) [Figure 7(b)]. In contrast, « norm (t) tends towards a steady state regardless of the evolution of probability density. Parametric analyses suggest that « norm is affected by the local depletion ofr but insensitive to the magnitudes of A [ Figure 7 Thus concludes our analysis of the potential term in equation (28) in Smoluchowski equation. What usually follows are calculations of survival probability, the probability that the colloidal particle has not reached the wall, the mean first passage time distribution and the mean rate of change of survival probability. Because they do not involve any quantum computation, they are outside the scope of this study. Interested readers are referred to Torres-Díaz et al. (2019).

Einstein-Smoluchowski equation
The general PDE introduced in equation (1) includes a non-homogeneous source term f, which is not admissible in Smoluchowski's description of colloidal probability density. To explore the effects of a source term, we switch over to the analogous Einstein-Smoluchowski equation (Cejas et al., 2019;Leong et al., 2023): where the operator F ¼ X n imposes a unit source in the far field, increasing the required number of quantum circuits by n p þ 1 [equation (17)] per time-step. The number of additional circuits scales with the number of unitaries required to express F . 3.2.1 Initialization. We seek a parameterized ansatz that encodes a Heaviside step function centred at j2 n-1 i: For a full circular ansatz [ Figure 1(d)], this can be encoded on a minimum of two R Y parameterized layers by setting the final layer as u n l 1;n ½ ¼ p 2 and the second R Y of the preceding layer as u n l À1 , where a reversal in the sign produces a step-down function instead.
3.2.2 Solutions and errors. We perform VQS on a 2 n ¼ 16 grid using time-step Dt ¼ 10 À4 as before, but on a full circular ansatz with five layers, which is already shown to yield high-  Figure 8(a) shows how the normalized concentration c evolves from the initial step function for the potential-free case (U ¼ 0). In the absence of an electric double layer (Z ¼ 0), strong attractive van der Waal's energy leads to fast convergence towards a steady-state profile [ Figure 8(b)]. Increasing Z shifts the steady-state concentration profile near wall [ Figure 8(c)], whereas decreasing k increases the depletion range [ Figure 8(d)]. Both trace and norm errors (Figure 8 insets) decay in time t towards convergence with « trace peaking earlier than « norm .

Conclusion
Currently, neither VQO nor simulation is capable of realizing an advantage for solving PDEs over classical methods (Anschuetz and Kiani, 2022), but that gap is closing fast (Tosti Balducci et al., 2022). For VQS, a significant progress has been made since the advent of imaginary time evolution (McArdle et al., 2019) notably in the field of quantum finance (Fontanela et al., 2021;Miyamoto and Kubo, 2021;Kubo et al., 2021).
Here, we list a formal approach to solving a 1D evolution PDE [equation (1)]: {@ t u(t), @ xx u(t)} terms handled using variational quantum imaginary time evolution. @ x u(t) term eliminated through substitution methods, such as u(t) ¼ e g v(t). Variational quantum simulation u(t) term included in the Hamiltonian H without additional complexity cost. f(t) term realized by an additional set of complementary circuits, whose complexity depends on F .
Superior performance of VQS is contingent on two factors: selection of ansatz and initialization of parameters. Comparing real-amplitude ansätze (Section 2.3), we found that the full circular ansatz significantly outperformed not only linear entangled ansätze but also the popular circularly entangled ansatz but with the final parameterized layer unentangled Alghassi et al., 2022). The advantage in solution fidelity persists over multiple parametric layers, which suggests that unentangled parameterized gates reduce overlap with quantum states that are characteristic of PDE solutions. For an initial state resembling a Dirac delta function (Section 2.4), we found that full circular ansatz can be mapped parametrically to a desired state jxi, thus reducing subsequent impulse encodings to only a trivial lookup.
As a proof-of-concept, we performed VQS to simulate the transport of colloidal particle to an absorbing wall as described by the Smoluchowski equation (Section 3.1) and found high solution fidelity during the initial spreading of the probability distribution. However, to satisfy the asymmetric boundary conditions, additional parameter layers are required, for example, up to -six to eight layers for a four-qubit problem. Higher-order time-stepping such as Runge-Kutta method can reduce norm errors more effectively than overparameterization for the same time complexity.
With near-wall DLVO potentials, we found that the van der Waal's interaction impacts VQS mainly through the potential w(A) of the Hamiltonian, whereas the electric double layer interaction affects the solution mainly through the factor e -U/2 obtained from change of variables. Simulations of colloidal concentration with unit boundary source in the far field (Section 3.2) require additional circuit evaluations equal to approximately half the number of parameters. Interestingly, this cost is offset by the fact that fewer parameters are required, here, for example, five layers for a four-qubit problem.
Overall, we find VQS an efficient tool for applications in colloidal transport because DLVO potentials do not incur additional costs in terms of quantum complexity. Compared to VQE (Leong et al., 2022), VQS enjoys significant advantages in that it does not require repeated encodings and iterative optimization loops. In terms of scalability, we found that the accuracy of quantum simulation not only depends on the number of qubits but also on the imposed boundary and the initial conditions. As with other gradient-based neural networks, VQS potentially suffers from barren plateau problems, which are exemplified by vanishing gradients on flat energy landscapes (McClean et al., 2018) and exacerbated by quantum circuits with high expressivity (Holmes et al., 2022). Mitigation strategies for barren plateaus remain an active area of research (Patti et al., 2021).
Future work can include extension to 2D model for non-spherical colloids (Torres-Díaz et al., 2019), optimal ansatz architecture (Tang et al., 2021) and initial state preparation Zoufal et al., 2020). Notes 1. For an introduction to quantum computation and Dirac notation, we refer the reader to Nielsen and Chuang (2010).
2. Not to be confused with the cyclic shift operator in equation (12), also denoted by S.