# Review and perspectives in quantum computing for partial differential equations in structural mechanics

^{1}Aerospace Structures and Computational Mechanics, Aerospace Structures and Materials, Faculty of Aerospace Engineering, Delft University of Technology, Delft, Netherlands^{2}Numerical Analysis, Applied Mathematics, Faculty of Electrical Engineering Mathematics and Computer Science, Delft University of Technology, Delft, Netherlands^{3}Aerodynamics, Flow Physics and Technologies, Faculty of Aerospace Engineering, Delft University of Technology, Delft, Netherlands

Structural mechanics is commonly modeled by (systems of) partial differential equations (PDEs). Except for very simple cases where analytical solutions exist, the use of numerical methods is required to find approximate solutions. However, for many problems of practical interest, the computational cost of classical numerical solvers running on classical, that is, silicon-based computer hardware, becomes prohibitive. Quantum computing, though still in its infancy, holds the promise of enabling a new generation of algorithms that can execute the most cost-demanding parts of PDE solvers up to exponentially faster than classical methods, at least theoretically. Also, increasing research and availability of quantum computing hardware spurs the hope of scientists and engineers to start using quantum computers for solving PDE problems much faster than classically possible. This work reviews the contributions that deal with the application of quantum algorithms to solve PDEs in structural mechanics. The aim is not only to discuss the theoretical possibility and extent of advantage for a given PDE, boundary conditions and input/output to the solver, but also to examine the hardware requirements of the methods proposed in literature.

## 1 Introduction

Structural mechanics is often modeled by means of partial differential equations (PDEs). However, it is only for a small set of simple problems, domains and boundary condition that an analytical solution is known. In all other cases, engineers must rely on numerical techniques to obtain an approximate solution.

Given the importance of PDEs, research on convergence and accuracy of numerical solvers has dominated in the past decades. Nevertheless, the computational demand is still high whenever the domain is large with respect to the physics’ scale or when nonlinearities require a high level of detail of the numerical solution.

Potentially, quantum computing can revolutionize the field of numerical computational mechanics, thanks to its promise of unprecedented theoretical speed-up. For instance, the Harrow, Hassidim, Lloyd (HHL) algorithm (Harrow et al., 2009) can solve sparse linear systems exponentially faster than any classical method. At a first glance, it may seem that HHL could be applied to the discretized Poisson equation and exponentially reduce the runtime to obtain the displacement field of a loaded structure.

As tantalizing as they sound, almost all promises of quantum advantage are currently bound to theory and simulations and likely will be for the next few decades. What has yet to catch up with the algorithms is the hardware, which is still far from the technological requirements for practical quantum advantage.

This limitation gave reason to a separate branch of research, which abandoned the idea of proving quantum speed-up, but focused on algorithms that take into account the limitations of current hardware. These are mostly hybrid methods, that use a classical computer in combination with a quantum one, assigned to solve classically hard tasks. As with the theoretical speed-up algorithms, these hybrid methods also found application in the field of numerical PDE solving.

Given the importance and richness of the field, this review wants to be the first work to present and compare different quantum PDE solvers. Specifically, this work is aimed at the applied mechanics community and thus limits its focus to equations in structural mechanics. However, it must be pointed out that the field of quantum algorithms for PDEs spans many other branches of science, such as fluid mechanics (Mezzacapo et al., 2015; Steijl and Barakos, 2018; Todorova and Steijl, 2020; Gaitan, 2020, 2021; Budinski, 2021a,b), finance (An et al., 2020; Chakrabarti et al., 2021; Fontanela et al., 2021), model discovery (Heim et al., 2021) and cosmology (Mocz and Szasz, 2021), which may also benefit from specialized reviews.

This paper is structured as follows. Section 2 presents some standard concepts of quantum computing and explains the quantum algorithms at the root of PDE-solving. Section 3 reviews the literature on quantum algorithms for linear PDEs, while Section 4 is devoted to nonlinear PDEs. Finally, Section 5 provides concluding remarks and discusses possible future prospects of quantum computation applied to structural problems.

## 2 Main concepts of quantum PDE solving

At the time of writing, almost no quantum algorithm applies to all possible combinations of partial differential equations, boundary conditions, discretizations, etc. However, one can trace out a generic PDE-solving workflow, as in Figure 1. Here it is important to notice that quantum computation takes place only in the so-called *quantum primitive*, i.e. an algorithm that acts on quantum states by means of quantum operators. Also, the quantum primitive does not alone solve the PDE, but rather its discrete form, which generally is the computational bottleneck in classical logic.

**FIGURE 1**. Workflow of the process of solving partial differential equations with a quantum primitive. The input is given, in the general case, by the PDE together with its boundary (BCs) and initial (ICs) conditions. After the differential problem has been discretized with schemes such as finite elements or finite differences, the input data is encoded into a quantum state (state preparation step). The quantum primitive then produces the solution to the discrete problem as a quantum state. However, the information in this state cannot be accessed, as this is generally in quantum superposition. Therefore, quantum PDE solver also needs to measure the state after the quantum primitives as many times as required by the user-defined solution precision. Depending of the application, such measurement are postprocessed to provide the final output.

Furthermore, in contrast to a classical algorithm, a quantum primitive does not operate in the same environment in which the data is stored and read-out. In other words, there is a state-preparation step prior to quantum computation, where classical data is encoded as a quantum state, as well as a measurement step afterwards, that provides output in classical form.

A similar workflow was previously proposed (Pesah, 2020), where the author points out that a discretized PDE can be mapped either to a Schrödinger equation or to one or more linear system(s). In the first case, one may use Hamiltonian simulation to obtain the solution in quantum state, while linear systems could be solved by the HHL algorithm or more efficient quantum linear solvers. However, this classification concerns only linear PDEs and does not consider algorithms that are not fully quantum (i.e. all computations from quantum-form data to quantum-form solution are done on a quantum computer) or quantum annealing.

This section gives the necessary background on the ‘quantum-block’ of PDE solving, that concerns the last three steps in Figure 1.

### 2.1 Quantum state preparation

The problem of quantum state preparation arises every time classical data needs to be encoded as a quantum state and it is by no means restricted to quantum PDE solvers. Taking the homogeneous heat equation as an example, one needs to encode an initial condition *u*_{0}(*x*) in a quantum register. First, space-discretization transforms *u*_{0}(*x*) into the discrete array of gridpoint values **u**_{0} of *N* entries. Then, this vector is normalized to turn it into the suitable quantum state

The vector on the left hand side of Eq. 1 is called *ket vector*, according to the so-called *braket* notation (Nielsen and Chuang, 2010).

Furthermore, *N* complex amplitudes defined with respect to a chosen basis, i.e.

where _{2} will be subsituted by log in the following, for simplicity of notation). By labelling the states of a qubit as 0 and 1, the basis states are written as

which constitute the so-called *computational basis*. It is easy to see that the computational basis forms an orthonormal basis, since, given two basis vectors *b*_{i} and *b*_{j}

where *δ*_{ij} is the Dirac delta.

One should keep in mind that the computational basis is just one of the infinite possible bases in the

Furthermore, the braket notation implies a unit norm vector, that is

Equivalently, Eq. 2 states that the squared amplitudes of a ket vector can be seen as probability amplitudes, modelling the fact that a *n*-qubits system will collapse to the basis state

Having access to the amplitues, the problem of state preparation becomes to evolve an initial state up to the state

where *U* is a unitary operator representing the quantum state evolution and

Clearly, the cost of implementing *U* influences the overall cost of solving the differential problem. In fact, preparing a generic quantum state over *n* qubits requires *O*(*N*) quantum gates (Nielsen and Chuang, 2010), where the “big *O*” notation represents the asymptotic upper bound on the number of computational resources. Therefore, even if the quantum primitive runs in

### 2.2 Hamiltonian simulation

One of the motivations that drove to the study of quantum computers was to simulate quantum systems, that are hard to simulate classically. Generally speaking, one aims at describing the evolution of a quantum state

where *H* is the system’s Hamiltionian, i.e. a Hermitian matrix. For the sake of explanation, the Hamiltonian is considered here as time-invariant.

Hamiltonian simulation consists in evolving an initial quantum state

Eq. 3 has the exact solution

where *e*^{−iHt} is a unitary operator, due to the fact that *H* is Hermitian. Therefore, *e*^{−iHt} is a valid quantum operation that could be applied to *H* is difficult to simulate directly.

However, one can see *H* as a sum of Hermitian matrices that are efficient to simulate. In that case

and *H*_{i} do not generally commute, i.e.

for some *i* and *j*. Anyway, by dividing the evolution time in *r* subsequent time intervals Δ*t* = *t*/*r*, the single *H*_{i} Hamiltonians can be evolved separately over Δ*t*. This evolution over the sub-interval is then repeated *r* times. In the end, the overall operator is an approximation of the actual evolution of the total Hamiltonian *H* over *t* as stated, for instance, by the Trotter’s formula (Trotter, 1959):

The importance of *H* and *H*_{i} being easy to simulate is clear from the Hamiltonian simulation runtime using product formulas. This is *O*(*f*(*n*)*t*/*ɛ*) (Nielsen and Chuang, 2010), where *n* is the number of qubits and *ɛ* is the desired Hamiltonian simulation error. While *f*(*n*) = poly(*n*) for simulatable Hamiltonians, *f*(*n*) = exp(*n*) in the general case, making for an exponential difference in the runtime.

In quantum physics, many natural systems are described by so-called *local Hamiltonians*, which are operators that act nontrivially on a few qubits and are known to be efficiently simulatable (Nielsen and Chuang, 2010). However and most importantly for PDEs, it was shown that also sparse Hamiltonians can be simulated in

Finally, we remark that recent work exponentially reduced the dependency on precision (Berry et al., 2014; 2015a) and also achieved an optimal dependency on the sparsity (Berry et al., 2015a). Most recently, an approach based on qubitization achieved an additive lower bound with respect to *t* and 1/*ɛ* (Low and Chuang, 2019).

### 2.3 Quantum linear solvers

Linear systems arise in several numerical methods for PDEs and their solution often dominates the cost of the overall classical PDE algorithm. After appropriate discretization, the PDE assumes the familiar form *A***x** = **b**, which is a linear system in the *N*-dimensional space, *N* being the number of unknowns, which can correspond to the number of gridpoints, particles or harmonic basis functions, depending on the method of discretization used.

If the matrix *A* is positive definite, the fastest general-purpose classical linear system solver is the conjugate gradient method (Shewchuk, 1994). This algorithm has asymptotic time complexity

where *s* is the matrix sparsity, *κ* is the conditioning number and *ɛ* is the desired error for a certain precision metric.

What motivated the study of linear solvers beyond classical computation is the conjecture, due to complexity arguments, that classical algorithms cannot invert *A* in time *N* (Harrow et al., 2009). In particular the time complexity of the HHL algorithm is

However, the HHL algorithm and subsequent quantum linear solvers, do not solve *A***x** = **b**, but rather

known as *quantum linear system problem* (QLSP). It is important to keep in mind that quantum linear systems algorithms (QLSAs) output the solution

The complexity of solving the QLSP was improved by later contributions. In terms of conditioning number, Ambainis et al. reduced the dependency from *O*(*κ*^{2}) to *O*(*κ*) using variable time amplitude amplification (VTAA) for post-selecting the solution (Ambainis, 2012). Also concerning the conditioning number, Clader et al. showed how to include the use of a preconditioner in quantum linear solvers in order to deal with ill-conditioned matrices (Clader et al., 2013). Later on, an approach based on linear combination of unitaries was proposed to replace the phase estimation step, reducing the dependency on precision to

Finally, several authors recently proposed to solve the QLSP with variational quantum algorithms (VQAs) (An and Lin, 2019; Bravo-Prieto et al., 2019; Huang et al., 2021; Xu et al., 2021; Patil et al., 2022), to be discussed later. These algorithms are appealing since they are suitable for implementation on near-term devices and require only shallow quantum circuits. However, VQAs are heuristics, meaning that the number of iterations needed for convergence depends on the choice of the optimization strategy. Thus, no theoretical runtime is known *a-priori* for VQAs. However, a recent study showed that VQA-based linear solvers include a quantum state verification step that requires at least ^{1} number of queries to the unitary preparing the state

### 2.4 Amplitude amplification and amplitude estimation

Amplitude amplification is an algorithm that iteratively promotes the probability of one or more amplitudes of a quantum state. This technique derives directly from Grover’s algorithm for database search (Grover, 1996, 1997). In this problem, a database and a function *f* are given. The database contains *N* = 2^{n} entries, each of which is a *n*-bits string, such as

The function *f* marks the target bitstring *m* as follows

Classical search methods require to check the bitstring individually and therefore scale as *O*(*N*) and require *N* − 1 queries in the worst case. On the other hand, Grover’s algorithm works with all database entries in superposition and progressively promotes the amplitude of the target state

where each *j* represents one of the *N* bitstrings.

The quantum search algorithm consists in applying the *Grover’s operator*

*k* times, where the outer product notation is used for ^{2}. It is easy to show that

Even though *G* contains the target state *G* uses an *oracle* for *f*, which is an operator such that

where the tensor product is implied between neighbouring ket vectors.

Therefore, the knowledge of *O*_{f}, which flips an ancilla qubit when the state register is in the target state. The specifics of the oracles depend on the different problems and go beyond the scope of this discussion, but suffices to say that they generally entail a unitary contolled by a *n* = ⌈ log *N*⌉ register, such that the unitary gets applied only when such register is in state

The concept of quantum search was extended to multiple target states with quantum amplitude amplification (QAA) (Brassard, et al. 2000). In this case, the Hilbert space spanned by the database is divided into a “good” subspace, containing the target states and a complementary “bad” subspace. Therefore, the Grover’s operator becomes

where *P*_{g} is the projector onto the “good” subspace. Also in this case a target state can be prepared quadratically faster than by using a classical search.

QAA is used as a subroutine in many quantum algorithms that require post-selection. This is usually a procedure used as the last step of a quantum computation, when a quantum register of interest contains the solution state with finite probability and it is entangled to an ancilla qubit in such a way that measuring the latter in one of its two states will collapse the register to the solution state. Therefore post-selecting requires to repeat the entire quantum routine until the correct measurement of the ancilla is obtained. For instance, in the case of many quantum linear solvers, the state before post-selection is

where

Quantum amplitude estimation (QAE) allows to compute the probability of one component of a quantum state. The main idea is to perform quantum phase estimation on Grover’s operator. In fact, the action of *G* on a state *θ* in the plane defined by the components of

where *Y* is the Pauli-*Y* operator, i.e.

**FIGURE 2**. Effect of one Grover iteration used in quantum amplitude amplification. The axes represent the components of

Therefore, quantum phase estimation on *G* returns an approximation of 2*θ*. In turn, this allows to compute the probability of the ‘good’ component of

where

Quantum amplitude estimation is useful, for instance, to estimate the integral of a PDE solution over a certain region *S*, i.e.

### 2.5 Variational quantum algorithms

Variational quantum algorithms (VQAs) are the main class of methods conceived to run on current or near-future devices, better known as Noisy Intermediate Scale Quantum (NISQ) (Preskill, 2018). These early-stage quantum computers use up to a few hundreds of the so-called *physical qubits*, which are affected by noise of different nature, as opposed to *logical qubits*, that instead are made of many physical qubits that act as a fully error-corrected qubit.

Besides being suitable for near-term applications, VQAs are thought to be candidates for quantum advantage, in areas such as quantum chemistry, nuclear physics, but also optimization and machine learning (Cerezo et al., 2021a).

The small circuit depth is achieved via a hybrid strategy internal to an optimization process. At every iteration, a cost function is evaluated by the quantum computer and fed to the classical one, which updates the design parameters according to an optimization algorithm. Figure 3 schematically shows the working principle of VQAs.

**FIGURE 3**. Working principle of variational quantum algorithms. The input are the hyperparameters, such as ansatz, optimizer and type of cost function and the initial value of the ansatz parameter. The cost is shown here as a linear combination of expectation values of unitaries, although other expressions are possible. In the quantum-classical optimization loop, the classical computer is in charge of updating the parameters ** θ** according to the optimization logic, while the quantum computer evaluates the cost function terms. At the end of the optimization, the optimal parameter set

*θ*_{opt}can be used to reconstruct the approximate solution as a wavefunction

*g*(

*u*).

A comprehensive review of VQAs, their main concepts, applications and future prospects is given in (Cerezo et al., 2021a). In what follows, only a succint explanation is given about the main elements of these algorithms.

#### 2.5.1 Ansatz

The ansatz is a parametrized unitary *U*(** θ**) that encodes the tentative solution to the problem, depending on the parameters

**. At the end of the optimization, the optimal parameters will define the solution’s approximation as**

*θ**U*(

*θ*_{opt}).

There exist many different families of ansatze, but a first important distinction is between the so-called *problem-inspired* ones, that incorporate information about the problem, and the *problem-agnostic* ones. For instance, a good example of the problem agnostic class is the so-called *hardware-efficient ansatz* (Kandala et al., 2017), which optimizes the use and distribution of gates according to hardware specification. On the other hand, examples of problem-inspired ansatze are the Unitary Coupled Clustered ansatz (Taube and Bartlett, 2006) and the Quantum Alternating Operator Ansatz (Farhi et al., 2014; Hadfield et al., 2019).

#### 2.5.2 Cost function

The cost function *C*(** θ**) is a metric of how far the tentative solution at the current iteration is from the actual solution of the problem. In VQAs, the cost function is computed by measuring one or more observables of the ansatz state

Some of the requirements for a good cost function apply to quantum as well as classical methods. For instance, the minimum of *C*(** θ**) must coincide with the solution of the problem and decreasing values of the cost function should correspond to better approximations of the solution. However, a good cost function for VQAs also needs to be hard to evaluate classically, since it would otherwise negate the possibility of quantum advantage (Cerezo et al., 2021a).

#### 2.5.3 Gradient

Many optimization routines use the cost function gradient to identify the direction of maximum descent to speed-up convergence. Often however, the cost function does not have an analytical expression and the gradient needs to be computed using finite differences. Nevertheless, a quantum-evaluated cost function is noisy in its nature, due to finite sampling and hardware noise, so the finite difference approximation can be highly imprecise.

Luckily, the so-called *parameter-shit rule* (PSR) allows to compute gradients analytically using two cost function evaluations, similar to what happens in finite differences. The derivative of the cost *C* with respect to the *l*th parameter is then

where *θ*^{±} = *θ* ± *α***e**_{l}, **e**_{l} is the vector with 1 in the *l*th entry and 0 in all others and *α* is the magnitude of the shift.

Even though Eq. 13 looks similar to the finite difference formula, the two differ by the term *α* = *π*/2, the difference in statistical error between finite differences and PSR is maximum (Mari et al., 2021).

#### 2.5.4 Optimizer

Optimizers for VQAs need to account for different complications that arise from a quantum-evaluated cost landscape. In fact, this is usually affected by stochastic and hardware noise or it can be flat with minima in narrow gorges, due to the phenomenon of *barren plateaus* (McClean et al., 2018; Cerezo et al., 2021b).

A main classification of VQA optimizers is between gradient-based and gradient-free algorithms. In principle, any classical optimizer can be used in VQAs, but the noisy character of the cost function makes some choices better suited than others to avoid stability and convergence issues. For instance, gradient-based techniques generally implement some form of stochastic gradient descent (SGD) methods, such as Adam (Kingma and Ba, 2015) or use natural gradients (Stokes et al., 2020). On the other hand, the simultaneous perturbation stochastic approximation (SPSA) method is probably the most popular gradient-free technique (Spall, 1992).

### 2.6 Quantum annealing

Quantum annealing (QA) is an algorithm that can solve combinatorial optimization (CO) problems. These require to search among a finite set of variables or choices with the aim is to find the “best” one according to a certain metric of merit. Importantly, many combinatorial optimization problems are NP-hard, meaning that finding their solution is difficult. More precisely, NP-hard problems are at least as difficult to solve as the most difficult problem in the NP complexity class, which includes all problems that can be verified in polynomial time by a nondeterministic algorithm. The NP class should be seen in contrast to the P class, which instead includes all problems that can be solved “efficiently”, i.e. in polynomial time using a deterministic algorithm. Given their hardness, CO problems are practially solved approximately, and the challenge is to produce accurate approximate solutions in short (polynomial) time.

However, QA cannot encode and solve a combinatorial optimization problem in generic form. Instead, the physical problem solved by QA is the one of finding the ground state of the Ising Hamiltonian (McGeoch, 2020).

where *q*_{i}, *q*_{j} represent the binary values associated to the qubits, e.g. *i* and *i* and *j*.

In order to reach to the ground state of *H*_{Ising}, the quantum annealing process starts from an initial Hamiltonian *H*_{0}, whose ground state is known and easy to prepare. A common choice for *H*_{0} is the transverse magnetic field

*X*_{i} being the Pauli-*X* operator acting on the *i*th qubit, where

From here, the system’s Hamiltonian is slowly changed according to an evolution law, such that

Here, *T* is the total evolution time, *f*(0) = 0, *f*(*T*) = 1. If *H*(*t*), known as the *total Hamiltonian*, is changed slowly enough, then the adiabatic principle ensures that the system’s configuration is at the ground state of *H*(*t*) at all times (Kato, 1950). At the end of adiabatic evolution, the system will therefore be in the ground state of the Ising Hamiltonian, which corresponds to the approximate solution of the CO problem.

However, quantum annealers operate with spin systems, where each spin *s*_{i} is a binary variable equal to ±1. Thus, QA programmers need to recast the CO problem to a binary optimization one. In D-Wave machines, which are the state-of-the-art quantum annealers, the problem is given as input in the Quadratic Unconstrained Binary Optimization (QUBO) form, that is

where *s*_{i} = 2*q*_{i} − 1.

### 2.7 Measurement

At the end of the quantum primitive the problem’s solution is encoded as a quantum state. This is a unit vector *O*(*N*) measurements, breaking any speed-up obtainable with a quantum primitive.

Therefore, quantum computing becomes appealing for PDEs when the aim is to compute a certain scalar function *g*(*u*) of the solution. This function is often taken as the expectation value of an observable *M*, i.e.

In principle, estimating ⟨*M*⟩ would require to measure *M*. However, this basis is not known in general and one needs to resort to methods other than direct measurement. One popular technique is the *Hadamard Test* (Aharonov et al., 2009), which assumes that *M* is a unitary operator. As showed in Figure 4, *M* is controlled on an ancilla qubit set in uniform superposition and later measured in the Pauli-*X* basis. It is easy to show that

where Pr(0) and Pr(1) are the probabilities of the ancilla qubit collapsing to

**FIGURE 4**. Hadamard Test circuit for computing ⟨*u*|*M*|*u*⟩. Oppositely to standard observable measurements, *M* needs to be unitary.

A recurrent case in applications is to estimate how “close” a state *SWAP Test* (Buhrman et al., 2001), which is a Hadamard test where *M* is equal to the *SWAP* gate. If *SWAP* gate acts between them as

In the SWAP Test, the *SWAP* operator is actually generalized between the quantum registers containing

In practice, all expectation values are determined statistically as averages. The number of samples *N*_{m} required to approximate ⟨*M*⟩ up to *ɛ* precision is given by the Chernoff-Hoeffding inequality, for which

where *O*(1/*ɛ*^{2}) measurements are necessary to esimtate ⟨*M*⟩ up to precision *ɛ*.

## 3 Quantum algorithms for linear PDEs

### 3.1 Poisson equation

The Poisson equation describes different phenomena in solid mechanics, such as the displacements of a solid undergoing loads, the elongation of a truss, etc. Let *d* is the number of spatial dimensions and

Furthermore, Eq. 21 must be complemented with Dirichlet, Neumann or Robin conditions on the boundary *δ*Ω in order to find a specific solution.

Assuming an hypercubic domain, the standard technique to solve Eq. 21 numerically is to discretize the domain in *N* grid points in each of the *d* spatial dimensions and either approximate the Laplacian with a central finite difference (FD) scheme or employ a finite element (FE) approximation. In both cases, one obtains a discrete equation, or Discrete Poisson Equation, that is a linear systems of

where *u*_{i} = *u*(*x*^{(i)}) is the solution at grid point *x*^{(i)} and *f*_{i} = *f*(*x*^{(i)}) for finite differences or **f** = ∫_{Ω}*f*(*x*)*φ*_{i}(*x*) for finite elements with basis functions *φ*_{i}(*x*). Eq. 22 can be solved classically using direct or iterative solvers. However, all classical linear solvers have a cost upper bound by a polynomial in *N* (*O*(poly(*N*)) complexity). Even the fastest iterative solvers such as conjugate gradient (Shewchuk, 1994) or multigrid techniques (Golub and van Loan, 2013) require *O*(*N*) iterations.

#### 3.1.1 HHL and preconditioning

It comes natural to think that QLSAs such as HHL (cfr Section 2) could solve the discrete Poisson problem exponentially faster than in classical logic. Indeed, several authors worked on this concept. Clader et al. (2013) were among the first ones to apply an improved version of HHL to a finite difference discretization of the stationary Maxwell equations. These form a system of PDEs that are not of the Poisson type, but still elliptical and that can be reduced to a linear system such as Eq. 22 upon discretization.

As mentioned, Clader et al. did not use the standard HHL, but a modified version of it. Most importantly, they noticed that the *O*(*κ*^{2}) in the cost of HHL nullifies the exponential speed-up in *N*, if *κ* = *O*(*N*). This is indeed the case for matrices induced by a finite element discretization of second order elliptical boundary value problems, for which *κ* = *O*(*N*^{2/d}) (Brenner and Scott, 2010). Therefore, Clader et al. proposed a quantum algorithm to reproduce the Sparse Preconditioner Approximate Inverse (SPAI) (Benzi and Tûma, 1999) operator. This technique uses a matrix *P* to achieve nearly optimal preconditioning, i.e. *PA* ≈ *I*, where *I* is the identity matrix. At the same time, *P* is such that *PA* preserves the sparsity of the *A*. If *s* is the sparsity of *A* in Eq. 22 and the preconditioner *P* has similar sparsity, then applying the *P* such that

can be done in *O*(*s*^{2}) queries to an oracle accessing *PA* and *O*(*s*^{3}) runtime (Clader et al., 2013). Then, if the SPAI procedure is applied succesfully, the conditioning will be independent of *N* or *s* and 1/*ɛ* are also constant or logarithmic in *N*, the complexity of Clader’s method would be *O*(poly log(*N*)), making for a true exponential speed-up. Unfortunately, in most if not all discretization schemes, 1/*ɛ* = *O*(*N*), as will be discussed later in this section.

Further than conditioning, Clader et al. discussed two other caveats of HHL (Aaronson, 2015), namely the state preparation and the measurement problems. Concerning the first, general state preparation is hard and could itself kill any polynomial speed-up. Therefore, Clader et al. proposed to use an oracle able to return

However, even though one would query this oracle a constant number of times to produce

About measurements, Clader et al. (2013) provide examples of classical quantities that can be extracted once *j*th component of

As mentioned, Clader et al. (2013) also describe how the preconditioned HHL can be applied to Maxwell’s equations and, if a scalar solution *N*. However, Clader’s method is very much problem-agnostic and the Maxwell discrete problem is seen more as a linear system for benchmarking their solver, rather than a discretized PDE. Moreover, all speed-up results are based on the existence and the efficiency of oracles, which are quantum black-boxes, instead of actual circuits. In the particular case of Maxwell’s equations, the authors only mention that an oracle implementing *A* and **f** would be ‘efficient’, but do not elaborate further on number of submodules, gates, etc (Clader et al., 2013). Even so however, the complexity analysis in Clader et al. (2013) is incomplete, as later noticed by Montanaro and Pallister (Montanaro and Pallister, 2016). In fact the preconditioned HHL still retains a *α* ∈ (0, 1), this solver still does not achieve an exponential speed-up for elliptical problems.

#### 3.1.2 Diagonalization with quantum fourier transform

While still using a QLSA, Cao et al. focused on the problem rather than the solver (Cao et al., 2013). In their work, they treated the *d*-dimensional Poisson equation (Eq. 22) and noticed that every classical linear solver requires at least

time, therefore suffering from the curse of dimensionality. Using the HHL algorithm, Cao et al. resolved the exponential dependency on *d*, by preparing the solution

quantum operations. It should be noticed how both Eqss 25, 26 are completely expressed in terms of 1/*ɛ* and *d*, thus not hiding any dependencies between complexity terms. Furthermore, the authors described the quantum circuits used down to submodules and gates.

Comparing Eq. 26 with the complexity of the HHL algorithm (Eq. 8), one notices that the *O*(1/*ɛ*) term in the cost of Cao’s algorithm is replaced by a logarithmic dependency. The exponential reduction does not derive from modifications to the linear solver, but from the specific structure of the Poisson matrix. This can be understood starting from the 1-dimensional case, defining *S* (Benedetto, 1997) as

where *ɛ* requires *d*-dimensional Poisson matrix

and its Hamiltonian simulation is

Therefore, simulating

Also, Eq. 26 should not mislead the reader in thinking that the cost of Cao’s method is independent of the conditioning number. Indeed Eq. 26 refers to a single run of the quantum circuit. However, *O*(*κ*^{2}) runs need to be performed, as required by the HHL algorithm to encode the solution state *κ* = *O*(*N*^{d/2}), the exponential advantage is lost, but Clader’s preconditioning technique (Clader et al., 2013) may be used to achieve constant or logarithmic dependency of the conditioning with respect to *N*.

Wang et al. built upon Cao’s work and simplified and improved the HHL-based quantum Poisson solver (Wang et al., 2020c). In particular, they presented a fully modular circuit of this algorithm, defining every module in its components, down to known quantum operations (such as addition and subtraction) and gates. In doing this, they noticed that several steps in Cao’s algorithm could be reduced to the evaluation of trascendental function. For this sake, they developed the so-called *quantum function-value binary expansion* (qFBE) algorithm (Wang et al., 2020a), which allows to replace more expensive power operations with arithmetic ones.

Moreover, the authors of Wang et al. (2020c) reduced the complexity of the controlled rotation operation in Cao’s method. In the latter, after the quantum phase estimation step is performed, the eigenvalue register *C* − *R*_{y}(*θ*_{j}) step is used to encode amplitude 1/*λ*_{j}, if the rotation angle is taken as *θ*_{j} = arcsin(1/*λ*_{j}). Cao et al. used a quantum implementation of the bisection method, to iteratively evaluate the arc sine function, requiring *O*(log^{4}(1/*ɛ*)) operations for a small-dimension, high-precision (*d* < log(1/*ɛ*)) problem. However, Wang et al. bypassed the computation of the reciprocal state completely and noticed that the angles for the controlled rotation could be estimated with minimum error through the following relation

Also, the arc-cotangent function could be estimated through the qFBE algorithm in *O*(log^{3}(1/*ɛ*)) operations.

Overall, Wang’s Poisson solver produces the solution state *O*(*κd* log^{3}(1/*ɛ*)) operations, where the usual consideration about the conditioning number applies. Comparing this scaling with one of Cao’s algorithm (cfr Eq. 26), Wang’s method is polynomially more efficient for low dimensional and high precision problems, which are usually the most interesting ones in structural analysis.

#### 3.1.3 Adaptive order scheme and spectral method

The methods of Cao (Cao et al., 2013) and Wang (Wang et al., 2020c) consider the case when the Laplacian in Eq. 21 is approximated with three points centered FD for all grid sizes. However, Childs et al. recently remarked that fixed finite difference, finite element and finite volume schemes require *O*(poly(1/*ɛ*)) time to bring the approximate solution *ɛ*-close to the actual solution on the grid (Childs et al., 2021). In fact, fixed schemes produce matrices with *κ* = *O*(poly(1/*ɛ*)) and all quantum linear solvers have polynomial time in the conditioning number.

Childs et al. accounted for this issue, by solving Poisson and general second order elliptical boundary value problems with two different numerical approximations, namely an adaptive order FD approximation and a spectral method (Childs et al., 2021). The adaptive FD approach is used to solve the *d*-dimensional Poisson problem. With periodic boundary conditions.

First, Childs shows that the conditioning number *κ* of the order-*k* Laplacian with periodic boundary conditions is *O*(*N*^{2}) if *N*, *k*, *d* and 1/*ɛ*. Choosing *N* is then automatically determined to ensure the total error is upper bound as *O*(*ɛ*):

By substituting Eq. 31 into the runtime of the complexity-optimal QLSA solver of Childs et al. (2017), the solution of the second order elliptical problem with periodic BCs is found in

runtime, meant as number of gate operations, which is polynomial in *d* (i.e. no curse of dimensionality) and has optimal dependency in 1/*ɛ*.

Furthermore, the authors show that the same runtime holds for Dirichlet or Neumann homogeneous boundary conditions. To achieve this, they use the method of images to extend the domain and symmetrize the solution *u* according to the BCs (Childs et al., 2021).

The second algorithm proposed by Childs et al. uses the spectral method (Childs et al., 2021). In this case, the solution is approximated globally and the discretized Laplacian is non-sparse. To obviate to this problem, two variations of the quantum Fourier transform, namely the *quantum shifted Fourier transform* (QSFT) and the *quantum* cosine *transform* (QCT) are used to induce sparsity. The algorithm makes use of an oracle that is queried

times for second-order elliptic problems with non-homogeneous Dirichlet boundary conditions and

times for the Poisson problem with homogeneous Dirichlet BCs.

#### 3.1.4 Full complexity analysis

All previous techniques are based on quantum linear system algorithms and demonstrate times that are exponentially better than classical solvers. Yet it is important to understand what output is prepared in such time and what input provides the starting point. To begin with, all algorithms for solving linear systems *A***x** = **b** requires the right hand side vector is given as an input in normalized form, i.e. **x** and therefore the full discrete solution.

In quantum computation, matters are not as trivial. The main point is that information is processed by quantum computers as quantum states, but it is only accessible and readable in classical form. For instance, the HHL algorithm requires the input vector **b** to be provided as a quantum state **b** might have been previously *stored* as a quantum state *O*(log(*N*)) time (Giovannetti et al., 2008b,a), it is uncertain whether they can be physically built and if they actually offer an advantage if a parallel computer with the same amount of resources is available (Ciliberto et al., 2018). On the other hand, the state *prepared* through a sequence of unitary operations. However, preparation of a generic quantum state from its amplitudes is a hard problem (Nielsen and Chuang, 2010) that would kill any potential exponential saving.

Also, a quantum algorithm does not compute the solution **x**, but rather prepares the state *N*, this will require *O*(*N*) measurements of

The full complexity of a QLE algorithm for differential problems, from input encoding to output readout was discussed by Montanaro and Pallister (Montanaro and Pallister, 2016). Specifically, they considered elliptical PDEs discretized with the finite element method, using piece-wise polynomial basis functions. Since the full quantum solution

Specifically Montanaro and Pallister assumed the final output to be the inner product *r*(*x*) defined over the domain Ω

where *φ*_{i} are the finite element basis functions and the vectors

The main findings of Montanaro and Pallister (2016) are summarized in Table 1. Here, the classical linear solver used is the conjugate gradient method, while the quantum one is that of Childs et al. (2017), which has logarithmic dependency on 1/*ɛ*. Furthermore, the basis functions used for the results in Table 1 are the linear ‘hat’ functions

defined over a uniform grid *x*_{i} ∈ [*x*_{0}, *x*_{1}, *…*, *x*_{N−1}] with equidistant spacing *h*. Montanaro and Pallister (2016) present scalings also for polynomial shape function of generic order *p*.

**TABLE 1**. Complexity results in Montanaro and Pallister (2016). The problem is a second order elliptical PDE, discretized with the finite element method, using piece-wise linear shape functions.

The complexity analysis was done both without preconditioning and with optimal preconditioning (i.e. *PA* = *I*), where a realistic preconditoning case lies in between these two extremes. The most important result in Table 1 is that, for fixed dimension *d*, no exponential quantum speed-up can be achieved, regardless of conditioning. The reason is that to compute *ɛ* requires *d* was allowed to vary, it may look as if the quantum solver is exponentially faster than classical (no curse of dimensionality). However, the authors warn that the *ɛ* but can vary with *d*, making it hard to say what kind of speed-up is achievable for variable dimensions. Furthermore, a classical random walk procedure can solve the Poisson equation in polynomial time, if the solution is required only at a particular point (Linden et al., 2020).

The runtimes in Table 1 scale with 1/*ɛ*, *l*-seminorm and *l*-norm of the analytical solution *u* and are indicative of the magnitudes up to the *l*th derivative of *u*. Therefore, for *d* relatively high, yet fixed, there can be a consistent polynomial speed-up for high precision problems and high second derivatives of the solution *u*(*x*). These considerations extend to higher-order finite elements, with the difference that higher-order Sobolev seminorms appear in the runtime. However, from a structural mechanics standpoint, elliptical problems reach at most *d* = 3, limiting the achievable polynomial speed-up.

It is important to carefully understand why, even considering a log(1/*ɛ*)-scaling algorithm, Montanaro and Pallister found out that quantum linear solvers cannot provide exponential speed-up. To do this, one should answer the questions identified by Aaronson in his ‘fine print’ for quantum linear algebra (Aaronson, 2015).

1. *Can* *be prepared in time O*(log(1/*ɛ*))*, starting from* **b***?* In general, this is a hard problem, but if *f*(*x*) in Eq. 21 is a polynomial or a function supported only on a few elements, the state *O*(log(1/*ɛ*)) (Montanaro and Pallister, 2016).

2. *Is there an algorithm for accessing the elements of A in time O*(log(1/*ɛ*))*?* Yes. In fact, quantum linear solvers require a sparse matrix and a sparse access to the matrix. The second point means that there should be an algorithm that, given a row index *r* and another index *i*, it returns the column index and value of the *i*th non-zero element in *A*. Finite element matrices satisfy both these requirements, since they are sparse and, if the mesh is regular, sparse access can be obtained by the knowledge of the element’s degrees of freedom and by the connectivity matrix.

3. *Is possible to apply efficient pre-conditioning to A?* Yes, for instance through the quantum-SPAI technique proposed by Clader et al. (Clader et al., 2013).

4. *Is it possible to measure the output in time O*(log(1/*ɛ*))*?* Not in general. Especially, it is not the case for estimating properties such as *ɛ*-close to each other require

Thus, the mere fact of having to sample a solution almost always kills any quantum exponential speed-up. Still, Montanaro and Pallister specify that there are cases in which some properties of the solution may be tested with logarithmic sampling. One such property is periodicity, which suggests the possibility of an exponential speed-up when using quantum linear solvers for the sake of finding the vibration frequencies of a structure.

Also, an even more interesting point in Montanaro’s analysis is the fact that a possible exponential speed-up will not arise from a QLE algorithm, but from a clever choice of state preparation or measurement. In fact, replacing a QLE routine with a classical linear solver implies at worst a polynomial slowdown (Montanaro and Pallister, 2016).

#### 3.1.5 NISQ solvers

So far, none of the work discussed concerns the implementation of algorithms on hardware. In the most abstract cases, parts of the algorithm are performed by oracles, which perform a certain function in given time complexity, but whose circuit structure is unknown. Other times the circuit is sketched at different levels of detail. Most notably for the Poisson case, Wang et al. showed the circuit of their algorithm and all its modules up to qubit operations (gates) (Wang et al., 2020c).

However, it is easy to realize that none of these algorithms can possibly run on quantum hardware either today or in the next few decades (Preskill, 2018). For instance, Wang demonstrates his algorithm to solve a minimal problem of a 4 × 4 Poisson matrix (Wang et al., 2020c) on the Sunway TaihuLight supercomputer, which acts as a quantum simulator (Chen et al., 2018). Interestingly, the authors also provide qubit and gate counts for this implementation, declaring 38 qubits and 800 gates. One should consider that the actual gate count is higher, since gates such as *TOFFOLI* and *SWAP* used in Wang et al. (2020c) must be expanded to native operations and all operations have to be mapped to actual hardware connectivity. Therefore, circuits of this size require quantum volumes that are beyond near-term capabilities by orders of magnitude. Consequently, other authors recently looked at the possibility to solve the Poisson equation with near-term techniques.

Wang et al. also proposed a way to solve the one-dimensional Poisson problem in NISQ, using circuits of *O*(poly log(1/*ɛ*)) operations (Wang et al., 2020b). Their main idea was that the gate-expensive Hamiltonian simulation can be bypassed if one is able to directly encode the inverse eigenvalues in the quantum state amplitudes. This is easily understood by recalling the expression of the solution of the

where *λ*_{j} and *j*th eigenvalue and eigenvector of *A* respectively and

Importantly, the Poisson problem with homogeneous Dirichlet boundary conditions results in a matrix whose eigenvalues have an analytical expression, i.e.

where *j* ranges from 1 *N* to 2.

Wang et al. noticed that this Poisson matrix is also a Cartan matrix. The eigenvalues of Cartan matrices are also sines and they are related by the following product relation (Damianou, 2014)

where *n* = log(*N*).

It is easy to see that the sine terms in the product are equal to the eigenvalues in Eq. 37 up to a constant term. Consequently, 1/*λ*_{j} can be written as a product of all other eigenvalues *λ*_{k}, for *k* ≠ *j*. Still, implementing Eq. 38 would require *O*(*N*) qubits, since every inverse eigenvalue depends on 2^{n} − 1 others. The authors however used the periodicity of the discrete sines and some trigonometric relations to prove that 1/*λ*_{j} can be computed from Eq. 38 as product of only *n* − 2 sine terms. This allows the algorithm to be implemented in 3*n* qubits. Moreover, the sine expressions in Eq. 38 can be peformed straightforwardly as controlled *R*_{y} rotations.

Wang et al. also present the circuit for their algorithm, stating that the total cost expressed in one and two qubits gates is *R*_{y} operations could further reduce the gate count to 10*n*^{2}, by adding *n* − 2 ancillary qubits. For a few (*n* < 10) qubit problems the required number of gates may fit the specifications of with NISQ devices. Nevertheless, one should keep in mind that both gate count and the declared time complexity of *O*(poly log(1/*ɛ*)) do not take into account the preparation of

A different NISQ approach to the discrete Poisson problem is to use variational quantum algorithms to solve the underlying linear system (Liu H.-L. et al., 2021). In case of linear systems, the variational state is *U*(** θ**) is the ansatz (see Section 2.5) and the loss functions needs to be zero when

In the Poisson case the matrix *A* is symmetric, therefore *A*^{†}*A* = *A*^{2} and Eq. 39 becomes

A necessary condition for advantage with VQAs is that the terms in Eq. 40 must be efficiently evaluated on a quantum computer and ideally must be hard to evaluate in classically. Bravo-Prieto et al. (2019) proved that the latter requirement is satisfied for loss functions such as that in Eq. 39. On the other hand, efficient quantum evaluation is possible only if the following two necessary conditions are met.

1. *A* and *A*^{2} can be decomposed in *O*(poly log(*N*)) operators *O*_{k}

2. These operators are observables and have a simple tensor-product form.

Notice that we will talk here about observables as Hermitian operators, whose eigenvectors form an orthonormal basis for measurement.

Liu et al. show that *A* and *A*^{2} for the 1-dimensional Poisson matrix fulfill both requirements. For the number of decomposition terms, one can recursively decompose *A*_{m}, with *m* being the number of qubits as

where

For example, the 1 and 2 qubit cases are

where *A*_{m−1} with the identity and adding the center off-diagonal terms for the *m*th case, *A*_{m} can be written as a sum of 2*m* + 1 terms.

The matrix *A*^{2} can instead be split into the following two submatrices

The matrix *B*_{m} in Eq. 43 is decomposed in the same fashion as *A*_{m}, while *C*_{m} is the sum of just two terms

Overall, *m* + 1 terms. Since both *A* and *A*^{2} are decomposed in *O*(log(*N*)) operators, the first requirement for advantage is satisfied. Furthermore, Eq. 28 states that the *d*-dimensional Poisson equation is the sum of tensor products of the identity tensor with the one-dimensional matrix. Therefore, the *d*-dimensional Poisson problem with Dirichlet boundary conditions can be handled with the same operators as in the *d* = 1 case and specifically with *d*(2*m* + 1) of them for *A* and (*d*(2*m* + 1))^{2} for *A*^{2}. However, while Liu et al. also give a decomposition for the Neumann, Robin and mixed boundary conditions in the *d* = 1 case, the authors claim that the {*σ*_{+}, *σ*_{−}} decomposition cannot be extended trivially to the multidimensional case with boundary conditions different than Dirichlet ones.

Still, the operators in *A*_{m} and *B*_{m} are not Hermitian and cannot be used as observables. Liu et al. notice however that they can be made symmetric by mapping them in the higher dimensional space of Bell states. For instance, one can consider the 2 × 2, 1 qubit case and build two new operators *O*_{11} and *O*_{12}, such that

where

it is possible to evaluate the terms appearing in Eq. 40, such as

For higher matrix dimensions, measuring the expectation values of *A*_{m} and *B*_{m} requires to use similar operators as those in Eq. 45, but whose eigenvalues are entangled states in more than two qubits.

By measuring in the Bell basis, the second requirement is also satisfied. Indeed, Figure 5 shows that the measurement circuit in Liu H.-L. et al. (2021) is shallow and made by only one and two qubits gates. However, this circuit assumes full qubits connectivity, which is generally not available in current hardware. Therefore, the actual implementation on hardware will likely require *SWAP* gates to meet the design connectivity.

**FIGURE 5**. Circuit for measuring the state

A final interesting remark regards the ansatz *U*(** θ**) chosen by Liu et al. for their simulations. This is the quantum alternating operator ansatz (QAOA) (Farhi et al., 2014; Hadfield et al., 2019), which consists of a layered circuit, each layer having only two parameters, corresponding to the evolution times of mixer and driver Hamiltonians. Liu H.-L. et al. (2021) chose the two Hamiltonians, such that their gate depth grows only linearly with the number of qubits. Furthermore, the results obtained on a quantum simulator show that the number of QAOA layers only needs to increase linearly with the number of qubits for fixed solution fidelity, resulting in a circuit that is overall suitable for NISQ hardware.

#### 3.1.6 Quantum annealing

Srivastava and Sundararaghavan provided the first example of solving differential equations on a quantum annealer (Srivastava and Sundararaghavan, 2018). The motivation behind their work was that the functional minimization form of the differential problem can be written in terms of the discrete solution and solved as a combinatorial optimization problem. More precisely, the problem becomes a binary graph-coloring one, which is NP-hard.

As seen in Section 2.5.5, a quantum annealer finds the ground state of the Ising Hamiltonian, hereby reported for clarity

The *q*_{i} binary variables encode the values of the discrete problem variables, while *H*_{i} and *J*_{ij} depend on problem data and boundary conditions.

For instance, consider the 1D Laplace problem with Dirichlet boundary conditions

where *L* is the length of the domain. The associated energy potential is

One can replace *u*(*x*) with the discrete solution

where *φ*_{i}(*x*) can be taken as the linear finite element shape functions in Eq. 36. In particular, finite element approximations are local, which fits well the fact that quantum annealers are only locally connected.

Considering only two elements (3 nodes) on the unit domain and substituting Eq. 50 in Eq. 49 leads to the discrete functional

where\

Srivastava and Sundararaghavan propose to associate every element to a local graph, to repeat this unit throughout the annealer’s graph and to form the global connections according to the finite element connectivity (Srivastava and Sundararaghavan, 2018). Every node *i* is associated with 3 qubits

In this way, *a*_{i} can in principle assume 9 possible values. However, some of these values lead to invalid solutions and need to be penalized when writing the discrete potential.

Figure 6A shows the node graph used in Srivastava and Sundararaghavan (2018). This has associated *H*_{l} and *H*_{l} and

**FIGURE 6**. Generic node graph **(A)** and element graph **(B)** for the one dimensional Laplace equation in Srivastava and Sundararaghavan (2018) with the corresponding weigths of the Ising Hamiltonian.

Figure 6B shows instead the element graph for a single element, that is characterized by the matrix

where *i* and *j* are the nodes connected by element *n*.

Once *H*_{l},

However, if a higher precision is required while spanning the same range of possible values, then more than 3 qubits would be required per node. In a realistic FE model with thousands of nodes this would require a high number of qubits and a high degree of connectivity, which may be beyond hardware capacity. Therefore Srivastava and Sundararaghavan proposed the so-called “box algorithm”. The main idea is to have the values *r*, thus

In this way, the nodal values are spaced around *r* and for *N* nodes, all possible admissible values will be distributed on the surface of a *N*-dimensional box. The procedure consists in doing subsequent annealings until the desired precision *r* is met. After every annealing, *r* are allowed to vary according to the following logic.

1. If *i* has lower potential than

2. If *r* is reduced.

The procedure continues until *r* is below a threshold defined by the required solution precision. Also, Srivastava and Sundararaghavan (2018) show that this is a convergent process.

The graph-coloring and box algorithm is benchmarked against two truss problems, described by the equation

where *E*(*x*) and *A*(*x*) are respectively the distributed axial stiffness and area, while *f*(*x*) is the distributed load over the truss. The first example is a truss with a discontinues cross-section at half length, while the second one is a tapered truss with distributed load. Both examples converge to the numerical FE solution using up to 6 elements. The authors also mention that a finer discretization would likely require a two-point version of the box technique, in order to better exploit the connectivity of the annealer.

### 3.2 Heat equation

The heat equation is another widely applicable mathematical model in solid mechanics, which can be used to predict temprerature profiles and heat concentrations in structural components.

Let the spatial domain be *u*(**x**, *t*) satisfies the equation

and the differential problem is completed once the boundary and initial conditions are also specified.

For problems relevant in structural mechanics, the problems are at most three-dimensional. Also, one needs to account for the material’s thermal diffusivity, which is in general a tensor **x**. From here onwards however, the material is assumed isotropic and homogeneous, therefore the thermal diffusivity becomes a constant scalar *α* and

#### 3.2.1 Comparison of classical and gate-based quantum methods

A comprehensive study on quantum solvers for the heat equation was conducted by Linden et al. (Linden et al., 2020). Specifically, the authors compared 5 classical and 5 quantum methods in terms of their theoretical runtimes to approximate the temperature integrated over a region *S* ⊆Ω, up to precision *ɛ* with 99% success probability of the algorithm, i.e.

where

The discretization in time consists in first-order forward finite differences, while second-order centered finite differences are used for the space grid. This scheme is also called *forward time centered space* or FTCS, for short. Therefore Eq. 57 becomes the finite difference equation

where Δ*t* = *T*/*M*, *M* being the number of time intervals, while Δ*x* = *L*/*N*, assuming equal length *L* for each dimension and division in *N* intervals. Furthermore, *u* due to the finite difference discretization.

By grouping the values of

Also, by defining *t*_{k}, *k* = {1, …*M*}, Eq. 60 can be written in vector form as

where

The classical solvers studied by Linden for this problem are the following.

1. *Single linear system approach with conjugate gradient.* Eq. 61 can be seen as a unique linear system, for the solution at subsequent times. In fact, if **u**_{0} is the initial condition in discrete space, then

where

Linden et al. (2020) suggest to use the conjugate gradient method to solve the sparse linear system in Eq. 62 in linear time. More accurately however, one should resort to a variation of this method such as the biconjugate gradient stabilized method (BiCGSTAB) (van der Vorst, 1992), which still runs in linear time, in order to account for the asymmetry of *A*.

2. *Time-stepping from initial condition.* This is a matrix-vector multiplication problem. In fact, Eq. 61 can be expanded up to **u**_{0} as

Then, the approximate solution at time *t*_{k} = *k*Δ*t* is obtained by *k* successive multiplications of **u**_{0}.

3. *Time-stepping using the Fast Fourier Transform (FFT).* For the FTCS scheme, the matrix *d* dimensions has the expression (Linden et al., 2020)

where *I*_{N} is the identity matrix of dimension *N*. Also,

which is a circulant matrix and can therefore be diagonalized by the Discrete Fourier Transform *F*. Since the circulant matrices operate on different dimensions, the matrix *F*, i.e. *F*^{⊗d}. Furthermore, *H* has eigenvalues

where Λ is the diagonal matrix with the eigenvalues of

Numerically, this consists in doing the Fast Fourier Transform (FFT) of **u**_{0}, multiplying the resulting vector by the *k*th powers of the eigenvalues of

4. *Random walk.* Eq. 60 can be seen in terms of stochastic quantities. In fact, introducing

It is easy to see that if *s* ≤ 1/(2*d*), Eq. 69 can be thought of as a stochastic process. In fact, the temperature *d*-dimensional lattice with probabilities determined by *s*. In this sense the FTCS equation is a random walk, where the position is the approximate temperature

5. *Fast random walk.* The standard random walk samples for all *m* time steps, each sample requiring *O*(*d* log(*N*)) time, for a total time of *O*(*Md* log(*N*)). A speed-up can be achieved with respect to this standard technique. First, sample from the intial distribution in *O*(*d* log(*N*)) time and then compute the number of steps in each dimension *d* and the number of positive/negative increments in every dimension, by sampling two binomial distributions in time *O*(log(*M*)) (Linden et al., 2020). The improved runtime is *O*(*d*(log(*M*) + log(*N*))).

For their comparison, Linden et al. discussed how quantum subroutines could speed-up the classical numerical algorithms for the heat equation. In the case of quantum algorithms, the final solution state **u** of the FTCS equation. Starting from *S* ⊆ Ω can be calculated up to precision *ɛ* using numerical quadrature, i.e.

where *G* is the *d*-dimensional grid in the domain Ω, *w*(**x**) are weights that depend on the specific numerical quadrature scheme.

The 5 quantum algorithms considered in (Linden et al., 2020) are listed below. The reader may notice that this list does not include the time-optimal quantum algorithm for ordinary differential equations proposed by Berry et al. (2017). However Linden et al. state that such method may have a cost higher than using a QLSA, even though the dependency with *d* is better for the quantum ODE solver (Linden et al. (2020); Appendix A).

1. *Quantum linear solver.* A quantum algorithm for linear systems can be used to solve *A***u** = **r**, with matrix *A* defined in Eq. 63. Linden et al. utilized the algorithm in (Chakraborty et al., 2019), that is logarithmically dependent on precision. Even though this method requires the matrix to be Hermitian, it can still be applied to the FTCS-discretized heat equation by solving for

2. *Fast forwarded quantum walk method.* Quantum walks are a form of random walks that can be performed with unitary operations on quantum states (Ambainis, 2004). The first results in quantum walks showed how these could simulate random walks quadratically faster in the limit behavior (for a number of time-steps going to infinity). However, one is generally interested in the dynamics of the heat equation, other than its steady state. Using the quantum walk fast-forwarding algorithm of Apers and Sarlette (Apers and Sarlette, 2018), allows to have quadratic speed-up even when simulating a random walk at intermediate times, making it applicable to the heat equation.

3. *Coherent diagonalisation.* The operator

where it is assumed that

where Λ is the same matrix appearing in Eq. 67. Therefore

a. Prepare

b. Apply QFT.

c. Apply the Λ^{k} operator. This is not unitary in general, but it can be implemented by using an ancilla qubit, applying a rotation controlled on the

d. Apply the inverse QFT.

4. R*andom walk with amplitude estimation.* The random walk technique can be ran *O*(1/*ɛ*) faster using amplitude estimation. In fact, approximating ∫_{S}*u*(**x**)d**x** requires *O*(1/*ɛ*^{2}) repetitions of the classical random walk due to the Chernoff bound. However, assume there is a Boolean function *f*(*s*) such that

at the end of the random walk. If *f* can be encoded as an oracle, then quantum amplitude estimation can estimate Pr(*f*(*s*) = 1) in *O*(1/*ɛ*) time. This can be directly used to compute ∫_{S}*u*(**x**)d**x**.

5. *Fast random walk with amplitude estimation.* In the same way as in standard random walks, amplitude estimation can be applied to the fast random walk technique (see classical methods).

Table 2 shows the time complexities of all classical and quantum methods analyzed in Linden et al. (2020). In general, the classical FFT diagonalization technique scores the best complexity for the one-dimensional heat equation, while fast random walks with quantum amplitude estimation have the lowest runtime for *d* ≤ 2. However, both of these techniques work only for hyper-rectangular domains, for which the heat equation has an analytical solution in terms of Fourier components (Haberman, 2014). Still, the amplitudes of these modes are integrals, often computed numerically. Depending on the initial condition **u**_{0} and its Fourier decomposition, one may need to estimate a high number of integrals, in which case the methods for rectangular regions may still be meaningful.

**TABLE 2**. Runtime comparison of classical and quantum methods for solving the FTCS heat equation (Linden et al., 2020). All runtimes are expressed in terms of terms of spatial dimensions *d* and error *ɛ* on the estimated temperature integral. The

On the other hand, quantum algorithms can still be faster even on generic domains. In fact, except for *d* = 1 where the classical time-stepping technique has lowest runtime, standard random walks with quantum amplitude estimation are as fast (*d* = 2) or slightly faster (*d* = 3).

However, Table 2 also shows that no quantum exponential speed-up is possible. This happens even when some of the underlying quantum subroutines are ‘exponentially faster’ than their classical counterparts, as in the case of linear solvers. However, as showed by Montanaro for elliptical problems discretized by FEM, obtaining a scalar quantity requires *O*(poly(1/*ɛ*)) samples of the final quantum state (Montanaro and Pallister, 2016).

Finally, Table 2 shows that methods using a quantum algorithms for linear system are never faster than the best classical algorithm for *d* < 5, be it for rectangular or generic domains. Thus, one should be aware of this limitation, if aiming for a speed-up to solve the heat equation in a 3- or lower-dimensional space.

#### 3.2.2 Quantum annealing

Another way to use quantum computing for solving the heat equation was proposed by Pollachini et al. (Pollachini et al., 2021). Their approach involved the use of quantum annealing in a quantum-classical loop, allowing the method to be implemented on DWave quantum annealers and solving the heat equation on a 9 × 9 grid.

The equation considered is

which describes the steady state of the temperature distribution on a 2D domain with source term *f*(*x*_{1}, *x*_{2}). Eq. 74 is actually a diffusion-reaction equation instead of a heat equation, where *u*(*x*_{1}, *x*_{2}) is the equilibrium temperature and *f*(*x*_{1}, *x*_{2}) is a distributed heat flux. Dirichlet conditions were used at the boundary.

Eq. 74 is discretized using centered finite differences in the usual way and the problem reduces to solving a linear system *A***u** = **f**, where *G* is the space grid.

Quantum annealers can solve problems that can be expressed as QUBO problems. In the case of linear systems, the quadratic Hamiltonian

has a ground state corresponding to the solution **u** = *A*^{−1}**f** (O’Malley and Vesselinov, 2016; Rogers and Singleton, 2020). Also, *u*_{i} can be restricted to the range [ − *d*_{i}, 2*c*_{i} − *d*_{i}) using the following mapping

where *d*_{i} and *c*_{i} are user-defined real numbers and *u*_{i} is associated to the binary string *d*_{i} and *c*_{i}.

Eq. 76 can be substituted in Eq. 75, providing the Ising Hamiltonian

where *N* is the number of nodes in the grid. Once the ground state of this Hamiltonian is obtained after annealing, the inverse of the mapping in Eq. 76 allows to reconstruct the solution.

Furthermore, Pollachini et al. provided a strategy to keep their algorithm hardware-feasible even for large problem dimensions. In fact, they proposed to use the iterative block Gauss-Seidel method, which consists in iteratively solving *D* blocks of dimension *N*/*D* instead of one *N*-dimensional linear system. For instance, taking *D* = 2

Eq. 78 can be solved iteratively, by making an initial guess for **u**_{2} as *k* + 1,

The quantum annealer takes care of solving each lower-dimensional linear systems in Eq. 79 and Gauss-Seidel iterations are repeated until convergence.

As mentioned, Pollachini et al. ran their algorithm on both DWave 2000Q and DWave Advantage quantum annealers (D-Wave, 2022). The source term *f*(*x*_{1}, *x*_{2}) was taken randomly and Eq. 74 was discretized on a 11 × 11 grid, corresponding to 9 internal points and a linear system with 81 unknowns. This is to date one of the largest linear systems solved (at least partially) on quantum hardware.

One of the issues that arised in the computations is a flattening of the error curve for increasing iterations of the block Gauss Seidel solver. This was attributed to saturating the floating-point precision achievable by a fixed number of qubits *R*. Indeed the authors fixed this issue by progressively shrinking the *d*_{i} range in Eq. 76, matching the same convergence curve as the classical Gauss-Seidel algorithm. However, increasing the number of qubits *R* per interval did not benefit the solution’s precision, likely due to an increase in hardware noise.

Despite the approach of Pollachini et al. being hardware-ready and verified, it is unclear whether it may provide an advantage. Furthermore, their method can only be applied to the steady-state problem. However, an alternative to solve the time-dependent within quantum annealing would be to perform a semi-discretization in space and then use the algorithm for systems of linear ODEs proposed in Zanger et al. (2021).

### 3.3 Wave equation

A third classical PDE is the wave equation, which describes the propagation of a perturbation through a medium. In solid mechaics, the wave equation can be used to model vibrations in structures or seismic wave propagation in soil.

The wave equation is written as

where *c* is the wave propagation speed in the medium. As usual, the problem is fully determined once the solution at the boundary and two initial conditions are given.

Except for a few instances, where an analytical solution can be found through separation of variables or using the method of characteristics, the FDM and FEM are generally used to find an approximate solution of the wave equation. For instance, in the case of FDM, the Laplacian on the right hand side of Eq. 80 is discretized with centered finite differences and then a Runge-Kutta scheme allows to find the solution at subsequent time steps, starting from the initial conditions

The exact runtimes of these classical methods depend on the order of discretization *r* of the Laplacian, the choice of the time-stepping technique, etc. However, their asymptotic behavior is bounded from below as

Nevertheless, replacing classical linear algebra subroutines with quantum algorithms can remove the exponential dependency on *d* also when solving the wave equation. This was proved by Costa et al. (Costa et al., 2019), who showed how to turn Eq. 80 into a Schrödinger equation and solve it using Hamiltonian simulation.

For the sake of explanation, let the wave equation be one dimensional and take *c* = 1. As usual, the domain can be reduced to a grid of spacing Δ*x* on which the Laplacian on the right hand side of Eq. 80 can be approximated. The number of gridpoints used for the finite difference approximation determines the *order r* of the discrete Laplacian *L*^{(r)} and the discretization error *r* = 2, such that

Furthermore, if one sees the grid that discretizes Ω as a graph *G*_{Δx} of *x*^{(i)} and *x*^{(i+1)} − *x*^{(i)}, the discrete Laplacian can be thought of as a matrix *L*^{(r)}(*G*_{Δx}) defined on this graph.

Keeping *r* = 2, Eq. 80 becomes

where **u** = [*u*_{1}, *u*_{2}, *…*, *u*_{N}] and *N* is the number of vertices.

Now, assume that a matrix *B* exists, such that *BB*^{†} = *L*. One can then write

where **u**_{V} = **u** and **u**_{E} are additional variables associated to the edges of the graph *G*.

Deriving Eq. 84 with respect to time, one obtains

which shows that, if *BB*^{†} = *L*, then **u**_{V} both evolves according to the Schrödinger equation (Eq. 84) and it is the solution of the original wave equation.

For an order 2 Laplacian, the *B* matrix is the graph *signed incidence matrix*. If one assigns random orientations to the edges of graph *G*_{Δx}, then

where *W*_{ij} are weights assigned to the edges of the graph. For instance, in case of the second order Laplacian, the graph is unweigthed (*W*_{ij} = 1 *∀i*, *j*).

If the Laplacian has order *r* > 2, the graph theoretical interpretation of *B* is not as straightforward. Yet, Costa et al. (2019) discusses a general algebraic procedure to determine the incidence matrices for these higher order Laplacians and provides the entries for *B* and *L*^{(r)} up to order 10.

In order to solve 84, one can perform Hamiltonian simulation to an initial state and determine **u**_{V} at time *t*. In particular, Costa et al. employ the algorithm of Berry et al. (2015b) for sparse Hamiltonian simulation, that is optimal with respect sparsity, error and simulation time.

It is shown that Hamiltonian simulation for time *t* = *T* requires a number of gates *g* that is *l* is a characteristic domain dimension. Most importantly, these runtimes show no exponential dependency on the dimension *d*, even though they do not include the time required to sample the output. Still, unless the full Hilbert space needs to be sampled, the measurement step would not reintroduce the curse of dimensionality, but just a *O*(1/*ɛ*) factor.

Thus, as seen for the Poisson and heat equations, the speed-up with respect to classical numerical solvers is exponential for variable dimensions, but at most polynomial if the dimension is fixed. What is interesting to notice though is that the homogeneous wave equation has the “quantum-appealing” characteristics of being interpreted as a Schrödinger equation and solved via Hamiltonian simulation. The same authors of Costa et al. (2019) notice that if Eq. 80 was treated as a second-order ODE, rather than a Schrödinger equation on an extended Hilbert space, it could be solved using the algorithm of Berry et al. (2017), but that would result in a quadratic slowdown with respect to using Hamiltonian Simulation.

The work of Costa had an important follow-up in Suau et al. (2021), where the authors studied the implementation, number of gates and actual runtime of Costa’s wave equation solver. As a benchmark problem, they took the simplest case of a 1-dimensional wave equation with homogeneous Dirichlet boundary conditions on the end points. However, Suau’s implementation slightly deviates from the original wave equation solver, since the authors replaced the optimal-complexity Hamiltonian simulation algorithm of Berry et al. (2015b) with the more common Lie-Trotter-Suzuki (LTS) product formula (Lloyd, 1996; Berry et al., 2006).

Suau et al. compute the number of gates and runtimes required for implementing their wave equation algorithm, choosing the following gate set

where.

The runtime is obtained by converting the Hamiltonian simulation circuit to the gates in Eq. 87 and using the gate execution times provided by the manufacturer (IBM, 2019). A first interesting point is that the circuit in Suau et al. (2021) represents one of the few instances of a quantum PDE algorithm specified in terms ‘common’ gates (i.e. directly translatable to hardware).

The gate counts of Suau’s algorithm match the asymptotic (big-*O*) behavior of the wave equation solver for variable error, simulation time and number of gridpoints. What is most interesting however, is that the constants hidden in the big-O scaling are huge (10^{5}–10^{8}). This results in extremely high gate counts even just for solving the simplest possible instance of the wave equation. For instance, given a moderately interesting grid size of 10^{6} nodes, the quantum wave equation solver requires 10^{17} gates. Furthermore, the total runtime required for such a problem size is almost 1,000 calendar years (Suau et al., 2021). In terms of number of qubits, the solver requires roughly 70 logical qubits, where ‘logical’ means fully error corrected. As mentioned by Suau et al. (2021) however, full error correction for a logical qubit requires between 1,000 and 10,000 physical qubits, therefore vastly overshooting the NISQ hardware characteristics.

## 4 Nonlinear PDEs

Possibly the most computationally demanding tasks in computational mechanics are those related to solving nonlinear problems. The non linearities may be characteristics of the material, such as hysteresis, plasticity or damage, but also arise in case of large displacements or in the presence of contact. Whatever the cause, the classical numerical solution is generally iterative and consists in solving large linear system of equations in possibly many iterations.

Quantum algorithms for nonlinear PDEs are scarce up to present date, and no work focuses specifically on structural mechanics. However, Lubasch et al. (2020) and Kyriienko et al. (2021) both proposed techniques to solve generic (or quasi-generic) nonlinear PDEs. Both approaches consist in variationally training a parametrized circuit and on using a hybrid stratregy, whereby the quantum computer estimates the cost function terms and the classical one implements the optimization update. However, the two methods have substantial differences in how to encode the nonlinearities and on how to compute the cost function.

The algorithm of Lubasch is schematically represented in Figure 7. The encoding of nonlinear term in the cost function is performed via the *quantum nonlinear processing unit* (QNPU), which is a circuit meant to compute nonlinear functions of polynomial form that appear in the cost function. These can be written as

**FIGURE 7**. Scheme of the variational circuit in (Lubasch et al., 2020). The quantum nonlinear processing unit (QNPU) takes as input *r* copies of the solution vector, generated by the ansatz *O*_{j} in Eq. 92 as quantum operators. The circuit’s output comes from the measurements of the ancilla qubit and corresponds to the required cost function term.

Where the terms *u*^{(i)} are copies of the solution function and *u*^{(1)}. Furthermore, the *O*_{j} terms are different linear operators that are applied on different copies.

One can consider *u*^{(j)}. As usual, *O*_{j} and performs point-by-point multiplication depending on the specific nonlinear terms in the PDE.

The introduction of repeated input and the QNPU complicates the quantum circuit with respect to VQAs for linear problems. However, Lubasch et al. face the problem of circuit depth by encoding the quantum ansatz and the *O*_{j} operators as matrix product states (MPS) of bond dimension *χ*. This ensures that the circuits have *n* and *χ* (Lubasch et al., 2020), preventing any curse of dimensionality.

The combination of multiple inputs and QNPU allows for an efficient way to reproduce nonlinearities. However, the QNPU block is strictly problem-dependent and it may not be trivial to implement, depending on the nonlinear expression.

A more versatile technique in this sense was proposed by Kyriienko et al, who considered the generic nonlinear problem

here written for the 1D case.

Similarly to other near–term methods, this algorithm uses the variational principle to find the parametrized solution. However, a key difference in Kyriienko et al. (2021) is the fact that the solution is not represented as a discrete set of values on a grid, but as a function of *x*, thanks to the so–called *quantum feature map*. The concept of feature map originates from the machine learning literature and it consists in embedding the data into the model as parameters. Translating to quantum circuits, this means that *x* can be mapped to a 2^{n}-dimensional space with a unitary operator *U*_{ξ}(*x*) parametrized through a nonlinear function *ξ*(*x*). Possibly the easiest instance of quantum feature map is

where *Y*_{i} is the Pauli Y gate applied to qubit *i*. A common choice is *ξ*(*x*) = arcsin(*x*), which means that *x* will be encoded in quantum amplitudes that are polynomials up to order 2^{n} in

As mentioned, the best approximation of *u*(*x*) is found variationally using the quantum circuit model. Therefore, Kyriienko’s technique also makes use of a circuit *U*_{θ}, whose parameters ** θ** are varied to minimized an appropriate loss function. Overall, the parametrized quantum state whose amplitudes embed the tentative solution is

In order to map between *u*(*x*), one needs also to specify an observable

Once, the parameters ** θ** have been optimized, one can then reconstruct the approximate solution at a specific point

*x*, by simply measuring the expectation value of

^{n}-dimensional Hilbert space, as it is necessary with the quantum algorithm based on amplitude encoding. Also, most of structural mechanics highly nonlinear problems are interesting with respect to the value of the solution or functions of it in just a small subset of points. For instance, this could be the case of probing the stress field at the crack tip of a fractured solid, by constructing

In general, the training of *U*_{θ} requires the minimization of a loss function ** θ**. Since the feature map is parametrized on

*x*, the derivative

*m*times (Crooks, 2019; Mari et al., 2021). For

*m*= 1

where “+” and “−” symbolically represent the positive/negative *x*-shift and the sum is among all the parametrized gates composing the feature map. Suitable choices for *L*_{θ} can be the residual in Eq. 93 or the mean squared difference with respect to the exact solution, if this is available.

One clearly sees a parallel between Kyriienko’s method and training of neural networks. In the PDE case, during training, the function *u* is evaluated on a grid, which can be considered as the *training dataset* of a machine learning routine. Afterwards, the solution *test dataset*, in order to assess the validity of the model.

## 5 Discussion

The previous sections reviewed the literature of partial differential equations pertinent to structural mechanics. This analysis was divided into linear and nonlinear PDEs, which is a standard classification for differential problems. The first group includes the works about Poisson, heat and wave equations, while the second one deals with the methods to solve general nonlinear problems.

Linear problems can be solved using all different quantum paradigms, i.e. full gate-based, hybrid quantum computing and quantum annealing. In terms of the full-quantum gate-based primitives, such as quantum linear solvers, quantum Hamiltonian simulation etc, the Poisson, heat and wave equations can be solved with these quantum algorithms and inherit their complexities. However, the different character of the equations and the specific discretization determine which quantum routines are applicable and the extent of the advantage. For instance, quantum linear solvers are applicable to every linear PDE, both stationary and time-dependent, if the latter are written as a single linear system spanning multiple time steps. Nevertheless, the Poisson equation (on rectangular domains) with periodic boundary conditions is a favourite candidate for QLSAs, since the finite difference approximation of the Laplacian results in a circulant matrix that is diagonalized by the QFT (Cao et al., 2013; Wang et al., 2020c,b; Childs et al., 2021). This allows to do Hamiltonian simulation in the QLSA solving the Poisson equation exponentially faster than with non-circulant matrices.

On the other hand, heat and wave equations benefit more from different quantum subroutines. One hint to this is the evolutionary character of both equations, which means that linear system dimensions scale multiplicatively with respect to time grid size. Also, the time dependency seems to suggest the approach of semi-discretizing in space and then solving systems of ODEs with some Hamiltonian simulation algorithm. Indeed, this approach is ideal for the wave equation where Hamiltonian simulation solves the related graph problem in the higher dimensional space (Costa et al., 2019). On the other hand, the heat equation does not to benefit as much from quantum ODE solvers and the useful analysis of Linden et al. (2020) proves that minimum runtimes are achieved when accelerating a classical method (classical random walks) with amplitude acceleration.

Of course, all previous considerations hold for quantum subroutines that require error-corrected hardware. Still, near-term quantum techniques for linear PDEs exist under the umbrella of quantum annealing and variational quantum computing, even though the efforts in this sense are in their infancy. For PDEs in structural mechanics, only two quantum anneling algoritms have been proposed, namely for elliptical FE problems (Srivastava and Sundararaghavan, 2018) and for the stationary heat equation (Pollachini et al., 2021). Clearly, a first gap in this branch of literature are quantum annealing algorithms for evolutionary problems, such as heat and wave equations.

Also VQAs have just recently been applied to linear PDEs. Of course, the general literature on variational quantum computing is vast, but their use for PDEs have been limited to generic nonlinear problems (Lubasch et al., 2020; Kyriienko et al., 2021). However, there is still a lack of works for specific PDEs, even though such specialization is critical. For instance, Liu H.-L. et al. (2021) showed the relevance of the Poisson matrix in the context of VQAs. In fact, the discrete Poisson matrix can be decomposed in a poly-logarithmic number of observables, which is a necessary condition for advantage of variational PDE solvers.

For nonlinear problems, the matter of choosing a quantum primitive is not as straightforward, because all quantum operations are ultimately linear. The ways forward seem ultimately two, i.e linearization and application of existing quantum techniques or variational algorithms. Even though there exist research on linearization of nonlinear ODEs and use of QLSAs at each step (Liu J.-P. et al., 2021), all works on nonlinear PDEs relied so-far only on variational quantum computing. As it seems likely, this paradigm together with quantum annealing will likely be the only quantum alternative viable in the near term to solve PDEs, linear and nonlinear alike. However, there is currently a gap in literature about what can be expected from quantum computation for nonlinear PDEs in the error-corrected era. Having more insight in this direction would be extremely valuable, since nonlinear problems represent the most expensive and thus most interesting problems from a computational standpoint.

A final remark concerns extent of the overlap between quantum PDE and structural mechanics literature. In essence, this is currently limited to a single work on quantum annealing for truss problems (Srivastava and Sundararaghavan, 2018). In fact, the vast majority of literature focuses on ‘academic’ PDEs (Poisson, heat and wave) on hypercubic domains and for favorable boundary conditions. The exceptions are the methods for nonlinear equations, which essentially provide a framework from data encoding to solution, but ultimately leave the choice of critical hyperparameters to the user (Lubasch et al., 2020; Kyriienko et al., 2021). Of course, both sides of the spectrum project onto structural mechanics, but quantum computation still has to be tested against the specific and interesting problems in this field. What is even more surprising is that other disciplines away from quantum physics, yet heavily relying on numerical calculus (fluid mechanics, finance, etc) already applied quantum algorithms to their own cost-intensive problems. For instance, several works in fluid mechanics field used quantum subroutines to solve both the lattice Boltzmann (Mezzacapo et al., 2015; Todorova and Steijl, 2020; Budinski, 2021a) and the Navier-Stokes (Steijl and Barakos, 2018; Gaitan, 2020; Budinski, 2021b; Gaitan, 2021) equations. The hope is that structual mechanics will also explore the use of quantum algorithms to support expensive simulations, such as those involving material nonlinearities and large structural deformations.

## Author contributions

GT reviewed the papers discussed in this work and wrote the manuscript. BC, MM, and MG all reviewed the work and provided critical feedback. RB followed the progress of the work and gave approval for publication.

## Funding

This review work was provided by the Delft University of Technology.

## Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## Footnotes

^{1}The “big Omega” notation expresses the asymptotic lower bound on the number of computational resources.

^{2}A state or a subspace corresponding to a subset of basis states. For example, given the state

## References

Aharonov, D., Jones, V., and Landau, Z. (2009). A polynomial quantum algorithm for approximating the jones polynomial. *Algorithmica* 55, 395–421. doi:10.1007/s00453-008-9168-0

Aharonov, D., and Ta-Shma, A. (2003). “Adiabatic quantum state generation and statistical zero knowledge,” in Proceedings of the Thirty-Fifth Annual ACM Symposium on Theory of Computing, San Diego CA USA, June 9 - 11, 2003. doi:10.1145/780542.780546

Ambainis, A. (2004). Quantum walks and their algorithmic applications. *Int. J. Quantum Inf.* 1, 507–518. doi:10.1142/s0219749903000383

Ambainis, A. (2012). “Variable time amplitude amplification and quantum algorithms for linear algebra problems,” in 29th International Symposium on Theoretical Aspects of Computer Science (STACS 2012), Paris, France, February 29th - March 3rd, 2012. doi:10.4230/LIPIcs.STACS.2012.636

An, D., and Lin, L. (2019). *Quantum linear system solver based on time-optimal adiabatic quantum computing and quantum approximate optimization algorithm*. United States: arXiv:1909.05500.

An, D., Linden, N., Liu, J.-P., Montanaro, A., Shao, C., Wang, J., et al. (2020). Quantum-accelerated multilevel Monte Carlo methods for stochastic differential equations in mathematical finance. *Quantum* 5, 481. doi:10.22331/q-2021-06-24-481

Apers, S., and Sarlette, A. (2018). *Quantum fast-forwarding: Markov chains and graph property testing*. United States: arXiv. *1804.02321*.

Benedetto, F. D. (1997). Preconditioning of block toeplitz matrices by sine transforms. *SIAM J. Sci. Comput.* 18, 499–515. doi:10.1137/s1064827595258335

Benzi, M., and Tûma, M. (1999). A comparative study of sparse approximate inverse preconditioners. *Appl. Numer. Math.* 30, 305–340. doi:10.1016/s0168-9274(98)00118-4

Berry, D. W., Ahokas, G., Cleve, R., and Sanders, B. C. (2006). Efficient quantum algorithms for simulating sparse Hamiltonians. *Commun. Math. Phys.* 270, 359–371. doi:10.1007/s00220-006-0150-x

Berry, D. W., Childs, A. M., Cleve, R., Kothari, R., and Somma, R. D. (2014). “Exponential improvement in precision for simulating sparse Hamiltonians,” in Proceedings of the Forty-Sixth Annual ACM Symposium on Theory of Computing, New York, NY, USA, 31 May 2014- 3 June 2014, 283–292. doi:10.1145/2591796.2591854

Berry, D. W., Childs, A. M., Cleve, R., Kothari, R., and Somma, R. D. (2015a). Simulating Hamiltonian dynamics with a truncated Taylor series. *Phys. Rev. Lett.* 114, 090502. doi:10.1103/PhysRevLett.114.090502

Berry, D. W., Childs, A. M., and Kothari, R. (2015b). “Hamiltonian simulation with nearly optimal dependence on all parameters,” in IEEE 56th Annual Symposium on Foundations of Computer Science, Berkeley, CA, USA, Oct. 17 2015 to Oct. 20 2015. doi:10.1109/focs.2015.54

Berry, D. W., Childs, A. M., Ostrander, A., and Wang, G. (20172017). Quantum algorithm for linear differential equations with exponentially improved dependence on precision. *Commun. Math. Phys.* 356, 1057–1081. doi:10.1007/s00220-017-3002-y

Brassard, G., Hoyer, P., Mosca, M., and Tapp, A. (2000). AMS Contemporary Mathematics in *Quantum computation and quantum information*. Editor J. SamuelLomonaco (Cambridge, New York: Cambridge University Press), 305, 53–74. doi:10.1090/conm/305/05215*Quantum amplitude amplification and estimation*

Bravo-Prieto, C., LaRose, R., Cerezo, M., Subasi, Y., Cincio, L., and Coles, P. J. (2019). *Variational quantum linear solver*. United States: arXiv:1909.05820v2. *1909*.

Brenner, S. C., and Scott, L. R. (2010). “The mathematical theory of finite element methods. No. 15,” in *Texts in applied mathematics* (New York, NY: Springer), 3. nachdr.

Budinski, L. (2021a). Quantum algorithm for the advection-diffusion equation simulated with the lattice Boltzmann method. *Quantum Inf. process.* 20, 57. doi:10.1007/s11128-021-02996-3

Budinski, L. (2021b). *Quantum algorithm for the Navier-Stokes equations*. United States: rXiv:2103.03804. *2103*.

Buhrman, H., Cleve, R., Watrous, J., and de Wolf, R. (2001). Quantum fingerprinting. *Phys. Rev. Lett.* 87, 167902. doi:10.1103/PhysRevLett.87.167902

Cao, Y., Papageorgiou, A., Petras, I., Traub, J., and Kais, S. (2013). Quantum algorithm and circuit design solving the Poisson equation. *New J. Phys.* 15, 013021. doi:10.1088/1367-2630/15/1/013021

Cerezo, M., Arrasmith, A., Babbush, R., Benjamin, S. C., Endo, S., Fujii, K., et al. (2021a). Variational quantum algorithms. *Nat. Rev. Phys.* 3, 625–644. doi:10.1038/s42254-021-00348-9

Cerezo, M., Sone, A., Volkoff, T., Cincio, L., and Coles, P. J. (2021b). Cost function dependent barren plateaus in shallow parametrized quantum circuits. *Nat. Commun.* 12, 1791. doi:10.1038/s41467-021-21728-w

Chakrabarti, S., Krishnakumar, R., Mazzola, G., Stamatopoulos, N., Woerner, S., Zeng, W. J., et al. (2021). A threshold for quantum advantage in derivative pricing. *Quantum* 5, 463. doi:10.22331/q-2021-06-01-463

Chakraborty, S., Gilyén, A., and Jeffery, S. (2019). “The power of block-encoded matrix powers: Improved regression techniques via faster Hamiltonian simulation,” in 46th International Colloquium on Automata, Languages, and Programming (ICALP 2019), Patras, Greece, 10 JULY 2019. doi:10.4230/LIPIcs.ICALP.2019.33

Chen, Z.-Y., Zhou, Q., Xue, C., Yang, X., Guo, G.-C., Guo, G.-P., et al. (2018). 64-qubit quantum circuit simulation. *Sci. Bull.* 63, 964–971. doi:10.1016/j.scib.2018.06.007

Childs, A. M., Kothari, R., and Somma, R. D. (2017). Quantum algorithm for systems of linear equations with exponentially improved dependence on precision. *SIAM J. Comput.* 46, 1920–1950. doi:10.1137/16m1087072

Childs, A. M., Liu, J.-P., and Ostrander, A. (2021). High-precision quantum algorithms for partial differential equations. *Quantum* 5, 574. doi:10.22331/q-2021-11-10-574

Ciliberto, C., Herbster, M., Ialongo, A. D., Pontil, M., Rocchetto, A., Severini, S., et al. (2018). Quantum machine learning: A classical perspective. *Proc. R. Soc. A* 474, 20170551. doi:10.1098/rspa.2017.0551

Clader, B. D., Jacobs, B. C., and Sprouse, C. R. (2013). Preconditioned quantum linear system algorithm. *Phys. Rev. Lett.* 110, 250504. doi:10.1103/PhysRevLett.110.250504

Costa, P. C. S., Jordan, S., and Ostrander, A. (2019). Quantum algorithm for simulating the wave equation. *Phys. Rev. A . Coll. Park.* 99, 012323. doi:10.1103/physreva.99.012323

Crooks, G. E. (2019). *Gradients of parameterized quantum gates using the parameter-shift rule and gate decomposition*. United States: arXiv:1905.13311. *1905*.

D-Wave (2022). *D-wave system documentation*. Burnaby, Canada: D-Wave Systems Inc. Available at: https://docs.dwavesys.com/docs/latest/index.php.

Damianou, P. A. (2014). A beautiful sine formula. *Am. Math. Mon.* 121, 120. doi:10.4169/amer.math.monthly.121.02.120

Farhi, E., Goldstone, J., and Gutmann, S. (2014). *A quantum approximate optimization algorithm*. United States: arXiv. *1411.4028*.

Fontanela, F., Jacquier, A., and Oumgari, M. (2021). *A quantum algorithm for linear pdes arising in finance*. United States: arXiv. *1912.02753*.

Gaitan, F. (2020). Finding flows of a Navier–Stokes fluid through quantum computing. *npj Quantum Inf.* 6, 61. doi:10.1038/s41534-020-00291-0

Gaitan, F. (2021). Finding solutions of the Navier-Stokes equations through quantum computing—Recent progress, a generalization, and next steps forward. *Adv. Quantum Technol.* 4, 2100055. doi:10.1002/qute.202100055

Giovannetti, V., Lloyd, S., and Maccone, L. (2008a). Architectures for a quantum random access memory. *Phys. Rev. A . Coll. Park.* 78, 052310. doi:10.1103/physreva.78.052310

Giovannetti, V., Lloyd, S., and Maccone, L. (2008b). Quantum random access memory. *Phys. Rev. Lett.* 100, 160501. doi:10.1103/physrevlett.100.160501

Golub, G. H., and van Loan, C. F. (2013). *Matrix computations*. fourth edn. Baltimore, Maryland, United States: JHU Press.

Grover, L. K. (1996). “A fast quantum mechanical algorithm for database search,” in Annual ACM Symposyum on Theory Of Computing, Rome, Italy, June 20-24, 2022, 212–219.

Grover, L. K. (19971997). Quantum mechanics helps in searching for a needle in a haystack. *Phys. Rev. Lett.* 79, 325–328. doi:10.1103/PhysRevLett.79.325

Haberman, R. (2014). *Applied partial differential equations with fourier series and boundary value problems*. Harlow: Pearson.

Hadfield, S., Wang, Z., O’Gorman, B., Rieffel, E., Venturelli, D., Biswas, R., et al. (2019). From the quantum approximate optimization algorithm to a quantum alternating operator ansatz. *Algorithms* 12, 34. doi:10.3390/a12020034

Harrow, A. W., Hassidim, A., and Lloyd, S. (2009). Quantum algorithm for linear systems of equations. *Phys. Rev. Lett.* 103, 150502. doi:10.1103/physrevlett.103.150502

Heim, N., Ghosh, A., Kyriienko, O., and Elfving, V. E. (2021). *Quantum model-discovery*. United States: arXiv. *2111.06376*.

Huang, H.-Y., Bharti, K., and Rebentrost, P. (2021). Near-term quantum algorithms for linear systems of equations with regression loss functions. *New J. Phys.* 23, 113021. doi:10.1088/1367-2630/ac325f

Huang, Y., and McColl, W. F. (1997). Analytical inversion of general tridiagonal matrices. *J. Phys. A Math. Gen.* 30, 7919–7933. doi:10.1088/0305-4470/30/22/026

Kandala, A., Mezzacapo, A., Temme, K., Takita, M., Brink, M., Chow, J. M., et al. (2017). Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets. *Nature* 549, 242–246. doi:10.1038/nature23879

Kato, T. (1950). On the adiabatic theorem of quantum mechanics. *J. Phys. Soc. Jpn.* 5, 435–439. doi:10.1143/jpsj.5.435

Kingma, D. P., and Ba, J. (2015). “Adam: A method for stochastic optimization,” in 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015.

Kyriienko, O., Paine, A. E., and Elfving, V. E. (2021). Solving nonlinear differential equations with differentiable quantum circuits. *Phys. Rev. A. Coll. Park.* 103, 052416. doi:10.1103/physreva.103.052416

Linden, N., Montanaro, A., and Shao, C. (2020). *Quantum vs. classical algorithms for solving the heat equation*. United States: arXiv:2004.06516v2. *2004*.

Liu, H.-L., Wu, Y.-S., Wan, L.-C., Pan, S.-J., Qin, S.-J., Gao, F., et al. (2021a). Variational quantum algorithm for the Poisson equation. *Phys. Rev. A* 104, 022418. doi:10.1103/PhysRevA.104.022418

Liu, J.-P., Kolden, H. Ø., Krovi, H. K., Loureiro, N. F., Trivisa, K., Childs, A. M., et al. (2021b). Efficient quantum algorithm for dissipative nonlinear differential equations. *Proc. Natl. Acad. Sci. U. S. A.* 118, e2026805118. doi:10.1073/pnas.2026805118

Lloyd, S. (1996). Universal quantum simulators. *Science* 273, 1073–1078. doi:10.1126/science.273.5278.1073

Low, G. H., and Chuang, I. L. (2019). Hamiltonian simulation by qubitization. *Quantum* 3, 163. doi:10.22331/q-2019-07-12-163

Lubasch, M., Joo, J., Moinier, P., Kiffner, M., and Jaksch, D. (2020). Variational quantum algorithms for nonlinear problems. *Phys. Rev. A . Coll. Park.* 101, 010301. doi:10.1103/physreva.101.010301

Mari, A., Bromley, T. R., and Killoran, N. (2021). Estimating the gradient and higher-order derivatives on quantum hardware. *Phys. Rev. A. Coll. Park.* 103, 012405. doi:10.1103/physreva.103.012405

McClean, J. R., Boixo, S., Smelyanskiy, V. N., Babbush, R., and Neven, H. (2018). Barren plateaus in quantum neural network training landscapes. *Nat. Commun.* 9, 4812. doi:10.1038/s41467-018-07090-4

McGeoch, C. C. (2020). Theory versus practice in annealing-based quantum computing. *Theor. Comput. Sci.* 816, 169–183. doi:10.1016/j.tcs.2020.01.024

Mezzacapo, A., Sanz, M., Lamata, L., Egusquiza, I. L., Succi, S., Solano, E., et al. (2015). Quantum simulator for transport phenomena in fluid flows. *Sci. Rep.* 5, 13153. doi:10.1038/srep13153

Mocz, P., and Szasz, A. (2021). Toward cosmological simulations of dark matter on quantum computers. *ApJ.* 910, 29. doi:10.3847/1538-4357/abe6ac

Montanaro, A., and Pallister, S. (2016). Quantum algorithms and the finite element method. *Phys. Rev. A. Coll. Park.* 93, 032324. doi:10.1103/physreva.93.032324

Nielsen, M. A., and Chuang, I. L. (2010). *Quantum computation and quantum information*. 10th anniversary ed edn. Cambridge; New York: Cambridge University Press.

O’Malley, D., and Vesselinov, V. V. (2016). Toq.jl: A high-level programming language for d-wave machines based on julia. In 2016 IEEE High Performance Extreme Computing Conference (HPEC), Waltham, MA, 13-15 Sept. 2016 1–7. doi:10.1109/HPEC.2016.7761616

Patil, H., Wang, Y., and Krstić, P. S. (2022). Variational quantum linear solver with a dynamic ansatz. *Phys. Rev. A . Coll. Park.* 105, 012423. doi:10.1103/physreva.105.012423

Pollachini, G. G., Salazar, J. P. L. C., Góes, C. B. D., Maciel, T. O., and Duzzioni, E. I. (2021). Hybrid classical-quantum approach to solve the heat equation using quantum annealers. *Phys. Rev. A . Coll. Park.* 104, 032426. doi:10.1103/physreva.104.032426

Preskill, J. (2018). Quantum computing in the NISQ era and beyond. *Quantum* 2, 79. doi:10.22331/q-2018-08-06-79

Rogers, M. L., and Singleton, R. L. (2020). Floating-point calculations on a quantum annealer: Division and matrix inversion. *Front. Phys.* 8. doi:10.3389/fphy.2020.00265

Shewchuk, J. R. (1994). *An introduction to the conjugate gradient method without the agonizing pain*. USA: Tech. rep.

Somma, R. D., and Subaşı, Y. (2021). Complexity of quantum state verification in the quantum linear systems problem. *PRX Quantum* 2, 010315. doi:10.1103/prxquantum.2.010315

Spall, J. (1992). Multivariate stochastic approximation using a simultaneous perturbation gradient approximation. *IEEE Trans. Autom. Contr.* 37, 332–341. doi:10.1109/9.119632

Srivastava, S., and Sundararaghavan, V. (2018). Box algorithm for the solution of differential equations on a quantum annealer. *Phys. Rev. A* 99, 052355. doi:10.1103/PhysRevA.99.052355

Steijl, R., and Barakos, G. N. (2018). Parallel evaluation of quantum algorithms for computational fluid dynamics. *Comput. Fluids* 173, 22–28. doi:10.1016/j.compfluid.2018.03.080

Stokes, J., Izaac, J., Killoran, N., and Carleo, G. (2020). Quantum natural gradient. *Quantum* 4, 269. doi:10.22331/q-2020-05-25-269

Suau, A., Staffelbach, G., and Calandra, H. (2021). Practical quantum computing: Solving the wave equation using a quantum approach. *ACM Trans. Quantum Comput.* 2, 1–35. doi:10.1145/3430030

Subaşı, Y., Somma, R. D., and Orsucci, D. (2019). Quantum algorithms for systems of linear equations inspired by adiabatic quantum computing. *Phys. Rev. Lett.* 122, 060504. doi:10.1103/physrevlett.122.060504

Taube, A. G., and Bartlett, R. J. (2006). New perspectives on unitary coupled-cluster theory. *Int. J. Quantum Chem.* 106, 3393–3401. doi:10.1002/qua.21198

Todorova, B. N., and Steijl, R. (2020). Quantum algorithm for the collisionless Boltzmann equation. *J. Comput. Phys.* 409, 109347. doi:10.1016/j.jcp.2020.109347

Trotter, H. F. (1959). On the product of semi-groups of operators. *Proc. Am. Math. Soc.* 10, 545–551. doi:10.1090/s0002-9939-1959-0108732-6

van der Vorst, H. A. (1992). Bi-Cgstab: A fast and smoothly converging variant of bi-cg for the solution of nonsymmetric linear systems. *SIAM J. Sci. Stat. Comput.* 13, 631–644. doi:10.1137/0913035

Wang, S., Wang, Z., Li, W., Fan, L., Cui, G., Wei, Z., et al. (2020b). *A quantum Poisson solver implementable on nisq devices*. United States: arXiv. *2005.00256*.

Wang, S., Wang, Z., Li, W., Fan, L., Cui, G., Wei, Z., et al. (2020a). Quantum circuits design for evaluating transcendental functions based on a function-value binary expansion method. *Quantum Inf. Process.* 19, 347. doi:10.1007/s11128-020-02855-7

Wang, S., Wang, Z., Li, W., Fan, L., Wei, Z., Gu, Y., et al. (2020c). Quantum fast Poisson solver: The algorithm and complete and modular circuit design. *Quantum Inf. process.* 19, 170. doi:10.1007/s11128-020-02669-7

Xu, X., Sun, J., Endo, S., Li, Y., Benjamin, S. C., Yuan, X., et al. (2021). Variational algorithms for linear algebra. *Sci. Bull.* 66, 2181–2188. doi:10.1016/j.scib.2021.06.023

Keywords: quantum computing, partial differential equations, quantum algorithms, linear, nonlinear, advantage, near-term

Citation: Tosti Balducci G, Chen B, Möller M, Gerritsma M and De Breuker R (2022) Review and perspectives in quantum computing for partial differential equations in structural mechanics. *Front. Mech. Eng* 8:914241. doi: 10.3389/fmech.2022.914241

Received: 06 April 2022; Accepted: 30 June 2022;

Published: 01 September 2022.

Edited by:

Rene Steijl, University of Glasgow, United KingdomReviewed by:

Luca Dede’, Politecnico di Milano, ItalySachin Satish Bharadwaj, New York University, United States

Copyright © 2022 Tosti Balducci, Chen, Möller, Gerritsma and De Breuker. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Boyang Chen, b.chen-2@tudelft.nl