Variational quantum simulation of partial differential equations: applications in colloidal transport

Abstract

Purpose

This study aims to assess the use of variational quantum imaginary time evolution for solving partial differential equations using real-amplitude ansätze with full circular entangling layers. A graphical mapping technique for encoding impulse functions is also proposed.

Design/methodology/approach

The Smoluchowski equation, including the Derjaguin–Landau–Verwey–Overbeek potential energy, is solved to simulate colloidal deposition on a planar wall. The performance of different types of entangling layers and over-parameterization is evaluated.

Findings

Colloidal transport can be modelled adequately with variational quantum simulations. Full circular entangling layers with real-amplitude ansätze lead to higher-fidelity solutions. In most cases, the proposed graphical mapping technique requires only a single bit-flip with a parametric gate. Over-parameterization is necessary to satisfy certain physical boundary conditions, and higher-order time-stepping reduces norm errors.

Practical implications

Variational quantum simulation can solve partial differential equations using near-term quantum devices. The proposed graphical mapping technique could potentially aid quantum simulations for certain applications.

Originality/value

This study shows a concrete application of variational quantum simulation methods in solving practically relevant partial differential equations. It also provides insight into the performance of different types of entangling layers and over-parameterization. The proposed graphical mapping technique could be valuable for quantum simulation implementations. The findings contribute to the growing body of research on using variational quantum simulations for solving partial differential equations.

Citation

Leong, F.Y., Koh, D.E., Ewe, W.-B. and Kong, J.F. (2023), "Variational quantum simulation of partial differential equations: applications in colloidal transport", International Journal of Numerical Methods for Heat & Fluid Flow, Vol. 33 No. 11, pp. 3669-3690. https://doi.org/10.1108/HFF-05-2023-0265

1. Introduction

Partial differential equations (PDEs) are fundamental to solving important problems in engineering and science. With the advent of nascent quantum computers, finding new efficient quantum algorithms and hardware for solving PDEs has become an active area of research (Tosti Balducci et al., 2022; Jin et al., 2022; Leong et al., 2022; Pool et al., 2022) in disciplines ranging from fluid dynamics (Budinski, 2021; Gaitan, 2020; Steijl, 2022; Steijl and Barakos, 2018; Griffin et al., 2019; Li et al., 2023), heat conduction (Liu et al., 2022) and electromagnetics (Ewe et al., 2022) to quantitative finance (Fontanela et al., 2021) and cosmology (Mocz and Szasz, 2021).

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 (Endo et al., 2021), 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., 2021; 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 (Kubo et al., 2021). 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.

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):

(1) $$\frac{\partial \mathit{u}\left(t\right)}{\partial t}=a\frac{{\partial }^2\mathit{u}\left(t\right)}{\partial x^2}+b\frac{\partial \mathit{u}\left(t\right)}{\partial x}+c\mathit{u}\left(t\right)+\mathit{f}\left(t\right),\mathit{u}\left(0\right)={\mathit{u}}_0,$$

where u(t) = u(x, t) is a function of space x and time t, and a, b and c are the coefficients to the second-, first- and zeroth-order derivative terms in x, respectively (bold symbols denote vectors in space x). f (t) = f (x, t) is a non-homogeneous source term and u₀ is the initial condition. Following (Kubo et al., 2021), we rewrite equation (1) in Dirac notation [1]:

(2) $$\frac{\partial \left|u\right(t)〉}{\partial t}=\mathcal{H}\left(t\right)\left|u\right(t)〉+\mathcal{F}(t\left)\right|0〉,\left|u\right(t=0)〉=|u_0〉,$$

where $\mathcal{H}$(t): = a∂_xx + b∂_x + c is the Hamiltonian operator, possibly non-Hermitian, and $\mathcal{F}$(t) is a linear operator satisfying $\mathcal{F}$(t)|0〉 = |f(t)〉. The non-homogeneous operator $\mathcal{F}$(t) can be expressed as a sum of unitaries. Using VQS (McArdle et al., 2019), the state |u(t)〉 can be approximated by an unnormalized trial state $\left|\tilde{u}\right(\theta \left(t\right))〉$ formed by a set of parameterized unitaries {R_k}_k∈[N] with N parameters:

(3) $$\left|\tilde{u}\right(\mathit{\theta }\left(t\right))〉:={\theta }_0(t\left)R_1\right({\theta }_1\left(t\right)\left)R_2\right({\theta }_2\left(t\right))\cdots R_N({\theta }_N\left(t\right)\left)\right|0〉,$$

where θ₀(t) is a normalization parameter. To minimize the distance $\Vert \left|u\right(t)〉-|\tilde{u}\left(\theta \right(t\left)\right)〉\Vert$, we apply the McLachlan’s variational principle (Yuan et al., 2019):

(4) $$\delta \left|\right|\frac{\partial }{\partial t}\left|\tilde{u}\right(\mathit{\theta }\left(t\right))〉-\mathcal{H}(t\left)\right|\tilde{u}\left(\mathit{\theta }\right(t\left)\right)〉-\left|\mathit{f}\right(t)〉\left|\right|=0,$$

where $\left|\right|v\left|\right|:=\sqrt{〈v|v〉}$ denotes the Euclidean norm and δ denotes infinitesimal variation. This yields a system of ordinary differential equations:

(5) $${\displaystyle \sum _{j=0}^N{A}_{ij}}{\dot{\theta }}_j\left(t\right)=C_i,i=1,\dots ,N,$$

where $\dot{\theta }\left(t\right):={\partial }_t\mathit{\theta }\left(t\right)$. The left-hand side matrix:

(6) $$A_{ij}=\Re \{\begin{array}{ll}\frac{\partial 〈\tilde{u}\left(\mathit{\theta }\right(t\left)\right)|}{\partial {\theta }_i}\frac{\partial \left|\tilde{u}\right(\mathit{\theta }\left(t\right))〉}{\partial {\theta }_j},\hfill & \text{\hspace{0.17em}if\hspace{0.17em}}0<i\le j\le N,\hfill \\ 〈\tilde{u}\left(\mathit{\theta }\right(t\left)\right)|\frac{\partial \left|\tilde{u}\right(\mathit{\theta }\left(t\right))〉}{\partial {\theta }_j},\hfill & \text{\hspace{0.17em}if\hspace{0.17em}}0=i\le j\le N,\hfill \\ 1,\hfill & \text{if\hspace{0.17em}}i=j=0\hfill \end{array}$$

and the right-hand side vector:

(7) $$C_i=\Re \{\begin{array}{ll}\frac{\partial 〈\tilde{u}\left(\mathit{\theta }\right(t\left)\right)|}{\partial {\theta }_j}\mathcal{H}\left(t\right)\left|\tilde{u}\right(\mathit{\theta }\left(t\right))〉+\frac{\partial 〈\tilde{u}\left(\mathit{\theta }\right(t\left)\right)|}{\partial {\theta }_n}\mathcal{F}(t\left)\right|0〉,\hfill & 0<i\le N,\hfill \\ 〈\tilde{u}\left(\mathit{\theta }\right(t\left)\right)\left|\mathcal{H}\right(t\left)\right|\tilde{u}\left(\mathit{\theta }\right(t\left)\right)〉+〈\tilde{u}\left(\mathit{\theta }\right(t\left)\right)\left|\mathcal{F}\right(t\left)\right|0〉,\hfill & i=0\hfill \end{array}$$

can be evaluated parametrically on quantum circuits (McArdle et al., 2019). See Appendix 1 for details.

With A and C specified, parameters θ are evolved in time using the forward Euler method as:

(8) $$\mathit{\theta }(t+\Delta t)\leftarrow \mathit{\theta }\left(t\right)+\Delta t\left[A{\left(t\right)}^{-1}·C\left(t\right)\right],$$

up to N_t time-steps in each Δt. 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⁻⁶ cutoff is sufficient for stable solutions (Fontanela et al., 2021).

2.2 Decomposition of Hamiltonian

The Hamiltonian operator $\mathcal{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^gv(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):

(9) $$\mathcal{H}=\mathcal{a}\frac{{\partial }^2}{\partial x^2}+{\mathit{\phi }}^T\circ \mathbb{I},$$

where $\mathbb{I}$ is the identity operator. The potential vector is:

(10) $$\mathit{\phi }=\mathcal{c}-\frac{{\mathcal{b}}^2}{4\mathcal{a}}-\frac{\mathcal{a}}{2}\frac{\partial }{\partial x}\frac{\mathcal{b}}{2\mathcal{a}},$$

where the last term can be neglected if { $a,b$} is independent of x. The Hamiltonian operator can be discretized in the space interval Δx, and decomposed into a linear combination of terms as:

(11) $$\begin{array}{l}\mathcal{H}={\left[\mathit{\phi }-\frac{\mathcal{a}}{\Delta x^2}1\right]}^T\circ \underset{H_1}{\underbrace{I^{\otimes n}}}+\frac{\mathcal{a}}{\Delta x^2}\{\underset{H_2}{\underbrace{I^{\otimes n-1}\otimes X}}\hspace{0.17em}+\hspace{0.17em}{S}^{†}[-\underset{H_3}{\underbrace{I^{\otimes n}}}+\underset{H_4}{\underbrace{I^{\otimes n-1}\otimes X}}-\underset{H_5}{\underbrace{I_0^{\otimes n-1}\otimes X}}+\underset{H_6}{\underbrace{I_0^{\otimes n-1}\otimes I}}]S\},\hfill \end{array}$$

where 1 = (1, 1, …, 1, 1) is the all-ones vector and $I_0=|0〉\text{<mo>〈</mo> <mn>0</mn> <mo>|</mo>}$. For the Neumann boundary condition, all six terms {H₁, …, H₆} are required. For the periodic boundary condition, only the first four terms {H₁, …, H₄} are required, and for the Dirichlet boundary condition, the first five terms {H₁, …, H₅} are required. Note that as an observable in the first term, the potential vector φ does not increase the quantum complexity; measurement of the existing H₁ suffices to evaluate the expectation value of the potential.

The operator S denotes the n-qubit cyclic shift operator (Sato et al., 2021):

(12) $$S={\displaystyle \sum _{i=0}^{2^n-1}|}(i+1)\text{mod}{2}^n〉〈i|,$$

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).

2.3 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 real-amplitude 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 (Kubo et al., 2021; 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:

(13) $$\frac{\partial }{\partial t}|u\left(t\right)〉=\frac{{\partial }^2}{\partial x^2}|u\left(t\right)〉+|f\left(t\right)〉,|u\left(0\right)〉=|u_0〉,$$

in space x ∈ [0, 1] and time t ∈ [0, T].

The initial trial state is set as a reverse step function (Sato et al., 2021):

(14) $$|u_0〉=H^{\otimes n}\left(X\otimes I^{\otimes n-1}\right)|0〉,$$

which can be implemented in practice by setting the final parameterized layer l as $R_Y{\left(-\frac{\pi }{2}\right)}^{\otimes n}$ with entanglement, or $R_Y\left(-\frac{\pi }{2}\right)\otimes R_Y{\left(\frac{\pi }{2}\right)}^{\otimes n-1}$ without. 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 Δt = 10⁻⁴ up to T = 10⁻².

We measure the fidelity of the VQS solution obtained from each ansatz compared to the classical solution and define the trace error as:

(15) $${\epsilon }_{\text{tr}}\left(t\right)=\sqrt{1-{\left|〈{\widehat{u}}_c\left(t\right)|u\left(t\right)〉\right|}^2},$$

Similarly, we define the norm error as:

(16) $${\epsilon }_{\text{norm}}\left(t\right)=\left|1-\frac{{\theta }_0\left(t\right)}{\Vert {\mathit{u}}_c\left(t\right)\Vert }\right|,$$

where θ₀(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 Δt = 10⁻⁴ up to T = 0.01, for periodic and Dirichlet boundary condition, the latter shown as closed symbols for |f(t)〉 = 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. 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.

2.4 Initialization

An initial quantum state |u(0)〉 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 sub-exponential encoding can be achieved only under limiting conditions (Nakaji et al., 2022; Mitsuda et al., 2022).

The Dirac delta function is a popular initial probability distribution found in Fokker–Planck equations (Kubo et al., 2021; Alghassi et al., 2022). To encode the state |x〉 in the computational basis $\left\{|\mathit{x}〉={\otimes }_i^n|x_i〉\right\}$ with x_i ∈ {0, 1}, one seeks a parameterized ansatz $\left|\tilde{u}\right(\mathit{\theta }\left(0\right))〉$ for an input state |0〉.

It turns out that for a full circular ansatz [Figure 1(d)], encoding |x〉 does not necessarily require costly optimization. To access a given state |x〉, one can search for a parameterized layer n_l – k such that a π bit-flip rotation on an $R_y\left({\theta }_{[1,n]}^{n_l-k}\right)$ gate (or gates) yields an input state |x₀〉, which transforms to the output state after k circular entangling layers, i.e. $C_n^k|x_0〉=|x〉$, 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⁴ – 1 = 15 |x〉 states can be encoded by a single π bit-flip rotation of an R_y gate within four parameterized layers.

2.5 Time complexity

To assess the time complexity of the VQS algorithm, we estimate the number of quantum circuits required per time-step as:

(17) $$n_q^{\text{VQS}}\approx \left(n_p+1\right)\left(\frac{n_p}{2}+n_h\right),$$

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:

(18) $$n_q^{\text{VQE}}\approx n_{\text{it}}\left(n_p+n_h\right),$$

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 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):

(19) $$t_q^{\text{VQS}}\sim \mathcal{O}\left(\frac{d_{\text{ansatz}}+d_{\text{shift}}}{{\epsilon }^2}\right),$$

where d_ansatz ∼ $\mathcal{O}$(n_l) is the depth of the ansatz, d_shift ∼ $\mathcal{O}$(n²) is the depth of the shift operator, and the denominator $\mathcal{O}$(ε²) reflects the number of shots required for estimated expectation values up to a mean squared error of ε². Another consideration is the depth for amplitude encoding d_enc, which can range from $\mathcal{O}$(n²) to $\mathcal{O}$(2ⁿ). 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ⁿ × 2ⁿ, compared to a n_p × n_p matrix for VQS, suggesting comparable performance at n_l ≈ 2ⁿ/n.

3. 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 study, we focus on colloidal transport as an application of choice, as the governing Smoluchowski equation involves deep interaction potential energy wells which can be modelled as a component of the Hamiltonian operator [equation (9)], an aspect oft-neglected in quantum simulations.

3.1 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:

(20) $$\frac{\partial p\left(h,t\right)}{\partial t}=\nabla ·D\left(\nabla +\nabla U\left(h\right)\right)p\left(h,t\right),$$

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:

(21) $$\mathit{U}=Z{e}^{-\kappa H}-\frac{A}{H},$$

where U = U(h) and H = (h – a)/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:

(22) $$\begin{array}{l}Z=\frac{64\pi \epsilon {\epsilon }_{0}a{k}_{\text{B}}T}{{\left(z_{v}e\right)}^2}\tanh \left(\frac{z_{v}e{\zeta }_p}{k_{\text{B}}T}\right)\tanh \left(\frac{z_{v}e{\zeta }_w}{k_{\text{B}}T}\right),\hfill \\ \kappa =a\sqrt{\frac{2{\left(z_{v}e\right)}^{2}I}{\epsilon {\epsilon }_0{k}_{\text{B}}T}},\hfill \\ A=\frac{A_{\text{H}}}{6k_{\text{B}}T},\hfill \end{array}$$

where ε is the relative permittivity of the medium, ε₀ 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. ζ_p and ζ_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:

(23) $$\begin{array}{l}\mathit{U}\prime =-Z\kappa e^{-\kappa H}+\frac{A}{H^2},\hfill \\ {\mathit{U}}^{″}=Z{\kappa }^2{e}^{-\kappa H}-\frac{2A}{H^3}.\hfill \end{array}$$

Rescaling time τ = tD/a², we rewrite equation (20) in dimensionless form, which gives the evolution of the probability p(τ) = p(H, τ) as:

(24) $$\frac{\partial \mathit{p}\left(\tau \right)}{\partial \tau }=\frac{{\partial }^2\mathit{p}\left(\tau \right)}{\partial H^2}+\mathit{U}\prime \frac{\partial \mathit{p}\left(\tau \right)}{\partial H}+\mathit{U}″\mathit{p}\left(\tau \right),$$

which follows from equation (1), where $a$ = 1, $b$ = U′, $c$ = U″ 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(∞)/∂H → 0 (Neumann) (Torres-Díaz et al., 2019).

Substituting p(τ) = ρ(τ)e^–U/2 (Fontanela et al., 2021), we express in Dirac notation:

(25) $$\frac{\partial }{\partial \tau }|\rho \left(\tau \right)〉=\mathcal{H}|\rho \left(\tau \right)〉,|\rho \left(0\right)〉=|{\rho }_0〉,$$

with the Hamiltonian operator:

(26) $$\mathcal{H}=\frac{{\partial }^2}{\partial H^2}+{\mathit{\phi }}^T\circ \mathbb{I},$$

and the potential term:

(27) $$\mathit{\phi }=\frac{{\mathit{U}}^{″}}{2}-{\left(\frac{\mathit{U}\prime }{2}\right)}^2,$$

which can be evaluated classically [equation (23)] and implemented as a quantum observable φ^T_° $\mathbb{I}$. With Dirichlet–Neumann boundary conditions enforced, the Hamiltonian operator is decomposed as:

(28) $$\begin{array}{ll}\mathcal{H}\hfill & ={\left[\mathit{\phi }-\frac{1}{\Delta x^2}(1-e_{2^n})\right]}^T\circ \underset{H_1}{\underbrace{I^{\otimes n}}}\hspace{0.17em}+\hspace{0.17em}\frac{1}{\Delta x^2}\{\underset{H_2}{\underbrace{I^{\otimes n-1}\otimes X}}\hspace{0.17em}+\hspace{0.17em}{S}^{†}[-\underset{H_3}{\underbrace{I^{\otimes n}}}+I^{\otimes n-1}\otimes X︸H_4-I_0^{\otimes n-1}\otimes X︸H_5\left]S\right\}\hfill \end{array},$$

where:

$$e_{2^n}=(0,0,\dots ,0,1).$$

3.1.1 Potential-free case (φ = 0).

Consider first the potential-free case where colloid–wall interactions are absent (φ = 0). The probability density state evolves in space H ∈ [0, 1] and time τ ∈ [0, T] as:

(29) $$\frac{\partial }{\partial \tau }|\rho (\tau )〉=\frac{{\partial }^2}{\partial H^2}|\rho \left(\tau \right)〉,\left|\rho \right(0)〉=|{\rho }_0〉,$$

where the initial impulse state is centred at |ρ₀〉 = |2^n–1〉. Using a full circular ansatz with six to eight layers, we evolve the initial pulse on a 2ⁿ = 16 grid using time-step Δτ = 10⁻⁴ for early times up to 10⁻² [Figure 4(a)] and late times up to T = 10⁻¹ [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⁻¹. 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⁻⁴; it is reduced only with further increase in the number of ansatz layers n_l > 6, leading to optimal scaling of ${\overline{\epsilon }}_{\text{tr}}\sim N_q^{-2}$. 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 \left({\epsilon }_{\text{tr}}\right)\sim N_q^{-0.8}$.

Using Euler time-stepping, the mean norm error scales as ${\overline{\epsilon }}_{\text{norm}}\sim N_q^{-0.8}$, regardless of n_l [Figure Figure 4(d)]. This cost scaling improves significantly up to ${\overline{\epsilon }}_{\text{norm}}\sim N_q^{-3.4}$ 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 \left({\epsilon }_{\text{norm}}\right)\sim N_q^{-0.8}$, regardless of time-stepping scheme.

3.1.2 DLVO potential φ(A, Z, κ).

In the presence of colloid–wall interactions, the DLVO potential term φ depends minimally on three parameters, specifically A, Z and κ [equation (22)]. Following the potential-free case (n = 4, H ∈ [0, 1], Δτ = 10⁻⁴), we perform VQS including φ using eight ansatz layers in time τ ∈ [0, T].

In the absence of the electric double layer (Z = 0), the DLVO potential φ(A) depends on only the van der Waal’s interaction energy, assumed here to be attractive. Figure 5(a) shows that the DLVO potential φ(H) profiles scaled by the square of the interval ΔH² for A ∈ {0.05, 0.1, 0.2, 0.5} is only short-ranged in H, so the quantum solution |ρ〉 is insensitive to A. Recall, however, the earlier substitution p(τ) = ρ(τ)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, τ) for A = 0.05 up to T = 0.1 is depleted near wall [Figure 5(c)] compared to the potential-free (φ = 0) case [Figure 5(b)]. Increasing A further increases the depletion range [Figure 5(d)].

Otherwise, the DLVO potential φ(A, 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 κ. 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, τ) for {Z, κ} = {10, 10} shows long-ranged influence of the electric double-layer interaction. Parametric analyses of {Z, κ} holding A = 0.5 show that Z depletes p(H, τ) near wall [Figure 6(c)], and a decrease in κ increases the deposition flux and depletion range [Figure 6(d)].

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(τ) [equation (15)] and the norm error ε_norm(τ) [equation (16)]. Figure 7 shows that ε_trace(τ) peaks and decreases during the early diffusion phase [Figure 4(a)], then peaks and decreases again as the normalized probability density $\widehat{\rho }$ 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(τ) [Figure 7(b)]. In contrast, ε_norm(τ) tends towards a steady state regardless of the evolution of probability density. Parametric analyses suggest that ε_norm is affected by the local depletion of $\widehat{\rho }$ but insensitive to the magnitudes of A [Figure 7(a)] and Z [Figure 7(b)].

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).

3.2 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):

(30) $$\frac{\partial c\left(h,t\right)}{\partial t}=\nabla ·D\left(\nabla +\nabla \mathit{U}\left(h\right)\right)c\left(h,t\right),$$

which describes the concentration of colloidal particles c(h, t) instead of probability density but is otherwise identical to the Smoluchowski equation [equation (20)]. The difference here is that a continuous concentration source can be imposed as a far-field Dirichlet boundary condition. Rescaling c(H, τ) in space H ∈ [0, 1] and time τ ∈ [0, T], we perform a change of variables c(τ) = ς(τ)e^–U/2 as before, and write:

(31) $$\frac{\partial }{\partial \tau }|\varsigma \left(\tau \right)〉=\mathcal{H}|\varsigma \left(\tau \right)〉+\mathcal{F}|0〉,|\varsigma \left(0\right)〉=|{\varsigma }_0〉,$$

where the operator $\mathcal{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 $\mathcal{F}$.

3.2.1 Initialization.

We seek a parameterized ansatz that encodes a Heaviside step function centred at |2^n–1〉:

(32) $$|{\varsigma }_0〉=\left(X\otimes H^{\otimes n-1}\right)|0〉.$$

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 ${\theta }_{[1,n]}^{n_l}=\frac{\pi }{2}$ and the second R_Y of the preceding layer as ${\theta }_2^{n_l-1}=\frac{{(-1)}^{n-1}\pi }{2}$, where a reversal in the sign produces a step-down function instead.

3.2.2 Solutions and errors.

We perform VQS on a 2ⁿ = 16 grid using time-step Δτ = 10⁻⁴ as before, but on a full circular ansatz with five layers, which is already shown to yield high-fidelity solutions [Figure 2(c) and 2(d)]. Figure 8(a) shows how the normalized concentration $\widehat{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 κ increases the depletion range [Figure 8(d)]. Both trace and norm errors (Figure 8 insets) decay in time τ towards convergence with ε_trace peaking earlier than ε_norm.

4. 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)]:

• {∂_tu(t), ∂_xxu(t)} terms handled using variational quantum imaginary time evolution.

• ∂_xu(t) term eliminated through substitution methods, such as u(t) = e^gv(t).

• u(t) term included in the Hamiltonian $\mathcal{H}$ without additional complexity cost.

• f(t) term realized by an additional set of complementary circuits, whose complexity depends on $\mathcal{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 (Kubo et al., 2021; 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 |x〉, 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 over-parameterization 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 φ(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 (Nakaji et al., 2022; Zoufal et al., 2020).

Figures

Figure 1.

Real-amplitude ansatz formed by repeating parameterized blocks with either linear (a, b) or circular (c, d) entanglement. (a, c): Final layer n_l is unentangled. (d): This full circular ansatz outperforms the other ansätze and is used throughout the study

Figure 2.

Initial step evolves under (a) periodic and (b) Dirichlet boundary condition on a 2ⁿ = 16 grid for a real-amplitude ansatz with four layers using time-step Δt = 10⁻⁴ plotted in increments of 2 × 10⁻³; mean (c) trace and (d) norm error plotted on a semi-log scale against the number of ansatz layers for various real-amplitude designs under periodic (open symbols) and Dirichlet (closed symbols) boundary conditions up to T = 10⁻² (insets show peak error)

Figure 3.

Mapping state |x〉 to an input state |x₀〉 at the n_l – k layer of a four-qubit full circular ansatz

Figure 4.

Normalized colloidal probability density $\widehat{\rho }\left(H\right)$ profiles without interaction (U = 0) on a 2ⁿ = 16 grid using time-step Δτ = 10⁻⁴ plotted in increments of (a) 2 × 10⁻³ and (b) 2 × 10⁻² up to T = 10⁻¹. Lines show VQS solutions based on full circular ansatz with six to eight layers, and circles show classical solutions based on the same discretization. Mean (c) trace and (d) norm errors plotted on log scale versus the total number of circuits N_q for run-time up to T = 10⁻¹ (insets show peak values). Data are grouped by the number of ansatz layers, each set using time-step Δτ ∈{10, 5, 2, 1} × 10⁻⁴. Closed symbols represent Runge–Kutta solutions using time-step Δτ ∈{10, 5, 2} × 10⁻⁴. Horizontal reference line: 0.05

Figure 5.

(a) DLVO potential φ(H) profiles scaled by ΔH² using a full circular ansatz with eight layers for Z = 0 varying A. Inset shows interaction energy U(H). (b–d) Space-time plots of colloidal probability density p(H, τ) ∈ [Δp, 0.2] with contour interval Δp = 0.01. Given Z = 0, (b) φ = 0, (c) A = 0.05 and (d) A = 0.5

Figure 6.

(a) DLVO potential φ(H) profiles scaled by ΔH² using a full circular ansatz with eight layers with A = 0.5 varying Z and κ. Inset shows interaction energy U(H). (b–d) Space-time plots of colloidal probability density p(H, τ) ∈ [Δp, 0.2] with contour interval Δp = 0.01. Given A = 0.5, {Z, κ} are (b) {10, 10}, (c) {20, 10} and (d) {10, 5}

Figure 7.

Trace error ε_trace and norm error ε_norm in time τ for (a) Z = 0, varying A, and (b) A = 0.5, varying Z and κ. Horizontal reference line: 0.05

Figure 8.

Normalized colloidal concentration $\widehat{c}\left(H\right)$ profiles on a 2ⁿ = 16 grid using time-step Δτ = 10⁻⁴ plotted in increments of 10⁻² up to T = 10⁻¹ (bold red). Lines show VQS solutions based on the full circular ansatz with five layers. (a) U = 0 and (b) A = 0.05, Z = 0. Given A = 0.05, (c) Z = 20, κ = 10 and (d) Z = 10, κ = 5.

Figure A1.

Quantum circuits to evaluate (a) $\Re \left(\frac{〈\tilde{u}\left(\mathbf{\theta }\right)|}{\partial {\theta }_i}\frac{\left|\tilde{u}\right(\mathbf{\theta })〉}{\partial {\theta }_j}\right)$ via Hadamard tests, (b) $\Re \left(\frac{〈\tilde{u}\left(\mathbf{\theta }\right)|}{\partial {\theta }_i}\mathcal{H}\right|\tilde{u}\left(\mathbf{\theta }\right)〉)$ and (c) $\Re \left(\frac{〈\tilde{u}\left(\mathbf{\theta }\right)|}{\partial {\theta }_i}\mathcal{F}\right|0〉)$ (Miyamoto and Kubo, 2021) via direct measurements with respect to observables $\mathcal{H}$ and $\mathcal{F}$, respectively (Zoufal et al., 2021)

Figure A2.

Directed graphs (digraphs) depicting the disjoint orbits arising from applying the C_n matrix to the computational basis states

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.

Appendix 1. Quantum circuits to evaluate A and C

The elements of matrix A [equation (6)] and vector C [equation (7)] can be evaluated via sampling the expectation of an observable Z using quantum circuits shown in Figure A1 (Zoufal et al., 2021). The derivative of the trial state $\left|\tilde{u}\right(\mathit{\theta }\left(t\right))〉$ with respect to θ_k is:

(33) $$\frac{\partial \left|\tilde{u}\right(\mathit{\theta }\left(t\right))〉}{\partial {\theta }_k}=f_k\left[{\theta }_0\left(t\right)R_1\left({\theta }_1\right(t\left)\right)R_2\left({\theta }_2\right(t\left)\right)\cdots {\sigma }_k{R}_k\left({\theta }_k\right)\cdots R_N\left({\theta }_N\right(t\left)\right)\right]|0〉,$$

such that for a single-qubit rotation gate $R_Y\left({\theta }_k\right)=e^{-i{\theta }_k{\sigma }_Y/2}$, the gate derivative $\partial R_k\left({\theta }_k\right)/\partial {\theta }_k=-(i/2){\sigma }_Y{e}^{-i{\theta }_k{\sigma }_Y/2}$ is measurable with coefficient f_k = –i/2 and a Pauli-Y gate σ_k = σ_Y inserted in the trial state. Accordingly, the quantum circuit may incur a global phase e^–iα, where α = 0 and π/2 for A and C, respectively, which may be rectified through an additional phase gate [2], $S={\sigma }_Z^{1/2}$ on the ancilla qubit.

We implemented Hadamard tests [Figure 1(a)] in IBM Qiskit using the aer_simulator backend with sampling count of 2¹² shots per circuit evaluation and direct measurements using statevector_simulator with respect to observables $\mathcal{H}$ [Figure 1(b)] and $\mathcal{F}$ [Figure 1(c)]. Note that the latter requires a controlled trial state formed by controlled parameterized unitaries (Miyamoto and Kubo, 2021).

Appendix 2. On the encoding of bit strings using full circular ansatz

In this appendix, we elaborate on the initialization procedure described in Section 2.4 that uses the full circular ansatz of Figure 1(d). Firstly, we note that the ansatz is diagonal in the computational basis and therefore preserves computational basis states. Hence, for the rest of this analysis, it suffices to just consider the action of the ansatz on bit strings. As a function C_n: {0, 1}ⁿ → {0, 1}ⁿ, the ansatz transforms bit strings as follows in little-endian:

(34) $$C_n(x_{n-1},\dots ,x_0)=\left({\displaystyle \sum _{i=0}^{n-1}{x}_i},{\displaystyle \sum _{i=0}^{n-2}{x}_i},{\displaystyle \sum _{i=0}^{n-3}{x}_i},\dots ,{\displaystyle \sum _{i=0}^1{x}_i},{\displaystyle \sum _{i=1}^{n-1}{x}_i}\right).$$

For example, when n = 6, one can check that:

(35) $$\begin{array}{l}{C}_6(x_5,x_4,x_3,x_2,x_1,x_0)=(x_0+x_1+x_2+x_3+x_4+x_5,x_0+x_1+x_2+x_3+x_4,x_0+x_1+x_2+x_3,x_0+x_1+x_2,x_0+x_1,x_1+x_2+x_3+x_4+x_5).\end{array}$$

Hence, C_n can be represented as the n × n matrix whose (i, j)-th entry is (C_n)_ij = 0 if j > i > 0 or i = j =0, and 1 otherwise. For example, when n = 6, this matrix is given by:

(36) $$C_6=\left[\begin{array}{cccccc}0& 1& 1& 1& 1& 1\\ 1& 1& 0& 0& 0& 0\\ 1& 1& 1& 0& 0& 0\\ 1& 1& 1& 1& 0& 0\\ 1& 1& 1& 1& 1& 0\\ 1& 1& 1& 1& 1& 1\end{array}\right].$$

Each layer of entangling CNOT gates (Figure 1) corresponds to an application of the C_n matrix to the input bit string, denoted |x₀〉. Successive application of C_n generates a sequence of bit strings $C_n|{\mathit{x}}_0〉,C_n^2|{\mathit{x}}_0〉,\dots ,C_n^k|{\mathit{x}}_0〉$, eventually resulting in the initial bit string for some period, p, i.e. $C_n^p=I_n$, thus forming a periodic sequence or orbit. We illustrate the sequences of bit strings generated by this procedure for the cases of n = 4, 5, 6 in Figure A2. We note that p ≤ 2ⁿ – 1, and in the cases where p < 2ⁿ – 1, there are disjoint orbits corresponding to the different irreducible representations of the full circular ansatz group. The trivial one-dimensional irrep corresponds to the all-zeroes bit string |00…0〉 set by θ = 0. For n = 4, the strings form two distinct orbits, whereas, for n = 5 and 6, the strings form four distinct orbits (see Figure A2). For n = 4 and 6, all the orbits, save the singleton orbit comprising the all-zeroes bit string, contains at least one string of positive Hamming weight, and hence all 2ⁿ – 1 states of positive Hamming weight are reachable from the strings of unit Hamming weight. This is not the case for n = 5, where all the strings of unit Hamming weight are in the same orbit (namely, the orbit of size 21); hence, only these 21 strings out of the 2⁵ – 1 = 31 strings of positive Hamming weight are reachable from the strings of unit Hamming weight.