We introduce an efficient variational hybrid quantum-classical algorithm designed for solving Caputo time-fractional partial differential equations. Our method employs an iterable cost function incorporating a linear combination of overlap history states. The proposed algorithm is not only efficient in terms of time complexity but also has lower memory costs compared to classical methods. Our results indicate that solution fidelity is insensitive to the fractional index and that gradient evaluation costs scale economically with the number of time steps. As a proof of concept, we apply our algorithm to solve a range of fractional partial differential equations commonly encountered in engineering applications, such as the subdiffusion equation, the nonlinear Burgers' equation, and a coupled diffusive epidemic model. We assess quantum hardware performance under realistic noise conditions, further validating the practical utility of our algorithm.

Differential equations play an important role in the modeling of many practical engineering problems. Fractional differential equations, which are generalizations of differential equations to an arbitrary noninteger order, involve fractional derivatives such as d α / d x α, where α is a noninteger real number. These models have been shown to be capable of modeling more complex processes that exhibit memory effects or nonlocal behavior. For example, fractional differential equations have found applications in fractional advection–dispersion equations1 and viscoelastic flow problems.2 Fractional-order derivatives in time have also been used to model cancer tumor growth, by numerically fitting the fractional order.3 

Here, we consider numerical solutions to time-fractional diffusive differential equations of order α, which derive from non-Markovian continuous-time random walks.4 Finite difference methods are effective and simple to implement for time-fractional problems,5 where implicit methods are preferred to explicit methods for numerical stability.6 However, the time-fractional derivative creates a global dependence problem7 where previous solutions are accessed from memory, rendering storage expensive for spatially large problems.

Quantum computers can offer new tools for efficiently solving differential equations. In particular, variational quantum algorithms (VQAs) have emerged as a leading strategy in the noisy intermediate-scale quantum (NISQ) era,8,9 featuring low-depth parameterized quantum circuits alongside classical optimizers used to minimize a cost function. These VQAs yield quantum states that encode approximate solutions to partial differential equations. They have been successfully applied to both linear10–12 and nonlinear differential equations13,14 in a diverse range of fields, including electromagnetics,15 colloidal transport,16 fluid dynamics,17–19 and finance.20 

In this paper, we propose and implement a variational quantum optimization scheme for solving time-fractional diffusive partial differential equations in space and time, either implicitly or semi-implicitly, using simple quantum overlap measurements for the time integral term. We show that our approach is not only efficient in time complexity but can also yield memory cost advantages in handling the global dependence problem. The rest of the paper is organized as follows: In Sec. II, we introduce the time-fractional diffusion equation and present our variational quantum optimization scheme for numerical integration. In Sec. III, we discuss the details of the implementation and evaluate time and memory complexities of our scheme. In Sec. IV, we present numerical experiments, including time-fractional Burgers' equation for hydrodynamics (Sec. IV B), and fractional diffusive epidemic model (Sec. IV C). Finally, we trial our algorithm on noise model simulators and hardware in Sec. V, and conclude with discussions in Sec. VI.

We consider the fractional diffusion equation6,21 in space x and time t,
D t α u ( t , x ) = x x u ( t , x ) ,
(1)
where D t α is the Caputo fractional derivative of the α th order with respect to t, defined as
D α u ( t , x ) t α = 1 Γ ( n α ) 0 t n u ( s , x ) s n d s ( t s ) α n + 1 , n 1 < α < n ,
(2)
where Γ denotes the Gamma function and n = α is the ceiling of α.

We set up a space–time rectangular domain with spatial extent L and temporal duration T, dividing it into a regularly spaced grid with integer size N × M. The spatial grid points are defined as x i = i h where h = L / N and i = 1 , 2 , , N, and the temporal grid points as t k = k τ where τ = T / M and k = 1 , 2 , , M.

The first-order finite difference of the Caputo derivative of order ( 0 < α 1 ) is approximated by5,22
D t α u ( t k , x n ) g α , τ j = 1 k w j ( α ) ( u n k j + 1 u n k j ) ,
(3)
where u n k u ( t k , x n ), g α , τ : = τ α / Γ ( 2 α ), and w j ( α ) : = j 1 α ( j 1 ) 1 α. Rearranging Eq. (3) yields a discrete sum of terms u n j expressed as
D t α u ( t k , x n ) g α , τ [ u n k w k u n 0 + j = 1 k 1 ( w j + 1 w j ) u n k j ] ,
(4)
where w 1 = 1 by definition.
The finite difference of the spatial derivative is approximated by
x x u ( t k , x n ) u n + 1 k 2 u n k + u n 1 k h 2 .
(5)
Together, the difference scheme for the fractional diffusion equation can be iterated in time as
( 1 a δ ) u n k = w k u n 0 j = 1 k 1 ( w j + 1 w j ) u n k j ,
(6)
where a = g α , τ 1 h 2 and δ u n : = u n + 1 2 u n + u n 1. In a more compact matrix form, this reads
A u k = w k u 0 j = 1 k 1 Δ w j u k j ,
(7)
where u k : = u ( t k , x ), Δ w j : = w j + 1 w j, and A is an N × N symmetric tridiagonal matrix with b = 1 + 2 a on the main diagonal, a on the adjacent off-diagonals, and terms in the corners of the matrix that are specified by the boundary conditions. More precisely, the matrix A is given by
A = [ c a 0 d a b a 0 0 0 a b a 0 0 0 a b a d 0 a c ] N × N ,
(8)
where ( c , d ) = ( b , a ) for the periodic boundary condition, ( c , d ) = ( b , 0 ) for the Dirichlet boundary condition, and ( c , d ) = ( 1 + a , 0 ) for the Neumann boundary condition. The right-hand side of Eq. (7) carries the history of solutions up to k 1.
Using Dirac notation, we rewrite the iterative fractional diffusion equation as
A | u ̃ k = | f ̃ k 1 ,
(9)
where we approximate the unnormalized quantum state | u ̃ k at time k by an ansatz or trial state formed by a set of parameterized unitaries R and unparameterized entangling unitaries W as
| u ̃ k : = r k | u k = r k i = 1 l W l i + 1 R ( θ l i + 1 k ) | 0 ,
(10)
where r k is the norm of the unnormalized state | u ̃ k , θ is a set of gate parameters, and l is the number of layers in the ansatz. The unnormalized source state | f ̃ k 1 is a linear combination of history states in [ 0 , k 1 ].
Following Ref. 11, we define a cost function based on potential energy at time k,
C k ( r k , θ k ) = 1 2 ( r k ) 2 u k | A | u k 1 2 r k ( u k | f ̃ k 1 + f ̃ k 1 | u k ) ,
(11)
where
| f ̃ k 1 = w k r 0 | u 0 j = 1 k 1 Δ w j r k j | u k j .
(12)
With that, u k | f ̃ k 1 expands into a list of k measurable overlap terms u k | u [ 0 , k 1 ] .23 
Introducing | f , u : = 2 1 / 2 ( | 0 | f + | 1 | u ), we have
C k ( r k , θ k ) = 1 2 ( r k ) 2 u k | A | u k r k ( w k r 0 u k , u 0 j = 1 k 1 Δ w j r k j u k , u k j ) .
(13)
Here, f , u : = f , u | X I n | f , u is used as a notation shorthand, where X = | 1 0 | + | 0 1 | is the Pauli-X matrix and I : = I 0 + I 1 = | 0 0 | + | 1 1 | is the 2 × 2 identity matrix. Also, n = log 2 N denotes the number of qubits representing N spatial grid points.
Taking the derivative of C k with respect to r k, the optimal norm
r k = u k , f ̃ k 1 u k | A | u k
(14)
can be eliminated from the cost function
C k ( r k , θ k ) = 1 2 ( u k , f ̃ k 1 ) 2 u k | A | u k ,
(15)
which can be measured using quantum circuits [Figs. 12(a) and 12(b)] and optimized classically as
arg min θ k C k ( r k , θ k ) ,
(16)
via either gradient-based (evaluate θ k C k) or gradient-free methods. A schematic of the proposed variational quantum algorithm is shown in Fig. 1.
Fig. 1.

Schematic diagram of the variational quantum algorithm for time-fractional differential equations. The initial solution, encoded by the vector θ 0 containing nl optimized parameters and the norm r 0, is stored in classical memory. Subsequent iterations add θ k and r k to the growing pool of classically stored parameters θ [ 0 , k 1 ] and norms r [ 0 , k 1 ], until the final time step M is reached. Solutions can be read out selectively via tomography using up to O ( n 2 n ) measurements on specific θ k parameterized ansatz corresponding to time k. Sampling complexity can be reduced using partial tomography or shadow tomography techniques.

Fig. 1.

Schematic diagram of the variational quantum algorithm for time-fractional differential equations. The initial solution, encoded by the vector θ 0 containing nl optimized parameters and the norm r 0, is stored in classical memory. Subsequent iterations add θ k and r k to the growing pool of classically stored parameters θ [ 0 , k 1 ] and norms r [ 0 , k 1 ], until the final time step M is reached. Solutions can be read out selectively via tomography using up to O ( n 2 n ) measurements on specific θ k parameterized ansatz corresponding to time k. Sampling complexity can be reduced using partial tomography or shadow tomography techniques.

Close modal
The Hamiltonian matrix A [Eq. (8)] can be decomposed into
A = b I n a { I n 1 X H 1 + S [ I n 1 X H 2 I 0 n 1 X H 3 + I 0 n 1 I H 4 ] S } .
(17)
Since expectation values of the identity operator I are equal to 1, i.e., ϕ | I n | ϕ = ϕ | I n | ϕ = 1,11 evaluating the expectation value of A requires only the evaluation of the expectation values of the simple terms ( H 1 2 for the periodic boundary condition, H 1 3 for the Dirichlet boundary condition, and H 1 4 for the Neumann boundary condition). The operator S denotes the n-qubit cyclic shift operator,
S = i = 0 2 n 1 | ( i + 1 ) mod 2 n i | ,
(18)
which can be implemented using relative-phase Toffoli gates, Toffoli gates, a cnot gate and an X gate.24 
For problems that admit only real solutions, it is preferable to use a real-amplitude ansatz, represented by Eq. (10) formed by l repeating blocks, each consisting of a parameterized R layer with one R Y ( θ ) gate on each qubit, followed by an unparameterized W layer with a cnot gate between consecutive qubits,25,
| u k = i = l 1 [ a = n 1 1 C a X a + 1 j = 1 n R Y ( θ i , j ) ] | 0 ,
(19)
where we have employed the product notation with an initial index greater than the stopping index to signify that the terms in the product are arranged in decreasing index order.26 In addition, we consider a real-amplitude ansatz with circular entanglement, such that
| u k = i = l 1 [ C n X 0 a = n 1 1 C a X a + 1 j = 1 n R Y ( θ i , j ) ] | 0 , n > 2 ,
(20)
which has been shown to be efficient for solving partial differential equations.16 See Fig. 11 for circuit diagrams of the ansätze.
To solve Eq. (9), we prepare the initial ansatz state | u ̃ 0 = r 0 | u ( θ 0 ) using classical optimization on parameters θ = ( θ i , j ) ( i , j ) { 1 , , l } × { 1 , , n } such that
θ 0 arg min θ 0 u ̂ 0 | u ( θ 0 ) ,
(21)
where | u ̂ 0 u 0 / u 0 is the normalized initial classical solution vector and r 0 = u 0 is the initial norm.
The gradient of the cost function Eq. (15) can be estimated on quantum computers using the parameter shift rule. Following Ref. 11, here we write down the partial derivative of the cost function with respect to parameter θ i , j k indexed by ( i , j ) [ 1 , n l ] at k as
C k θ i k = 1 2 ( u k , f ̃ k 1 u k | A | u k ) [ w k r 0 u k θ i k , u 0 j = 1 k 1 Δ w j r k j u k θ i k , u k j ] + 1 2 ( u k , f ̃ k 1 u k | A | u k ) 2 u k θ i k , u k A ,
(22)
where f , u A : = f , u | X A | f , u is used as a notational shorthand. The parametric shift is implemented by a π rotation on the ith R Y gate as
u k θ i k = u k ( θ i % n i / n + π ) ,
(23)
where i % n is the remainder r { 0 , 1 , , n 1 } obtained when i is divided by n and i / n is the integer part of i / n. Thus, for each index i at time k, gradient estimation requires at least k overlap measurements.

Time complexity: Here, we briefly estimate the time complexity of the algorithm based on the number of time steps, quantum circuits, gates, and shots required, neglecting classical overheads such as parameter update, iteration, functional evaluation, and initial encoding.

For a number of time steps M, the overall time complexity is
T M N it N g [ M l + ( l + n 2 ) ϵ 2 ] ,
(24)
where n is the number of qubits, l is the number of ansatz layers, N it is the mean number of iterations required per time step, and N g is the mean number of evaluations required for gradient estimation per iteration. Specifically, the terms in square brackets are the gate complexities for the fractional derivative Ml based on M / 2 overlap circuits containing 2l parametric gates in depth and for the Hamiltonian ( l + n 2 ) based on O ( 1 ) circuits [ 2 4 Eq. (17)] containing l parametric gates in depth for the ansatz and O ( n 2 ) nonparametric gates for the shift operator. Each circuit is sampled repeatedly for O ( ϵ 2 ) shots per measurement, where ϵ is the desired precision. We remark that with quantum amplitude estimation,27 the number of shots required may be quadratically reduced from O ( ϵ 2 ) to O ( ϵ 1 ). This decrease, however, entails increased circuit depth, as exemplified by the Ω ( ϵ 1 ) circuit depths utilized in Refs. 27–29. Given these challenges in reducing circuit depth in quantum amplitude estimation and its extensions,30–33 we do not consider this approach in the present study, deferring it for future exploration.
For n l, the gradient-based optimizer N g O ( n l ) yields an overall time complexity
T O ( max { N it M 2 n l 2 ϵ 2 , N it M n 3 l ϵ 2 } )
(25)
scaling as O ( M 2 ) for overlap circuits, or O ( M ) for shift operators, the latter scaling applicable to the nonfractional diffusion/heat equation ( α = 1,11). Using the simultaneous perturbation stochastic approximation (SPSA),34 only two circuit executions are required per gradient evaluation N g O ( 1 ), independent of the number of parameters, similar to gradient-free optimization.

The present algorithm is efficient in time complexity with respect to the number of spatial grid points N ( = 2 n ),11,35 neglecting initial encoding costs. Quantum advantage with respect to time T may be achievable in higher dimensions (Remark 2,35 Remark 336).

Memory scaling: To solve the Caputo derivative Eq. (4), a classical algorithm assesses O ( k N ) parameters from memory at time k, compared to O ( knl ) parameters for variational quantum algorithm Eq. (12). Assuming n l, the memory cost of the Caputo sums scales subexponentially as O ( k log N ) (Fig. 1), demonstrating efficient algorithmic space complexity.

The circuits are implemented in the quantum software framework Pennylane37 using noiseless statevector emulation on lightning qubit backend. Parametric updates are performed using the limited-memory Broyden–Fletcher–Goldfarb–Shanno boxed (L-BFGS-B) optimizer with a relative tolerance of 10 6 and a gradient tolerance of 10 6. The L-BFGS-B optimizer evaluates gradients using finite difference without explicit gradient inputs from Eq. (22). Initial encoding Eq. (21) is performed using SciPy's minimize function with an initial null parameter set as ( 0 , , 0 ).

We solve for the time-fractional ( 0 < α 1 ) subdiffusion equation Eq. (1) as an initial value problem with an initial condition u ( 0 , x ) = x ( 1 x ) and a Dirichlet boundary condition u ( t , 0 ) = u ( t , 1 ) = 0.6 

Figure 2 shows that the variational quantum solutions agree with classical subdiffusion solutions in the 32 × 32 space–time domain with qubit-layer count ( n , l ) = ( 5 , 4 ), for (a) α = 1 and (b) α = 0.5. Figure 3 shows that low fractional-order solutions tend to diffuse faster for short times, but slow down for long times, both characteristics of subdiffusive behavior.

Fig. 2.

Variational quantum solutions to time-fractional diffusion equation in N × M = 32 × 32 space–time domain with qubit-layer count ( n , l ) = ( 5 , 4 ), for (a) α = 1 and (b) α = 0.5. Initial and boundary conditions are u ( 0 , x ) = x ( 1 x ) and u ( t , 0 ) = u ( t , 1 ) = 0. The colorbar denotes a relative deviation from the classical finite difference solution on the same domain.

Fig. 2.

Variational quantum solutions to time-fractional diffusion equation in N × M = 32 × 32 space–time domain with qubit-layer count ( n , l ) = ( 5 , 4 ), for (a) α = 1 and (b) α = 0.5. Initial and boundary conditions are u ( 0 , x ) = x ( 1 x ) and u ( t , 0 ) = u ( t , 1 ) = 0. The colorbar denotes a relative deviation from the classical finite difference solution on the same domain.

Close modal
Fig. 3.

Variational quantum solutions (as Fig. 2) plotted on (a) α x at t = 0.5 and (b) t α at x = 0.5. The colorbar denotes a relative deviation from the classical finite difference solution on the same domain.

Fig. 3.

Variational quantum solutions (as Fig. 2) plotted on (a) α x at t = 0.5 and (b) t α at x = 0.5. The colorbar denotes a relative deviation from the classical finite difference solution on the same domain.

Close modal
The trace error at time k,
ϵ tr k = 1 u ̂ k | u ( θ k ) 2 ,
(26)
can be averaged over the number of time steps M for 0 < α 1, where | u ̂ k u k / u k is the normalized classical solution vector at time k.

1. Trace error and optimization

Figure 4(a) compares the results between the linear real-amplitude hardware-efficient ansatz Eq. (19) and the circular ansatz Eq. (20), suggesting an improved solution fidelity for the circular ansatz over the linear ansatz ( n > 4). According to the unitary dependence theory,38 the effect of an additional circular entangling cnot gate ( C n X 0) on a single layer ansatz ( l = 1) is to transform the unitary dependence of the first qubit from q 1 : { R Y ( θ 1 ) } to q 1 : { R Y ( θ 2 ) R Y ( θ n ) }, which effectively increases the connectivity of q 1 by n 2 qubits. This may improve the performance of an ansatz for certain problems.38 

Fig. 4.

(a) Time-averaged trace error ϵ tr plots for linear positive-amplitude hardware-efficient ansatz Eq. (19) (dotted lines) vs circular ansatz Eq. (20) (continuous lines), for 0 < α 1. Error bars denote standard deviation. (b) Number of function evaluations k = 1 M N eval k required by L-BFGS-B gradient-free optimizer for α = 0.5 (dashed lines) and α = 1 (continuous lines) scales with number of time steps M to a fractional power of 0.6 in time t = [ 0 , 0.5 ]. Indicated ( n , l ) pairs are selected based on the minimum number of layers l required to encode | u ( θ 0 ) for n { 3 , 4 , 5 }.

Fig. 4.

(a) Time-averaged trace error ϵ tr plots for linear positive-amplitude hardware-efficient ansatz Eq. (19) (dotted lines) vs circular ansatz Eq. (20) (continuous lines), for 0 < α 1. Error bars denote standard deviation. (b) Number of function evaluations k = 1 M N eval k required by L-BFGS-B gradient-free optimizer for α = 0.5 (dashed lines) and α = 1 (continuous lines) scales with number of time steps M to a fractional power of 0.6 in time t = [ 0 , 0.5 ]. Indicated ( n , l ) pairs are selected based on the minimum number of layers l required to encode | u ( θ 0 ) for n { 3 , 4 , 5 }.

Close modal

Figure 4(b) shows that the number of function evaluations k = 1 M N eval k, required by the L-BFGS-B gradient-free optimizer scales sublinearly with the number of time steps O ( M 0.6 ) in time t = [ 0 , 0.5 ]. The overhead costs for classical gradient estimation remains economical up to M = 2 6.

Figure 5 shows that increasing the number of parameters nl does not improve the time-averaged trace error ϵ tr, but instead increases the number of function evaluations which scales as O ( n l ). The linear scaling suggests that the problem of vanishing gradients due to excessive parameterization, commonly known as barren plateaus,39 may be mitigated by short time steps, such that the initial solution at each time step is closer to the optimal solution.

Fig. 5.

(a) Time-averaged trace error ϵ tr is insensitive to over-parameterization, shown here for α = 0.5 (dashed lines) and α = 1 (continuous lines). Error bars denote standard deviation. (b) Number of function evaluations k = 1 M N eval k scales linearly with number of parameters nl.

Fig. 5.

(a) Time-averaged trace error ϵ tr is insensitive to over-parameterization, shown here for α = 0.5 (dashed lines) and α = 1 (continuous lines). Error bars denote standard deviation. (b) Number of function evaluations k = 1 M N eval k scales linearly with number of parameters nl.

Close modal
The fractional Burgers' equation40,41 has found applications in shallow water waves and waves in bubbly liquids.42 We consider the following 1D time-fractional Burgers' equation43 as an extension to the time-fractional diffusion equation Eq. (1),
D t α u ( t , x ) = ν x x u ( t , x ) u ( t , x ) x u ( t , x ) , 0 < α 1 ,
(27)
where D t α is the Caputo fractional derivative of the α th order, as previously defined, and ν is the fluid kinematic viscosity. Readers interested in quantum computing algorithms for nonfractional Burgers' equation ( α = 1) are directed to Refs. 13, 44, and 45.

1. Semi-implicit scheme

Using explicit source terms, we apply central difference on the nonlinear advection term, such that
( 1 a δ ) u n k = w k u n 0 j = 1 k 1 Δ w j u n k j b ( u n + 1 k 1 u n 1 k 1 ) u n k 1 ,
(28)
where a = ν / ( g α , τ h 2 ) and b = 1 / ( 2 g α , τ h ). In Dirac notation, the time-fractional Burgers' equation can be iterated as in Eq. (9) and expressed as
A | u ̃ k = | f ̃ k 1 b ( Λ ̃ + k 1 Λ ̃ k 1 ) | u ̃ k 1 ,
(29)
where Λ ̃ : = diag ( u ) is a diagonal matrix with u along its diagonal, and the subscript ± denotes either an incremental or a decremental shift in each eigenstate under periodic boundary conditions, i.e., | 2 n = | 0 .
With that, the cost function for Burgers' equation is
C k ( r k , θ k ) = 1 2 [ u k , f ̃ k 1 b ( u k , Λ ̃ + k 1 u ̃ k 1 u k , Λ ̃ k 1 u ̃ k 1 ) ] 2 u k | A | u k ,
(30)
which includes two nonlinear terms u k , Λ ̃ ± k 1 u ̃ k 1 evaluated using redundant copies of k 1 quantum states13,14 on a quantum circuit [Fig. 12(c)].

2. Test case

We test a 1D bi-directional flow using time-fractional Burgers' equation Eq. (27) with an initial condition u ( 0 , x ) = sin ( 2 π x ) and Dirichlet boundary conditions u ( 0 , x ) = x ( 1 x ) and u ( t , 0 ) = u ( t , 1 ) = 0.

Figures 6(a), 6(c), and 6(e) show time series solutions in the space–time 32 × 32 domain of size ( L , T ) = ( 1 , 1 ) with qubit-layer count ( n , l ) = ( 5 , 5 ) for fractional order α = { 1 , 0.8 , 0.6 } and kinematic viscosity ν = 0.02. The Reynolds number is R e 2 u c L / ν = 100, where the characteristic velocity u c = 1. In the inviscid case, an initial sine-wave profile tends toward an N-wave shock solution. As with Fig. 3, reducing the fractional order leads to fast decay at short times and slow convergence at long times,46 shown here for α { 1 , 0.8 , 0.6 }. Insets show the relative deviations from the classical finite difference solution are less than 2%, without any clear dependence on α.

Fig. 6.

Time-series solutions of 1D Burgers' equation plotted in time intervals of 1 / 32 in (a), (c), and (e) N × M = 32 × 32 and (b), (d), and (f) N × M = 32 × 128 space–time domain of size ( L , T ) = ( 1 , 1 ) for qubit-layer count ( n , l ) = ( 5 , 5 ), kinematic viscosity ν = 0.02 ( R e = 100), and fractional order α (a) and (b) 1, (c) and (d) 0.8, and (e) and (f) 0.6. Insets (a), (c), and (e) show the time-averaged relative deviation from the classical finite difference solution on the same domain, and insets (b), (d), and (f) show the relative norm convergence with M at times t { 1 / 32 , 1 / 16 }. Initial and boundary conditions are u ( 0 , x ) = sin ( 2 π x ) and u ( t , 0 ) = u ( t , 1 ) = 0.

Fig. 6.

Time-series solutions of 1D Burgers' equation plotted in time intervals of 1 / 32 in (a), (c), and (e) N × M = 32 × 32 and (b), (d), and (f) N × M = 32 × 128 space–time domain of size ( L , T ) = ( 1 , 1 ) for qubit-layer count ( n , l ) = ( 5 , 5 ), kinematic viscosity ν = 0.02 ( R e = 100), and fractional order α (a) and (b) 1, (c) and (d) 0.8, and (e) and (f) 0.6. Insets (a), (c), and (e) show the time-averaged relative deviation from the classical finite difference solution on the same domain, and insets (b), (d), and (f) show the relative norm convergence with M at times t { 1 / 32 , 1 / 16 }. Initial and boundary conditions are u ( 0 , x ) = sin ( 2 π x ) and u ( t , 0 ) = u ( t , 1 ) = 0.

Close modal

For semi-implicit time-fractional Burgers' equation, the Courant number C u c τ α / h depends on the fractional index α. To verify the dependence of the numerical stability on t α, the number of time steps M is quadrupled from 32 to 128 in the time series solutions, as shown in Figs. 6(b), 6(d), and 6(f). By inspection, time-fractional solutions are sensitive to time step τ. Specifically, for α = 0.6, convergence in solution accuracy requires M N [see Fig. 6(f) inset].

3. Trace error and optimization

Figure 7(a) shows for ν 0.04, the maximum trace error max ϵ tr for normal Burgers' equation ( α = 1) is significantly greater than for the time-fractional Burgers' equation ( α < 1), the latter which converges to steady state solutions faster (Fig. 6). Figure 7(b) shows the required number of function evaluations N eval scales with both ν and α, due to the difficulty of convergence toward unstable solutions driven by the explicit advection term.

Fig. 7.

α ν contour plots of (a) the maximum trace error max ϵ tr and (b) the number of function evaluations N eval for variational quantum solutions of 1D Burgers' equation (same as Fig. 6).

Fig. 7.

α ν contour plots of (a) the maximum trace error max ϵ tr and (b) the number of function evaluations N eval for variational quantum solutions of 1D Burgers' equation (same as Fig. 6).

Close modal

Together, these results suggest a trade-off between solution fidelity and gradient evaluation costs, depending on fractional index α and diffusion parameter ν.

A natural extension to the nonlinear differential equations involves exploring coupled systems of differential equations, as exemplified by the susceptible–exposed–infectious–recovered (SEIR) epidemic model.44 Fractional-order epidemic models have garnered significant attention in the wake of the COVID-19 pandemic.47–49 These models include spatiotemporal models that account for the geographical spread of the disease.50–52 A minimal representation of a fractional diffusive SEIR model is given by [Ref. 53, Eqs. (24)–(29)]
D t α 1 S = ν 1 x x S β I S μ S + Π ,
(31)
D t α 2 E = ν 2 x x E + β I S ( σ + μ ) E ,
(32)
D t α 3 I = ν 3 x x I + σ E ( ρ + μ ) I ,
(33)
D t α 4 R = ν 4 x x R + ρ I μ R ,
(34)
which partition a population into four distinct cohorts: susceptible individuals S, exposed individuals E, infectious individuals I, and those who have recovered R, without accounting for vaccination and quarantine. The symbols α i and ν i, where i { 1 , 2 , 3 , 4 }, are the fractional order and the diffusion coefficient corresponding to each respective { S , E , I , R } cohort in that order, Π is the population influx rate, β is the infective rate, μ is the death rate, σ is the progression rate, and ρ is the recovery rate.

1. Coupled semi-implicit scheme

The coupled system of Eqs. (31)–(34) is discretized via finite difference and expressed in the following iterative form with common time step τ:
( β I k 1 + μ + g α 1 , τ ν 1 h 2 δ ) S k = g α 1 , τ f α 1 ( S [ 0 , k 1 ] ) + Π ,
(35)
( σ + μ + g α 2 , τ ν 2 h 2 δ ) E k = g α 2 , τ f α 2 ( E [ 0 , k 1 ] ) + β I k 1 S k 1 ,
(36)
( ρ + μ + g α 3 , τ ν 3 h 2 δ ) I k = g α 3 , τ f α 3 ( I [ 0 , k 1 ] ) + σ E k 1 ,
(37)
( μ + g α 4 , τ ν 4 h 2 δ ) R k = g α 4 , τ f α 4 ( R [ 0 , k 1 ] ) + ρ I k 1 ,
(38)
where the operator f α ( u [ 0 , k 1 ] ) = w k [ α ] u 0 j = 1 k 1 Δ w j [ α ] u k j. Since the β I S term is strictly positive, it can be iterated as β I k 1 S k [Eq. (35)]. With only an O ( 1 ) number of cohorts, the coupled system of equations can be solved using a separate cost function for each cohort at time k, i.e.,
C S k = ( r S k ) 2 2 ( S k | A S | S k + β S k | Λ ̃ I k 1 | S k ) r S k ( g α 1 , τ S k , f ̃ S , α 1 k 1 + Π S k , + ̃ n ) ,
(39)
C E k = ( r E k ) 2 2 E k | A E | E k r E k ( g α 2 , τ E k , f ̃ E , α 2 k 1 + β E k , Λ ̃ I k 1 S ̃ k 1 ) ,
(40)
C I k = ( r I k ) 2 2 I k | A I | I k r I k ( g α 3 , τ I k , f ̃ I , α 3 k 1 + σ I k , E ̃ k 1 ) ,
(41)
C R k = ( r R k ) 2 2 R k | A R | R k r R k ( g α 4 , τ R k , f ̃ R , α 4 k 1 + ρ R k , I ̃ k 1 ) ,
(42)
where | + ̃ n 2 n | + n is the (unnormalized) equal superposition over all computational basis states, and each subscript label S { S , E , I , R } corresponds to the indicated cohort. The Hamiltonian matrices N × N are A S = ( μ + g α 1 , τ ) I ν 1 / h 2 L, A E = ( σ + μ + g α 2 , τ ) I ν 2 / h 2 L, A I = ( ρ + μ + g α 3 , τ ) I ν 3 / h 2 L, and A R = ( μ + g α 4 , τ ) I ν 4 / h 2 L, where L is a symmetric tridiagonal matrix with 2 along its main diagonal and 1 along the adjacent off-diagonals, and I = I n is an n-qubit identity matrix. Terms such as S ( θ k ) | Λ I k 1 | S ( θ k ) can be measured using the circuit shown in Fig. 12(d).

2. Basic reproduction number

The spread of a transmissible disease is governed by an important quantity known as the basic reproduction number,54 
R = β σ Π μ ( μ + σ ) ( μ + ρ ) ,
(43)
which is defined, in epidemiological terms, as the average number of secondary cases produced by one infected individual in a susceptible population. Depending on the value of R, the number of infectious individuals may increase resulting in an epidemic ( R > 1), or decrease resulting in the disease dying out ( R < 1).

3. Test case

Consider a set of 2 n spatially connected populations, each with a small infected cohort and a much larger susceptible cohort. Following Ref. 55, we assume the following set of SEIR parameters: Π = 750,53, μ = 0.033 25, β = 5.1 × 10 6, σ = 0.17, and ρ = 0.1109, leading to R 0.66. Individuals are assumed to move around within a region [0, 1] according to Fick's law, with diffusion coefficients ν 1 = ν 2 = 10 3 and ν 3 = ν 4 = 5 × 10 4. Initial cohort sizes are S ( 0 , x ) = 22 500 ( Π / μ for near-steady populations55), I ( 0 , x ) = 20 e 2 x (for spatial dependence53), and E ( 0 , x ) = R ( 0 , x ) = 1. Neumann boundary conditions are applied, i.e., x S ( t , 0 ) = x S ( t , L ) = 0, where S { S , E , I , R }.

In a 16 × 32 domain of size ( L , T ) = ( 1 , 100 ) with ( n , l ) = ( 4 , 5 ), Fig. 8 verifies that R < 1 results in a monotonic decrease (a), (c), and (e) in the infectious cohort I, whereas R > 1 results in an increase (b), (d), and (f) in I before eventual decrease at long times (not shown). A decrease in the fractional index α flattens the spatiotemporal profile of I, in such a way that I decreases more slowly for R < 1, and increases more slowly for R > 1.

Fig. 8.

Infectious I cohort (actual values) in the 16 × 32 space–time domain of size ( L , T ) = ( 1 , 100 ) with qubit-layer count ( n , l ) = ( 4 , 5 ) for each cohort, depending on basic reproduction number (a), (c), and (e) R = 0.66 and (b), (d), and (f) R = 1.33, and fractional index (a) and (b) α = 1, (c) and (d) α = 0.8, and (e) and (f) α = 0.6.

Fig. 8.

Infectious I cohort (actual values) in the 16 × 32 space–time domain of size ( L , T ) = ( 1 , 100 ) with qubit-layer count ( n , l ) = ( 4 , 5 ) for each cohort, depending on basic reproduction number (a), (c), and (e) R = 0.66 and (b), (d), and (f) R = 1.33, and fractional index (a) and (b) α = 1, (c) and (d) α = 0.8, and (e) and (f) α = 0.6.

Close modal

Similar time-fractional memory effects also apply to other cohorts, namely susceptible S, exposed E and recovered R, as shown in Fig. 9 at the start of an epidemic ( R > 1). The apparent smoothness of the solutions, as presented, affirms the potential of the iterative solver in other applications involving coupled fractional-order derivatives, such as predator–prey population models56 and Oldroyd-B viscoelastic flows.2 

Fig. 9.

(a) and (b) Susceptible S, (c) and (d) exposed E, and (e) and (f) recovered R cohorts (actual values) in 16 × 32 space–time domain of size ( L , T ) = ( 1 , 100 ) with qubit-layer count ( n , l ) = ( 4 , 5 ) for each cohort at the start of an epidemic ( R = 1.33 > 1). Fractional index (a), (c), and (e) α = 1 and (b), (d), and (f) α = 0.8.

Fig. 9.

(a) and (b) Susceptible S, (c) and (d) exposed E, and (e) and (f) recovered R cohorts (actual values) in 16 × 32 space–time domain of size ( L , T ) = ( 1 , 100 ) with qubit-layer count ( n , l ) = ( 4 , 5 ) for each cohort at the start of an epidemic ( R = 1.33 > 1). Fractional index (a), (c), and (e) α = 1 and (b), (d), and (f) α = 0.8.

Close modal

Near-term quantum devices are subjected to hardware noise, which significantly impairs the performance of variational quantum algorithms,57,58 despite their inherent resilience to coherent errors.59 To explore the effect of noise, we employ a realistic noise model sampled from the IBM Nairobi quantum hardware device (Falcon r5.11H, version 1.3.3) applied to the Aer Simulator backend provided by the Pennylane-Qiskit plugin. Circuits are transpiled using default settings with level 3 optimization and SABER routing. The optimizer of choice is the simultaneous perturbation stochastic approximation (SPSA), which is known for its exceptional scalability and performance under noisy conditions.57,58 We apply the SPSAOptimizer class on Pennylane using default hyperparameters for 200 iterations, equivalent to 400 circuit evaluations. Data readout is assumed to be noiseless via statevector emulation.

Here, we solve minimally for the subdiffusion problem (Sec. IV A) in N × M = 4 × 2 using qubit-layer count ( n , l ) = ( 2 , 1 ) for α = 1. Initial conditions are u ( 0 , x ) = x ( 1 x ) and boundary conditions are periodic instead of Dirichlet. Figure 10 shows the effects of the hardware noise model on Aer Simulator based on 40 instances with fixed 10 000 shots per circuit evaluation. Of note here is the effect of noise on the evaluation of the norm, which directly affects the weights of the cost function C for the subsequent time step. For fair evaluation of each time step, we restore the norm to the classical error-free value before each time step.

Fig. 10.

Hardware noise model tests across 40 instances on Aer Simulator using the SPSA optimizer set at 200 iterations per time step. (a) Effect of the noise model on overlap and Hamiltonian measurements results in shifting of (b) the optimized cost function C. (c) Box plots of norm fidelity φ norm and trace error ϵ tr with the noise model in time steps k = { 1 , 2 }. Reference values without the noise model are shown as φ norm 1 and ϵ tr 0. (d) displays plots of the norm fidelity φ norm k and the trace error ϵ tr k for time steps k { 1 , 2 }. The trace error ϵ tr k at time step k is defined in Eq. (26), while the norm fidelity is defined as φ norm k=|rk/|uk|| where rk is the measured norm form Eq. (14) and | u k| is the classical norm.

Fig. 10.

Hardware noise model tests across 40 instances on Aer Simulator using the SPSA optimizer set at 200 iterations per time step. (a) Effect of the noise model on overlap and Hamiltonian measurements results in shifting of (b) the optimized cost function C. (c) Box plots of norm fidelity φ norm and trace error ϵ tr with the noise model in time steps k = { 1 , 2 }. Reference values without the noise model are shown as φ norm 1 and ϵ tr 0. (d) displays plots of the norm fidelity φ norm k and the trace error ϵ tr k for time steps k { 1 , 2 }. The trace error ϵ tr k at time step k is defined in Eq. (26), while the norm fidelity is defined as φ norm k=|rk/|uk|| where rk is the measured norm form Eq. (14) and | u k| is the classical norm.

Close modal
Figure 10(a) shows that the averaged overlap measurements are shifted by the noise model to u 1 , u 0 = 0.686 ± 0.049 and u 2 , u 1 = 0.682 ± 0.049, from approximately 1 for measurements without the noise model. The observed relative noise error is approximately 30 %, which is significantly greater than the error from Hamiltonian measurements u k | A | u k with the noise model [less than 10 % by inspection, Fig. 10(b)]. This can be attributed to greater gate errors from the Toffoli (controlled cnot) gates found in controlled entangling layers of the overlap circuit. In addition, the fractional noise error of the cost function in Eq. (15),
( η C ) 2 = ( 2 η O ) 2 + ( η H ) 2 ,
(44)
depends on η O, the fractional noise error of the overlap measurements u k , u k 1 , and η H, the fractional noise error of the Hamiltonian measurements u k | A | u k . Given η O / η H 3, the proportion of noise error of the cost function due to the noise of the overlap measurement is 36 / 37 0.986. This means that the overall performance of the algorithm depends almost entirely on the noise error of the overlap measurements and very little on that of the Hamiltonian measurements.
Figure 10(c) shows the norm fidelity,
φ norm k = | r k u k | ,
(45)
where r k is the measured norm Eq. (14) and u k is the classical norm, and trace error ϵ tr k Eq. (26) at time step k { 1 , 2 }. Here, the mean norm fidelities φ ¯ norm k 0.65 are nearly identical for k { 1 , 2 } due to the classical norm reset, whereas the mean trace errors are ϵ ¯ tr 1 0.039 and ϵ ¯ tr 2 0.052.

Figure 10(d) shows the effects of hardware noise on the quantum solution states | u k over 40 instances. Although the degree of data scatter at each node (width of box plot) is of a similar magnitude to the difference between consecutive classical solutions, the mean readout states | u ¯ k are in general agreement with the respective classical solutions. Further tests were conducted on actual hardware, in this case, the IBMQ Mumbai 27 qubit device across 20 instances (see inset) using Qiskit Runtime sampler and estimator primitives for overlap and Hamiltonian measurements, respectively, set at 200 SPSA iterations for a single time step.

In this study, we proposed and implemented a hybrid variational quantum algorithm for solving time-fractional diffusive differential equations in engineering applications, including nonlinear equations such as Burgers' equation (Sec. IV B) and coupled systems of equations such as those in epidemic models (Sec. IV C). Our approach is not only efficient in time complexity but also offers advantages in terms of memory utilization. While the classical memory cost of numerically solving these equations scales linearly with the number of spatial grid points N, the quantum memory cost scales only logarithmically with N. This helps to mitigate the global dependence problem in fractional calculus, where previous solutions must be repeatedly accessed from memory, causing the classical solver to have prohibitive memory costs when the spatial range is large. In terms of time complexity, the present iterative scheme is limited by the quadratic scaling with the number of time steps; this has some implications for the numerical stability of nonlinear fractional equations (see Fig. 6).

The choice of the ansatz in the form of a parameterized circuit is reflected in its expressibility and in the entangling capability performance60 of a variational quantum algorithm.9 Hence, as a matter of heuristics, the real-amplitude ansätze used in this study [Eqs. (19) and (20)] are by no means the most efficient designs for such applications. Designing and optimizing ansatz architecture, including utilizing adaptive methods61 or leveraging geometric tools,62–64 remains an active area of research.65–67 

While the hybrid approach offers notable advantages, it is not exempt from the challenges that are commonly encountered in variational quantum algorithms. First, unlike a classical computer, the cost of loading initial conditions as classical data as an n-qubit quantum state could be prohibitive.68 Encoding schemes, such as parametric optimization [Eq. (21)], may incur an exponential cost in primitive operations since the encoded n-qubit state occupies a space of dimension O ( 2 n ). More efficient encoding techniques are currently being developed; for example, solutions may be represented in the form of finite Fourier series, using the quantum Fourier transform as part of the Fourier series loader introduced by Moosa et al.69 and employed by Sarma et al.14 The use of the Walsh series instead of the Fourier series may yield an even more efficient subexponential scaling in encoding.70 In addition, extracting the complete solution at each of the M time steps would require full readout of the statevector via quantum state tomography. This experimental procedure is notorious for its exponential cost,71 which is a bottleneck to quantum advantage as pointed out in Ref. 68. Therefore, where possible, applications should consider a partial readout of the (spatial) solutions only at a few selected time points of interest to minimize the expensive overhead due to tomography, or measuring only integrated quantities like mean, variance, or other moments. In these cases, robust classical shadow tomography techniques72–76 can be employed to further reduce the required sample complexity.

In closing, the hybrid algorithm investigated in this work presents a nascent but promising pathway to solving fractional partial differential equations of the integro-differential form, subjected to potentially prohibitive space or memory requirements. More work needs to be done to address issues related to the trainability of variational quantum algorithms, such as barren plateaus and narrow gorges,39,77 as well as practical hardware implementation, including error mitigation78 and correction.

We acknowledge the use of IBM Quantum services for this work. The views expressed are those of the authors, and do not reflect the official policy or position of IBM or the IBM Quantum team. This research is supported by the National Research Foundation, Singapore and the Agency for Science, Technology and Research (A*STAR) under the Quantum Engineering Programme (NRF2021-QEP2-02-P03), and A*STAR C230917003. This project is supported by the Singapore Ministry of Health's National Medical Research Council through the Programme for Research in Epidemic Preparedness and Response (PREPARE), under Environmental Transmission & Mitigation Co-operative (PREPARE-CS1-2022-004). DEK acknowledges funding support from the A*STAR Central Research Fund (CRF) Award for Use-Inspired Basic Research. DP acknowledges support from the Ministry of Education Singapore, under Grant No. MOE-T2EP50120-0019.

The authors have no conflicts to disclose.

Fong Yew Leong: Conceptualization (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Validation (lead); Writing – original draft (lead); Writing – review & editing (equal). Dax Enshan Koh: Conceptualization (supporting); Formal analysis (supporting); Investigation (supporting); Methodology (supporting); Writing – original draft (supporting); Writing – review & editing (equal). Jian Feng Kong: Formal analysis (supporting); Methodology (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Siong Thye Goh: Investigation (supporting); Validation (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Jun Yong Khoo: Investigation (supporting); Supervision (equal); Writing – original draft (supporting); Writing – review & editing (supporting). Wei-Bin Ewe: Investigation (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Hongying Li: Methodology (supporting); Resources (equal); Writing – review & editing (supporting). Jayne Thompson: Project administration (equal); Supervision (equal); Writing – review & editing (supporting). Dario Poletti: Project administration (equal); Resources (equal); Supervision (equal); Writing – review & editing (supporting).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

In this appendix, we give explicit circuit diagrams for some of the circuits we used in this paper. Real-amplitude hardware-efficient ansätze are shown in Fig. 11. Circuits used for measuring expectations and overlaps are shown in Fig. 12.

Fig. 11.

(a) Real-amplitude hardware-efficient ansatz [Eq. (19)] and (b) circular real-amplitude ansatz [Eq. (20)].

Fig. 11.

(a) Real-amplitude hardware-efficient ansatz [Eq. (19)] and (b) circular real-amplitude ansatz [Eq. (20)].

Close modal
Fig. 12.

Circuits for measuring (a) overlap u ( θ k ) , u j , (b) expectation U ( θ k ) | A | U ( θ k ) [Eq. (15)], (c) nonlinear overlap u k , Λ ± k 1 u k 1 [Eq. (30)], (d) expectation S ( θ k ) | Λ I k 1 | S ( θ k ) [Eq. (39)], and (e) higher order overlap u k , B u k 1 [Eq. (B4)]. U ( θ k ) denotes an ansatz with parameters θ optimized at time k and U j denotes an ansatz with parameters fixed at time j. I / S denotes an optional shift operator [Eq. (18)] and S / S denotes either an incremental or a decremental shift operator. Operators S and I are the ansätze representing the respective variables.

Fig. 12.

Circuits for measuring (a) overlap u ( θ k ) , u j , (b) expectation U ( θ k ) | A | U ( θ k ) [Eq. (15)], (c) nonlinear overlap u k , Λ ± k 1 u k 1 [Eq. (30)], (d) expectation S ( θ k ) | Λ I k 1 | S ( θ k ) [Eq. (39)], and (e) higher order overlap u k , B u k 1 [Eq. (B4)]. U ( θ k ) denotes an ansatz with parameters θ optimized at time k and U j denotes an ansatz with parameters fixed at time j. I / S denotes an optional shift operator [Eq. (18)] and S / S denotes either an incremental or a decremental shift operator. Operators S and I are the ansätze representing the respective variables.

Close modal
The quadratic scaling of quantum complexity in time steps in Eq. (25) motivates higher order time-stepping methods in order to limit time discretization and its errors. The Crank–Nicolson method is a second-order finite-difference scheme that retains the stability advantage associated with implicit methods.79 Using the trapezoidal rule, the finite difference of the spatial derivative in Eq. (5) takes the midpoint between forward and backward Euler schemes, yielding
x x u ( t k , x n ) 1 2 h 2 [ ( u n + 1 k 2 u n k + u n 1 k ) + ( u n + 1 k 1 2 u n k 1 + u n 1 k 1 ) ] .
(B1)
Equating to the finite difference of the Caputo derivative, we obtain
( 1 a 2 δ ) u n k = a 2 δ u n k 1 + w k u n 0 j = 1 k 1 Δ w j u n k j ,
(B2)
or in matrix form,
A u k = B u k 1 + w k u 0 j = 1 k 1 Δ w j u k j ,
(B3)
where A = I B N × N, B = ( a / 2 ) L N × N, L is a symmetric tridiagonal matrix with 2 along its main diagonal and 1 along the adjacent off-diagonals, and I is the identity matrix. This can be iterated using variational quantum optimization in Eq. (9) as before, with an updated cost function
C k = 1 2 ( r k 1 u k , B u k 1 + u k , f ̃ k 1 ) 2 u k | A | u k ,
(B4)
which includes a new term u k , B u k 1 to be measured as a sum of observables [2–4, Eq. (17)], using the quantum circuit [Fig. 12(e)].
Since the coefficient Δ w j of the ( k j ) th term decreases monotonically with j for constant fractional power index α, the O ( M 2 ) time complexity scaling Eq. (25) can be limited to O ( M ) by setting an upper bound to the limits of the Caputo integral. This is done by replacing the Caputo finite difference derivative in Eq. (4) with
D t α u ( t k , x n ) g α , τ [ w 1 u n k + j = 1 min { k , ξ } ( w j + 1 w j ) u n k j ] ,
(C1)
where ξ [ 1 , M ] and w k + 1 = 0. For ξ < k, the truncation error thus incurred is
ε tr ( α , ξ ) O [ g α , τ ( Δ w ξ + 1 ) r k ( ξ + 1 ) ] ,
(C2)
where Δ w ξ + 1 = w ξ + 2 w ξ + 1. A similar trade-off between solution fidelity and time complexity applies to classical finite difference schemes.
1.
T.
Liu
and
M.
Hou
,
Adv. Math. Phys.
2017
,
1
.
2.
F.
Wang
and
Y.
Wang
,
Adv. Math. Phys.
2023
,
1
.
3.
O. S.
Iyiola
and
F. D.
Zaman
,
AIP Adv.
4
,
107121
(
2014
).
4.
V. E.
Lynch
,
B. A.
Carreras
,
D.
Del-Castillo-Negrete
,
K. M.
Ferreira-Mejias
, and
H. R.
Hicks
,
J. Comput. Phys.
192
,
406
(
2003
).
5.
C.
Li
and
F.
Zeng
,
Numer. Funct. Anal. Optim.
34
,
149
(
2013
).
6.
D. A.
Murio
,
Comput. Math. Appl.
56
,
1138
(
2008
).
7.
Y.
Lin
and
C.
Xu
,
J. Comput. Phys.
225
,
1533
(
2007
).
8.
K.
Bharti
,
A.
Cervera-Lierta
,
T. H.
Kyaw
,
T.
Haug
,
S.
Alperin-Lea
,
A.
Anand
,
M.
Degroote
,
H.
Heimonen
,
J. S.
Kottmann
et al,
Rev. Mod. Phys.
94
,
015004
(
2022
).
9.
M.
Cerezo
,
A.
Arrasmith
,
R.
Babbush
,
S. C.
Benjamin
,
S.
Endo
,
K.
Fujii
,
J. R.
McClean
,
K.
Mitarai
,
X.
Yuan
et al,
Nat. Rev. Phys.
3
,
625
(
2021
).
10.
H.-M.
Li
,
Z.-X.
Wang
, and
S.-M.
Fei
,
Phys. Rev. A
108
,
032418
(
2023
).
11.
Y.
Sato
,
R.
Kondo
,
S.
Koide
,
H.
Takamatsu
, and
N.
Imoto
,
Phys. Rev. A
104
,
052409
(
2021
).
12.
M.
Ali
and
M.
Kabel
,
Phys. Rev. Appl.
20
,
014054
(
2023
).
13.
M.
Lubasch
,
J.
Joo
,
P.
Moinier
,
M.
Kiffner
, and
D.
Jaksch
,
Phys. Rev. A
101
,
010301
(
2020
).
14.
A.
Sarma
,
T. W.
Watts
,
M.
Moosa
,
Y.
Liu
, and
P. L.
McMahon
, “
Quantum variational solving of nonlinear and multi-dimensional partial differential equations
,” arXiv:2311.01531 (
2023
).
15.
W.-B.
Ewe
,
D. E.
Koh
,
S. T.
Goh
,
H.-S.
Chu
, and
C. E.
Png
,
IEEE Trans. Microwave Theory Tech.
70
,
2517
(
2022
).
16.
F. Y.
Leong
,
D. E.
Koh
,
W.-B.
Ewe
, and
J. F.
Kong
,
Int. J. Numer. Methods Heat Fluid Flow
33
,
3669
(
2023
).
17.
D.
Jaksch
,
P.
Givi
,
A. J.
Daley
, and
T.
Rung
,
AIAA J.
61
,
1885
(
2023
).
18.
P.
Rigas
, “
Variational quantum algorithm for measurement extraction from the Navier-Stokes, Einstein, Maxwell, Boussniesq-type, Lin-Tsien, Camassa-Holm, Drinfeld-Sokolov-Wilson, and Hunter-Saxton equations
,” arXiv:2209.07714 (
2022
).
19.
F. Y.
Leong
,
W.-B.
Ewe
, and
D. E.
Koh
,
Sci. Rep.
12
,
10817
(
2022
).
20.
K.
Kubo
,
K.
Miyamoto
,
K.
Mitarai
, and
K.
Fujii
,
IEEE Trans. Quantum Eng.
4
,
3100717
(
2023
).
21.
F.
Mainardi
,
G.
Pagnini
, and
R.
Gorenflo
,
Appl. Math. Comput.
187
,
295
(
2007
).
22.
S.
Taha Abdulazeez
and
M.
Modanli
,
Alexandria Eng. J.
61
,
12443
(
2022
).
23.
Alternatively, | f ̃ k 1 can be prepared separately through postselection of linear combination of unitaries80 using a quantum multiplexer, to be measured directly as one single overlap term u k | f ̃ k 1 .
24.
D.
Maslov
,
Phys. Rev. A
93
,
022311
(
2016
).
25.
H.
Alghassi
,
A.
Deshmukh
,
N.
Ibrahim
,
N.
Robles
,
S.
Woerner
, and
C.
Zoufal
,
Quantum
6
,
730
(
2022
).
26.
For example, in decreasing index order, we have Π i = 3 1 A i = A 3 A 2 A 1.
27.
G.
Brassard
,
P.
Hoyer
,
M.
Mosca
, and
A.
Tapp
,
Contemp. Math.
305
,
53
(
2002
).
28.
Y.
Suzuki
,
S.
Uno
,
R.
Raymond
,
T.
Tanaka
,
T.
Onodera
, and
N.
Yamamoto
,
Quantum Inf. Process.
19
,
75
(
2020
).
29.
S.
Aaronson
and
P.
Rall
, “
Quantum approximate counting, simplified
,” in
Symposium on Simplicity in Algorithms
(
SIAM
,
2020
), pp.
24
32
.
30.
T.
Giurgica-Tiron
,
I.
Kerenidis
,
F.
Labib
,
A.
Prakash
, and
W.
Zeng
,
Quantum
6
,
745
(
2022
).
31.
D. E.
Koh
,
G.
Wang
,
P. D.
Johnson
, and
Y.
Cao
,
J. Math. Phys.
63
,
052202
(
2022
).
32.
G.
Wang
,
D. E.
Koh
,
P. D.
Johnson
, and
Y.
Cao
,
PRX Quantum
2
,
010346
(
2021
).
33.
P. D.
Johnson
,
A. A.
Kunitsa
,
J. F.
Gonthier
,
M. D.
Radin
,
C.
Buda
,
E. J.
Doskocil
,
C. M.
Abuan
, and
J.
Romero
, “
Reducing the cost of energy estimation in the variational quantum eigensolver algorithm with robust amplitude estimation
,” arXiv:2203.07275 (
2022
).
34.
J. C.
Spall
,
IEEE Trans. Automat. Control
37
,
332
(
1992
).
35.
Y.
Sato
,
R.
Kondo
,
I.
Hamamura
,
T.
Onodera
, and
N.
Yamamoto
, “
Hamiltonian simulation for time-evolving partial differential equation by scalable quantum circuits
,” arXiv:2402.18398 (
2024
).
36.
J.
Hu
,
S.
Jin
,
N.
Liu
, and
L.
Zhang
, “
Quantum Circuits for partial differential equations via Schrödingerisation
,” arXiv:2403.10032 (
2024
).
37.
V.
Bergholm
,
J.
Izaac
,
M.
Schuld
,
C.
Gogolin
,
M. S.
Alam
,
S.
Ahmed
,
J. M.
Arrazola
,
C.
Blank
,
A.
Delgado
et al, “
Pennylane: Automatic differentiation of hybrid quantum-classical computations
,” arXiv:1811.04968 (
2018
).
38.
Z.
Hu
and
S.
Kais
,
Commun. Phys.
6
,
68
(
2023
).
39.
M.
Cerezo
,
A.
Sone
,
T.
Volkoff
,
L.
Cincio
, and
P. J.
Coles
,
Nat. Commun.
12
,
1791
(
2021
).
40.
Y.
Xu
and
O. P.
Agrawal
,
Fractional Calculus Appl. Anal.
16
,
709
(
2013
).
41.
N.
Sugimoto
,
J. Fluid Mech.
225
,
631
(
1991
).
42.
R.
Garra
,
Phys. Rev. E
84
,
036605
(
2011
).
43.
V.
Mukundan
and
A.
Awasthi
,
Nonlinear Eng.
5
,
219
(
2016
).
44.
J. P.
Liu
,
H. Ø.
Kolden
,
H. K.
Krovi
,
N. F.
Loureiro
,
K.
Trivisa
, and
A. M.
Childs
,
Proc. Natl. Acad. Sci. U. S. A.
118
,
e2026805118
(
2021
).
45.
F.
Oz
,
R. K. S. S.
Vuppala
,
K.
Kara
, and
F.
Gaitan
,
Quantum Inf. Process.
21
,
30
(
2022
).
46.
T.
Akram
,
M.
Abbas
,
M. B.
Riaz
,
A. I.
Ismail
, and
N. M.
Ali
,
Alexandria Eng. J.
59
,
2201
(
2020
).
47.
Z.
Ali
,
F.
Rabiei
,
M. M.
Rashidi
, and
T.
Khodadadi
,
Eur. Phys. J. Plus
137
,
395
(
2022
).
48.
S.
Rezapour
,
H.
Mohammadi
, and
M. E.
Samei
,
Adv. Differ. Equation
2020
,
490
.
49.
M. A.
Khan
and
A.
Atangana
, “
Modeling the dynamics of novel coronavirus (2019-nCov) with fractional derivative
,”
Alexandria Eng. J.
59
,
2379
(
2020
).
50.
B.
Barnes
,
M.
Anokye
,
M. M.
Iddrisu
,
B.
Gawu
, and
E.
Afrifa
,
Adv. Math. Phys.
2023
,
1
.
51.
N.
Nazia
,
Z. A.
Butt
,
M. L.
Bedard
,
W.-C.
Tang
,
H.
Sehar
, and
J.
Law
,
Int. J. Environ. Res. Public Health
19
,
8267
(
2022
).
52.
B.
Hu
,
P.
Ning
,
J.
Qiu
,
V.
Tao
,
A. T.
Devlin
,
H.
Chen
,
J.
Wang
, and
H.
Lin
,
Int. J. Infect. Dis.
110
,
247
(
2021
).
53.
A. I.
Alaje
and
M. O.
Olayiwola
,
Healthcare Anal.
4
,
100230
(
2023
).
54.
P.
van den Driessche
,
Infect. Dis. Modell.
2
,
288
(
2017
).
55.
B.
Kammegne
,
K.
Oshinubi
,
O.
Babasola
,
O. J.
Peter
,
O.
Babatope Longe
,
R.
Bosede Ogunrinde
,
E. O.
Titiloye
,
R. T.
Abah
, and
J.
Demongeot
,
Pathogens
12
,
88
(
2023
).
56.
H.
Jafari
,
R. M.
Ganji
,
N. S.
Nkomo
, and
Y. P.
Lv
,
Results Phys.
27
,
104456
(
2021
).
57.
A.
Pellow-Jarman
,
I.
Sinayskiy
,
A.
Pillay
, and
F.
Petruccione
,
Quantum Inf. Process.
20
,
202
(
2021
).
58.
X.
Bonet-Monroig
,
H.
Wang
,
D.
Vermetten
,
B.
Senjean
,
C.
Moussa
,
T.
Bäck
,
V.
Dunjko
, and
T. E.
O'Brien
,
Phys. Rev. A
107
,
032407
(
2023
).
59.
E.
Fontana
,
N.
Fitzpatrick
,
D. M.
Ramo
,
R.
Duncan
, and
I.
Rungger
,
Phys. Rev. A
104
,
022403
(
2021
).
60.
S.
Sim
,
P. D.
Johnson
, and
A.
Aspuru-Guzik
,
Adv. Quantum Technol.
2
,
1900070
(
2019
).
61.
H. R.
Grimsley
,
S. E.
Economou
,
E.
Barnes
, and
N. J.
Mayhall
,
Nat. Commun.
10
,
3007
(
2019
).
62.
A.
Katabarwa
,
S.
Sim
,
D. E.
Koh
, and
P.-L.
Dallaire-Demers
,
Quantum
6
,
782
(
2022
).
63.
T.
Haug
,
K.
Bharti
, and
M. S.
Kim
,
PRX Quantum
2
,
040309
(
2021
).
64.
T.
Haug
and
M. S.
Kim
,
Phys. Rev. A
106
,
052611
(
2022
).
65.
J.
Qin
,
J. Phys.: Conf. Ser.
2634
,
012043
(
2023
).
66.
J.-B.
You
,
D. E.
Koh
,
J. F.
Kong
,
W.-J.
Ding
,
C. E.
Png
, and
L.
Wu
, “
Exploring variational quantum eigensolver ansatzes for the long-range XY model
,” arXiv:2109.00288 (
2021
).
67.
J.
Tilly
,
H.
Chen
,
S.
Cao
,
D.
Picozzi
,
K.
Setia
,
Y.
Li
,
E.
Grant
,
L.
Wossnig
,
I.
Rungger
et al,
Phys. Rep.
986
,
1
(
2022
).
68.
S.
Aaronson
,
Nat. Phys.
11
,
291
(
2015
).
69.
M.
Moosa
,
T. W.
Watts
,
Y.
Chen
,
A.
Sarma
, and
P. L.
McMahon
,
Quantum Sci. Technol.
9
,
015002
(
2024
).
70.
J.
Zylberman
and
F.
Debbasch
,
Phys. Rev. A
109
,
042401
(
2024
).
71.
R.
O'Donnell
and
J.
Wright
, “
Efficient quantum tomography
,” in
Proceedings of the forty-eighth annual ACM symposium on Theory of Computing
,
2016
.
72.
H.-Y.
Huang
,
R.
Kueng
, and
J.
Preskill
,
Nat. Phys.
16
,
1050
(
2020
).
73.
K.
Bu
,
D. E.
Koh
,
R. J.
Garcia
, and
A.
Jaffe
,
npj Quantum Inf.
10
,
6
(
2024
).
74.
S.
Chen
,
W.
Yu
,
P.
Zeng
, and
S. T.
Flammia
,
PRX Quantum
2
,
030348
(
2021
).
75.
D. E.
Koh
and
S.
Grewal
,
Quantum
6
,
776
(
2022
).
76.
D.
Grier
,
H.
Pashayan
, and
L.
Schaeffer
, “
Sample-optimal classical shadows for pure states
,”
Quantum
8
,
1373
(
2024
).
77.
A.
Arrasmith
,
Z.
Holmes
,
M.
Cerezo
, and
P. J.
Coles
,
Quantum Sci. Technol.
7
,
045015
(
2022
).
78.
Z.
Cai
,
R.
Babbush
,
S. C.
Benjamin
,
S.
Endo
,
W. J.
Huggins
,
Y.
Li
,
J. R.
McClean
, and
T. E.
O'Brien
,
Rev. Mod. Phys.
95
,
045005
(
2023
).
79.
U.
Ali
,
F. A.
Abdullah
, and
A. I.
Ismail
,
J. Interpolation Approximation Sci. Comput.
2017
,
18
.
80.
A. M.
Childs
and
N.
Wiebe
,
Quantum Inf. Comput.
12
,
901
(
2012
).