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.

## I. INTRODUCTION

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 \alpha / d x \alpha $, where $\alpha $ 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 equations^{1} 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 $\alpha $, 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* problem^{7} 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 linear^{10–12} and nonlinear differential equations^{13,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.

## II. FRACTIONAL DIFFUSION EQUATIONS

^{6,21}in space

*x*and time

*t*,

*t*, defined as

### A. Difference scheme

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 \xd7 M$. The spatial grid points are defined as $ x i = i h$ where $ h = L / N$ and $ i = 1 , 2 , \u2026 , N$, and the temporal grid points as $ t k = k \tau $ where $ \tau = T / M$ and $ k = 1 , 2 , \u2026 , M$.

^{5,22}

*A*is an $ N \xd7 N$ symmetric tridiagonal matrix with $ b = 1 + 2 a$ on the main diagonal, $ \u2212 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

### B. Variational quantum optimization

*k*by an ansatz or trial state formed by a set of parameterized unitaries

*R*and unparameterized entangling unitaries

*W*as

*l*is the number of layers in the ansatz. The unnormalized source state $ | f \u0303 k \u2212 1 \u27e9$ is a linear combination of history states in $ [ 0 , k \u2212 1 ]$.

*k*,

*k*measurable overlap terms $ \u27e8 u k | u [ 0 , k \u2212 1 ] \u27e9$.

^{23}

*N*spatial grid points.

## III. IMPLEMENTATION

### A. Hamiltonian decomposition

*A*[Eq. (8)] can be decomposed into

^{11}evaluating the expectation value of

*A*requires only the evaluation of the expectation values of the simple terms ( $ H 1 \u2212 2$ for the periodic boundary condition, $ H 1 \u2212 3$ for the Dirichlet boundary condition, and $ H 1 \u2212 4$ for the Neumann boundary condition). The operator $S$ denotes the

*n*-qubit cyclic shift operator,

^{24}

### B. Ansatz initialization

*l*repeating blocks, each consisting of a parameterized

*R*layer with one $ R Y ( \theta )$ gate on each qubit, followed by an unparameterized

*W*layer with a cnot gate between consecutive qubits,

^{25}

^{,}

^{26}In addition, we consider a real-amplitude ansatz with circular entanglement, such that

^{16}See Fig. 11 for circuit diagrams of the ansätze.

### C. Gradient estimation

*k*as

*i*th $ R Y$ gate as

*i*is divided by

*n*and $ \u230a i / n \u230b$ is the integer part of $ i / n$. Thus, for each index

*i*at time

*k*, gradient estimation requires at least

*k*overlap measurements.

### D. Time complexity and memory scaling

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

*M*, the overall time complexity is

*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 2

*l*parametric gates in depth and for the Hamiltonian $ ( l + n 2 )$ based on $ O ( 1 )$ circuits [ $ 2 \u2212 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 ( \u03f5 \u2212 2 )$ shots per measurement, where $\u03f5$ is the desired precision. We remark that with quantum amplitude estimation,

^{27}the number of shots required may be quadratically reduced from $ O ( \u03f5 \u2212 2 )$ to $ O ( \u03f5 \u2212 1 )$. This decrease, however, entails increased circuit depth, as exemplified by the $ \Omega ( \u03f5 \u2212 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.

^{11}). Using the simultaneous perturbation stochastic approximation (SPSA),

^{34}only two circuit executions are required per gradient evaluation $ N g \u223c 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 3^{36}).

*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 \u223c l$, the memory cost of the Caputo sums scales subexponentially as $ O ( k \u2009 log \u2009 N )$ (Fig. 1), demonstrating efficient algorithmic space complexity.

## IV. NUMERICAL EXPERIMENTS

The circuits are implemented in the quantum software framework *Pennylane*^{37} 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 \u2212 6$ and a gradient tolerance of $ 10 \u2212 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 , \u2026 , 0 )$.

### A. Subdiffusion equation

We solve for the time-fractional $ ( 0 < \alpha \u2264 1 )$ subdiffusion equation Eq. (1) as an initial value problem with an initial condition $ u ( 0 , x ) = x ( 1 \u2212 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 \xd7 32$ space–time domain with qubit-layer count $ ( n , l ) = ( 5 , 4 )$, for (a) $ \alpha = 1$ and (b) $ \alpha = 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.

*k*,

*M*for $ 0 < \alpha \u2264 1$, where $ | u \u0302 k \u27e9 \u2261 u k / \Vert u k \Vert $ 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 ( \theta 1 )}$ to $ q 1 : { R Y ( \theta 2 ) \u223c R Y ( \theta n )}$, which effectively increases the connectivity of $ q 1$ by $ n \u2212 2$ qubits. This may improve the performance of an ansatz for certain problems.^{38}

Figure 4(b) shows that the number of function evaluations $ \u2211 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 $ \u03f5 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.

### B. Time-fractional Burgers' equation

^{40,41}has found applications in shallow water waves and waves in bubbly liquids.

^{42}We consider the following 1D time-fractional Burgers' equation

^{43}as an extension to the time-fractional diffusion equation Eq. (1),

#### 1. Semi-implicit scheme

*u*along its diagonal, and the subscript $\xb1$ denotes either an incremental or a decremental shift in each eigenstate under periodic boundary conditions, i.e., $ | 2 n \u27e9 = | 0 \u27e9$.

^{13,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 \pi x )$ and Dirichlet boundary conditions $ u ( 0 , x ) = x ( 1 \u2212 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 \xd7 32$ domain of size $ ( L , T ) = ( 1 , 1 )$ with qubit-layer count $ ( n , l ) = ( 5 , 5 )$ for fractional order $ \alpha = { 1 , 0.8 , 0.6}$ and kinematic viscosity $ \nu = 0.02$. The Reynolds number is $ R e \u2261 2 u c L / \nu = 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 $ \alpha \u2208 { 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 $\alpha $.

For semi-implicit time-fractional Burgers' equation, the Courant number $ C \u223c u c \tau \alpha / h$ depends on the fractional index $\alpha $. To verify the dependence of the numerical stability on $ t \alpha $, 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 $\tau $. Specifically, for $ \alpha = 0.6$, convergence in solution accuracy requires $ M \u226b N$ [see Fig. 6(f) inset].

#### 3. Trace error and optimization

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

Together, these results suggest a trade-off between solution fidelity and gradient evaluation costs, depending on fractional index $\alpha $ and diffusion parameter $\nu $.

### C. Fractional diffusive 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)]

*S*, exposed individuals

*E*, infectious individuals

*I*, and those who have recovered

*R*, without accounting for vaccination and quarantine. The symbols $ \alpha i$ and $ \nu i$, where $ i \u2208 { 1 , 2 , 3 , 4}$, are the fractional order and the diffusion coefficient corresponding to each respective $ { S , E , I , R}$ cohort in that order, $\Pi $ is the population influx rate, $\beta $ is the infective rate, $\mu $ is the death rate, $\sigma $ is the progression rate, and $\rho $ is the recovery rate.

#### 1. Coupled semi-implicit scheme

*k*, i.e.,

*n*-qubit identity matrix. Terms such as $ \u27e8 S ( \theta k ) | \Lambda I k \u2212 1 | S ( \theta k ) \u27e9$ can be measured using the circuit shown in Fig. 12(d).

#### 2. Basic reproduction number

^{54}

#### 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: $ \Pi = 750$,^{53}^{,} $ \mu = 0.033 \u2009 25$, $ \beta = 5.1 \xd7 10 \u2212 6$, $ \sigma = 0.17$, and $ \rho = 0.1109$, leading to $ R \u2248 0.66$. Individuals are assumed to move around within a region [0, 1] according to Fick's law, with diffusion coefficients $ \nu 1 = \nu 2 = 10 \u2212 3$ and $ \nu 3 = \nu 4 = 5 \xd7 10 \u2212 4$. Initial cohort sizes are $ S ( 0 , x ) = 22 \u2009 500$ ( $ \u223c \Pi / \mu $ for near-steady populations^{55}), $ I ( 0 , x ) = 20 e \u2212 2 x$ (for spatial dependence^{53}), and $ E ( 0 , x ) = R ( 0 , x ) = 1$. Neumann boundary conditions are applied, i.e., $ \u2202 x S ( t , 0 ) = \u2202 x S ( t , L ) = 0$, where $ S \u2208 { S , E , I , R}$.

In a $ 16 \xd7 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 $\alpha $ 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$.

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 models^{56} and Oldroyd-B viscoelastic flows.^{2}

## V. HARDWARE NOISE

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 \xd7 M = 4 \xd7 2$ using qubit-layer count $ ( n , l ) = ( 2 , 1 )$ for $ \alpha = 1$. Initial conditions are $ u ( 0 , x ) = x ( 1 \u2212 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.

Figure 10(d) shows the effects of hardware noise on the quantum solution states $ | u k \u27e9$ 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 \xaf k \u27e9$ 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.

## VI. DISCUSSION

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 performance^{60} 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 methods^{61} 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 techniques^{72–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 mitigation^{78} and correction.

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

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

## DATA AVAILABILITY

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

### APPENDIX A: QUANTUM CIRCUIT DIAGRAMS

### APPENDIX B: HIGHER ORDER CRANK–NICOLSON SCHEME

^{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

### APPENDIX C: TRUNCATED CAPUTO INTEGRAL

*j*for constant fractional power index $\alpha $, 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

## REFERENCES

^{80}using a quantum multiplexer, to be measured directly as one single overlap term $\n\u27e8\n\nu\nk\n|\n\n\n\nf\n\u0303\n\nk\n\u2212\n1\n\u27e9$.