A quantum computing concept for 1-D elastic wave simulation with exponential speedup

SUMMARY Quantum computing has attracted considerable attention in recent years because it promises speedups that conventional supercomputers cannot offer, at least for some applications. Though existing quantum computers are, in most cases, still too small to solve significant problems, their future impact on domain sciences is already being explored now. Within this context, we present a quantum computing concept for 1-D elastic wave propagation in heterogeneous media with two components: a theoretical formulation and an implementation on a real quantum computer. The method rests on a finite-difference approximation, followed by a sparsity-preserving transformation of the discrete elastic wave equation to a Schr¨odinger equation, which can be simulated directly on a gate-based quantum computer. An implementation on an error-free quantum simulator verifies our approach and forms the basis of numerical experiments with small problems on the real quantum computer IBM Brisbane. The latter produce simulation results that qualitatively agree with the error-free version but are contaminated by quantum decoherence and noise effects. Complementing the discrete transformation to the Schr ¨ odinger equation by a continuous version allows the replacement of finite differences by other spatial discretization schemes, such as the spectral-element method. Anticipating the emergence of error-corrected quantum chips, we analyze the computational complexity of the best quantum simulation algorithms for future QCs. This analysis suggests that our quantum computing approach may lead to wave field simulations that run exponentially faster than simulations on classical computers.


INTRODUCTION
The progress of seismological research has always been closely linked to advances in computer technology.As early as 1968, the thennovel IBM 360/65 enabled the random sampling of more than 200,000 seismic models per hour, leading to one of the first Monte Carlo inversions for whole-Earth structure (Press, 1968).The Connection Machine CM-5 was the workhorse for large-scale finite-difference wave field simulations (Igel et al., 1995) and early 2-D acoustic full-waveform inversion (Igel et al., 1996), which could be extended to massive 3-D elastic inversions when GPUs became available (e.g., Bozdag et al., 2016).Novel chip designs that integrate CPUs and GPUs substantially reduce code complexity, thereby facilitating the efficient development and execution of seismic applications on personal computers (Gebraad and Fichtner, 2023).The list of prominent examples is endless.
Emerging quantum computers (QCs) are opening a new chapter in scientific computing, even though available QCs are still relatively small and inefficient compared to the latest supercomputers that are based on large numbers of CPUs and GPUs.The current stage of the QC development is called Noisy Intermediate Scale Quantum (NISQ; coined by Preskill, 2018), since their efficiency is mainly hampered by elevated noise (e.g., Knill, 2005;Nielsen and Chuang, 2010;Preskill, 2018).Recent advances in quantum error correction (e.g., Sivak et al., 2023;Kim et al., 2023;Bravyi et al., 2023) and the design of larger, more modular QCs (e.g., Gibney, 2019;Gambetta, 2020;Sevilla and Riedel, 2020) are promising developments, suggesting the feasibility of quantum computation for numerous algorithmic challenges within the next decade.Presently, QCs are already available in cloud environments, like IBM Quantum, and can be accessed freely for experimental use.This accessibility facilitates the testing and validation of novel quantum algorithms.

arXiv:2312.14747v2 [physics.geo-ph] 7 May 2024
QCs harness quantum phenomena such as entanglement and superposition for problem-solving algorithms, thereby promising up to exponential reductions in both runtime and memory complexity compared to some classical algorithms (e.g., Nielsen and Chuang, 2010;Benenti et al., 2004).As a result, quantum computation could make it feasible to tackle classes and sizes of problems considered intractable on traditional computing hardware (e.g., Shor, 1994;Harrow et al., 2009;Montanaro, 2016).As of today seismic imaging and related uncertainty quantification problems continue to be primarily limited by computational resources.It is therefore natural and timely to ask whether the potential of QCs can be harnessed for such inverse problems (Moradi et al. (2018(Moradi et al. ( , 2019))).Recent studies have used adiabatic quantum computers (O'Malley (2018); Golden and O'Malley (2021); Souza et al. (2022); Dukalski et al. (2023)) and variational quantum algorithms (Trahan et al. (2023)) for solving geoscience-related problems (see Section 5.5 for more details).In all studies, the quantum resources are used to solve linearised inverse problems.Here, we intend to take a first step toward non-linear Full Wave Form Inversion (FWI) (Virieux and Operto (2009); Fichtner (2010); Liu and Gu (2012)) by tackling the forward problem of solving the elastic wave equation, using a different algorithmic approach, known as Hamiltonian simulation (HS).HS has been invented to simulate quantum systems and is among the most promising applications of QCs with the potential for exponential speedup over classical computers (Low and Chuang (2019)).It has recently been demonstrated in a handful of studies that classical wave problems can be mapped to a Schrödinger equation so that HS can be harnessed to simulate Maxwell's equations (Jin et al. (2022(Jin et al. ( , 2023))), homogeneous acoustic waves, and mass-spring models (Suau et al. (2021); Costa et al. (2019); Babbush et al. (2023)).This study exploits HS to solve the 1-D elastic wave equation in heterogeneous media.Additionally, we perform initial experiments on actual QCs.
More specifically, we extend numerical wave propagation via HS (Costa et al., 2019;Suau et al., 2021) to heterogeneous media, a necessary step toward FWI on QCs.For this, we employ an analytical Cholesky decomposition, which facilitates the transformation of the elastic wave equation with heterogeneous medium properties to a Schrödinger equation.We proceed with the implementation of our method on the QC IBM Brisbane, in order to solve a small-scale problem and explain current challenges.These include systematic errors caused by noise and decoherence rates, which are characteristic for current NISQ QCs.Anticipating future developments, we finally delve into the potential scalability of an optimal HS algorithm on future quantum hardware equipped with error correction.By recognizing our problem as a d-sparse HS (Low and Chuang (2019); Babbush et al. ( 2023)), we establish that the algorithm is capable of providing exponential speedup compared to computations on classical computers.
To conclude the study, we briefly elaborate on perspectives to resolve the read-out problem, arising from the probabilistic nature of quantum mechanics.We further derive a continuous version of the classical-to-quantum map and discuss how this paves the way for employing other discretization methods.Finally, we discuss the extensions and limitations of our method and how it relates to other quantum algorithms.

A BRIEF SUMMARY OF QUANTUM COMPUTING
To set the stage, we begin with a condensed review of the tensor product formalism and the quantum circuit model that should be accessible to a geophysical audience.This framework is inherently dictated by the postulates of quantum mechanics and forms the fundamental programming model of quantum computing.A more general formalism of density matrices, as well as more in-depth introductions, can be found in the literature (Nielsen and Chuang (2010); Watrous (2018); Adedoyin et al. (2022); Benenti et al. (2019)).

Quantum particles and state space
Consider a single quantum particle characterized by its quantized properties, such as the atomic energy level.This particle can occupy N different orthogonal basis states, which span the state space of the particle.To the basis states, we assign labels Σ = {0, . . ., N − 1}, collectively referred to as the label set of size |Σ| = N .A particle with N basis states is called an N -level particle.The current state of the observable is represented by the state vector |ψ⟩ = [α0, . . ., αN−1] T ∈ C N , where α k ∈ C is the probability amplitude corresponding to the k-th basis state |k⟩, written in the bra-ket or Dirac notation that is ubiquitous in quantum mechanics.It encapsulates all information about the individual particles and their number and concisely represents exponentially large state vectors.The ket |a⟩ is a column vector, whereas the bra ⟨b| is a row vector with complex conjugate transpose ⟨b| † = |b⟩.Bra and ket form the bra-ket ⟨b|a⟩, which is the scalar product |b⟩ † • |a⟩.The main difference, compared to the common vector notation in linear algebra, is that the content of the bra/ket does not represent individual vector entries but some convenient collective label.The state vector |ψ⟩ can be expressed in terms of the basis states |k⟩ with k = 0, ..., N − 1 as the linear combination The state vector of a system composed of n particles is given by where |ψj⟩ ∈ C N j is the state of the j th particle, and Nj is the number of states that it can occupy.We reuse the symbols N and ψ, as a single-particle system is a special case of a multi-particle system.The symbol ⊗ denotes the Kronecker or tensor product.For two arbitrary 1-D elastic wave propagation on quantum computers 3 states u and v, it is defined as where Nu and Nv are the lengths of the vectors u and v, respectively.Now, we introduce the kets with integer strings.If we assume basis vectors |j⟩ and |k⟩ having label sets Σ and Γ, respectively, we have where jk is a string representing (j, k) ∈ Σ×Γ with the index I = j •|Γ|+k in the alphabetical ordering of Σ×Γ, and |jk⟩ is a vector with 1 at the I th position and zeros otherwise.The tensor product of two or more superposition states can be deduced from (3) using multi-linearity.
Quantum computing typically uses particles with two basis states, which is a natural choice because it is the minimum number for a non-trivial system and the easiest to distinguish and manipulate in a hardware implementation.Such a particle is referred to as a qubit in analogy to the classical bit.Just like a bit, a qubit is an abstract object that can be mapped to a concrete physical object in many different ways.The two standard basis states |0⟩ = [1, 0] T and |1⟩ = [0, 1] T can be represented, for instance, by the spin-up and spin-down of a spinor or two disjoint subsets of atomic energy levels.A single qubit |ψ⟩ can be in an arbitrary superposition of the two states, |ψ⟩ = β0 |0⟩ + β1 |1⟩ = [β0, β1] T ∈ C 2 , in contrast to a classical bit, which can only be in either of two states, 0 or 1.A state of an n-qubit system with qubits |ψj⟩, j = 0, ..., n − 1, is again given by (1) with N = 2 n , because Nj = 2 for every qubit.We can see that the quantum state space dimension N grows exponentially with the number of qubits n.This is the main source of the quantum computing potential.
Any quantum state must be normalized, i.e., satisfy the normalization condition Combining normalized states with the ⊗ product, the resulting composite state is normalized as well.If α k ̸ = 0 for at least two different indices k, then the state |ψ⟩ simultaneously holds multiple potential outcomes Ω = {|k⟩ : α k ̸ = 0} and is called a superposition state.
In contrast to classical systems, a measurement of the quantum system causes the state |ψ⟩ to collapse from the superposition into one of the possible basis states |k⟩ ∈ Ω randomly with probability P (k) = |α k | 2 , which explains the term probability amplitude used for the coefficients α k .The type of measurement just introduced is a full-system measurement in the standard basis, which is the simplest case; we do not deal with more general cases, which can be found in the literature at the top of this section.Measurements allow us to read the label of the collapsed state in the form of classical bits, forming the only interface between quantum and classical information.

Quantum gates and circuits
Quantum computations consist of reversible, information-preserving transformations, represented by matrices U that are unitary, i.e., that satisfy U −1 = U † , where † represents a conjugate transpose.A gate-based QC implements a set of gates, which are hardware implementations of unitary operations forming a quantum-Turing complete set.These operations can simulate an arbitrary quantum algorithm, given a sufficient time and number of qubits.The only non-unitary QC operation is the measurement.In practical quantum programming, a somewhat larger, redundant set of quantum gates is used for more convenient programming but then automatically translated into the minimum set supported by the hardware as needed before execution.These basic gates form the "quantum assembly language".They include single-qubit gates of size 2 × 2 such as the identity gate GI, the Hadamard gate GH, the NOT gate GX or rotation gates, and multi-qubit gates such as GCNOT, GSWAP ∈ C 4×4 or the Toffoli gate GCCNOT ∈ C 8×8 .
Any quantum algorithm is formed by a combination of basic gates (possibly used repeatedly) that are composed into a quantum circuit, represented by a unitary matrix W formed by a product of D layers Wj, where D is called the circuit depth.The layers are evaluated from right to left in time, meaning that the indices reflect the order of evaluation.Further, each layer Wj is formed by a tensor product of basic gates, where ti ∈ {I, H, X, CNOT, ...} is the type of the (basic) gate and ⊗ is the extension of the Kronecker product from (2) to arbitrary matrices A ∈ C M A ×N A with entries denoted as ai,j, and The number of factors in the tensor multi-product ( 6) is at most n (n if every basic gate of the layer acts on a single qubit), but it can vary between layers based on the basic gates used; the only requirement is that the total size is N so that (5) makes sense.The art of finding the most efficient circuit for the given task, i.e., the representation with the lowest n and D, is a subject of today's research in quantum computing.

Hamiltonian simulation
The dynamics of any quantum system are governed by the Schrödinger equation where i denotes the imaginary unit, and H ∈ C N ×N is a Hermitian linear operator, called the Hamiltonian, which controls the time evolution of the system.The solution of (8) can symbolically be written using the matrix exponential, with the initial state |ψ(0)⟩.The time evolution operator U(t) = exp(−iHt) is unitary.Our approach to elastic wave propagation on a QC rests on the quantum or Hamiltonian simulation (HS) approach (e.g., Georgescu et al., 2014;Nielsen and Chuang, 2010), which amounts to the direct evaluation of (9).HS was originally introduced to model multi-particle quantum dynamic systems using QCs.Due to the aforementioned exponential space growth, quantum dynamics modeling becomes prohibitively complex on classical computers.In contrast, by being quantum systems themselves, QCs naturally feature the same state space expansion, potentially making multi-particle simulations possible.
For initialization, HS requires the application of a unitary operation V that sets the desired initial state, i.e. |ψ(0)⟩ = V |00..00⟩.Subsequently, the unitary time evolution operator U(t), constructed using basic gates, is applied to |ψ(0)⟩, The main challenge of HS is that the matrix exponential U(t) = exp(−iHt) is notoriously difficult to evaluate exactly.Hence, the application of U(t) is replaced by some unitary UA(t) that evolves the initial state to the final one with the desired error ε, ||U(t) |ψ(0)⟩ − UA(t) |ψ(0)⟩ || < ε, and can be implemented efficiently using basic quantum gates, with circuit depth growing along with 1/ε.
Since the matrix exponential can be expressed in terms of a Taylor series, , we can approximate it by truncating it to a finite Taylor sum, , with some integer m < ∞.The finite sum is then further approximated by decomposing H into a sum of unitary matrices.The number of terms in this sum, determined by the given H and desired accuracy, dictates the efficiency of this particular approach (Berry et al., 2015).
Other approaches can be found in the literature, with efficiency varying based on properties of the given H.These include methods based on the Trotter-Suzuki Decomposition (e.g., Trotter, 1959;Suzuki, 1976Suzuki, , 1991;;Hatano and Suzuki, 2005;Dhand and Sanders, 2014;Yi and Crosson, 2022) and the recent qubitization approach (Low and Chuang, 2019), discussed more in detail in appendix B.
It has been shown that HS can achieve up to exponential speedup compared to analogous classical algorithms (e.g., Berry et al., 2007;Childs et al., 2018;Daley et al., 2022;Babbush et al., 2023;Low and Chuang, 2019).We will discuss this matter as well as the question of utilizing the final state in section 5.
In addition to the simulation of quantum systems, HS may also be harnessed to solve classical equations of motion, provided that an invertible transformation to the Schrödinger equation can be found (e.g., Costa et al., 2019;Suau et al., 2021;Babbush et al., 2023).Such transformations exist for classical wave equations (Süsstrunk and Huber, 2016), which share key properties with the Schrödinger equation.Both describe wave-like phenomena and are linear.In the following, we will describe such a transformation for the elastic wave equation with heterogeneous medium parameters.

ELASTIC WAVEFORM PROPAGATION USING HS
We consider a 1-D elastic wave equation in the spatial domain x ∈ [0, W ] subject to Dirichlet and Neumann boundary conditions on the right and left, respectively, with time t, the wave field u(x, t), density ρ(x) and elastic modulus µ(x).Furthermore, the wave field is subject to the initial conditions u0 = u(x, 0) and v0 = ∂tu(x, 0).We introduce a first-order finite-difference discretization in space with 2 n−1 grid points, each with position and velocity, and a grid spacing of ∆x = W/(2 n−1 − 1), thereby approximating the continuous wave field u(x, t) by the discrete wave field vector u(t) = [u(x0, t), . . ., u(x 2 n−1 −1 , t)] T .In section 5.3, we discuss the potential use of other discretization methods.We further define the positive definite mass matrix M that serves as a discrete representation of the density distribution, 1-D elastic wave propagation on quantum computers 5 and the positive definite elasticity matrix E, which is the discrete version of µ(x), To obtain the stiffness matrix K, we use a second-order stencil with the stress given by σi = µ(xi)[u(xi+1) − u(xi)]/∆x and the elastic forces by ∂xσ(xi) = [σ(xi) − σ(xi−1)]/∆x (e.g., Moczo et al., 2014).We define the first-order accurate forward finite-difference matrix D as The symmetric negative-definite stiffness matrix K is consequently obtained as This implementation naturally results in a Dirichlet and a Neumann boundary condition on the right and left, respectively, assuming two ghost nodes.In the interest of a smooth derivation, we defer the case of arbitrary boundary conditions to appendix C. With these definitions, the space-discretized wave equation reads To solve ( 16) on a QC, it needs to be transformed into the Schrödinger equation ( 8).For this, we first apply a transformation that utilizes the positive definite, diagonal square root matrix M 1/2 , which produces the mass-transformed wave field ũ = M 1/2 u and the masstransformed stiffness matrix K = M −1/2 KM −1/2 .The latter inherits the property of being symmetric negative-definite from K. We proceed by transforming the second-order system into a first-order system where ṽ = ∂t ũ = M 1/2 ∂tu ∈ R 2 n−1 is the mass-transformed velocity vector and Q the impedance matrix.While ( 17) is already closer to the Schrödinger equation than ( 16), it remains to symmetrize Q and obtain the imaginary unit i on the left-hand side.For this, we multiply (17) by an arbitrary invertible transformation T ∈ R 2 n ×2 n , which results in If T can be constructed such that H = i TQT −1 is Hermitian, we can employ H as a quantum Hamiltonian that emulates elastic wave field dynamics.The transformed wave field vector is identified as |ϕ⟩ = T[ũ, ṽ] T , which is subsequently encoded in a quantum state.The final Schrödinger equation for 1-D elastic wave propagation therefore becomes subject to the initial conditions which are typically sparse in seismic applications thanks to the spatial localization of sources.Furthermore, owing to the localized interactions of the discretized degrees of freedom, both the mass and stiffness matrices are sparse, which implies that Q is sparse, too.Eq. ( 19) can, in principle, be solved by a QC, and it is isomorphic to the discretized elastic wave equation ( 16) because T is invertible.The actual elastic wave field at all times can, therefore, be reconstructed via It remains to design a practically useful transformation T. Several transformations exist that achieve this classical-to-quantum mapping (e.g., Süsstrunk and Huber, 2016;Kane and Lubensky, 2014;Babbush et al., 2023).In the subsequent analysis, we choose a transformation that retains sparsity in both the initial conditions (20) and the Hamiltonian H.This preservation of sparsity is pivotal for scalable HS algorithms.We will revisit this topic in appendix B. To construct the desired T, we define which is an analytical Cholesky decomposition, with U being the upper triangular matrix.Since the negative-definite K is invertible, U is also invertible.This follows from det( K) = det(−U T U) = det(−U T )det(U).With this, we finally define which leads to where H is Hermitian, meaning that it can be employed as a Hamiltonian in QC simulations.This step concludes the process of mapping the elastic wave equation with heterogeneous medium parameters to the Schrödinger equation.It produces a Hamiltonian with two entries per row (2-sparse) for a finite-difference scheme of second-order accuracy in the 1-D case.We discuss higher space dimensions in section 5.

IMPLEMENTATION AND EXPERIMENTS ON QUANTUM SIMULATORS AND COMPUTERS
The transformations discussed in the previous sections are pre-processing steps executed on a classical computer.In the following, we present the implementation of the transformed elastic wave equation on both an actual QC and an ideal (error-free) quantum simulator.A quantum simulator is a software running on classical computers that emulates the behavior of actual QCs.We create a quantum circuit that defines all operations that must be physically conducted on qubits during the circuit's runtime.All of these operations, except for the measurements of the qubit states, are necessarily unitary.Hence, the quantum circuit itself can be represented by a single unitary matrix.
The quantum circuit that solves the elastic wave equation through quantum time evolution is composed of 3 sub-circuits; (1) one that prepares the initial condition |ϕ(0)⟩, (2) one that performs the quantum simulation through implementing the Hamiltonian H, and (3) one that utilises the evolved state |ϕ(t)⟩.This utilization can, for example, be a measurement of the qubits yielding binary results.In this section, we employ the open-source software development kit QISKIT by IBM Research (Aleksandrowicz et al., 2019) to create and run the quantum circuit on quantum simulators and actual QCs.The representation of the initial conditions |ϕ(0)⟩ on a QC requires normalisation, because quantum states must have a unit norm.This normalization has to be undone at the end of the simulation via Due to the linearity of ( 19), no information is lost in this procedure.We begin the simulation procedure with the initialization of an n-qubit quantum register, |00..00⟩ = |0⟩ ⊗ • • • ⊗ |0⟩ ∈ C 2 n , using the QISKIT constructor qc = QuantumCircuit(n).We must then apply the initialisation unitary transformation V to obtain the initial state, This is done by the QISKIT method qc.prepare state(psi0).Constructing and applying the quantum circuit corresponding to V can be challenging (e.g., Kak, 1999).However, in the context of seismic applications and using the above-derived transformation T, the initial state |ψ(0)⟩ is sparse and can be initialised efficiently as discussed in appendix B. Next, the quantum wave simulation circuit is realized by appending the time-evolution operator to the quantum circuit, yielding We currently perform the time evolution, i.e., the action of U, using the QISKIT class PauliEvolutionGate.It implements the scaling and squaring method (Al-Mohy and Higham, 2010).The unitary time evolution for a particular evolution time, U(t), is constructed classically and decomposed into quantum gates.These gates are implemented as part of the overall quantum circuit.Therefore, the time evolution is equivalent to a single matrix-vector product between the time evolution operator U(t) and the initial state |ψ(0)⟩.We have chosen this method because it is implementable on today's freely accessible QCs, such as IBM Brisbane, with a moderate circuit overhead for smallscale problems.For large-scale problems, which will only become tractable on QCs with higher qubit counts and better error correction, this approach becomes inapplicable due to the exponential scaling; for n qubits, U is a 2 n × 2 n dense matrix that needs to be stored and applied to |ψ(0)⟩.As already explained in section 2.3, the action of U must be carried out only approximately to obtain a scalable approach.Appendix B discusses scaling if the best available approaches are applied to our settings.Finally, the resulting quantum state |ψ(t)⟩ is utilized further within the QC or is read out using quantum state tomography (QST) (Smolin et al., 2012).Here, we opt for the latter to compare the QC simulation with the solution obtained on a classical computer.As discussed shortly, this approach is not practical for larger problems but is used here to validate the theory and implementation.
The squared entries {|ψi(t)| 2 } 2 n i=1 define a probability mass function which can be sampled by running and measuring the circuit multiple times.QST allows us to go beyond sampling the probability mass function but instead infer the actual amplitudes {ψi(t)} 2 n i=1 .
It includes two steps: Firstly, estimating the expectation value of 3 n Pauli observables, which in our case reduces to 2 n because |ψ(t)⟩ is purely real.Secondly, after all measurements are acquired, they are used in a least-squares inversion process to fit the amplitudes of |ψ(t)⟩ by determining the density matrix |ψ(t)⟩ ⟨ψ(t)| and its eigenvectors and eigenvalues (Smolin et al., 2012).Finally, it should be noted that the error in the estimated amplitudes {ψ est.i (t)} 2 n i=1 , given by ϵ = |ψi(t) − ψ est.i (t)|, decreases up to prefactors with 1/ √ S, where S is the number of samples.
QST is prohibitively expensive for large problems, due to the exponential, here 3 n (2 n ) number of observables that need to be measured for a general complex (real) state.Additionally, to obtain a single sample with respect to an observable, a quantum measurement needs to be performed, which collapses the state (Nielsen and Chuang, 2010), therefore requiring a new HS run for each sample.Reading out the final state therefore requires exponentially many HS passes, nullifying any speedup from HS.We discuss this problem and perspectives to resolve it in the context of imaging in section 5.
In summary, the quantum time evolution circuit proceeds as follows: (i) Initialization of the quantum circuit with the appropriate number of bits, each in |0⟩ state.
(ii) Preparation of the state corresponding to the initial wave field per ( 28).
(iv) Quantum state tomography (QST) via repeated HS runs and measurements.
Details on the implemented quantum circuit using QISKIT can be found in Appendix D.
The quantum state tomography output, |ψ(t)⟩, represents the quantum state after the simulation.To retrieve the classical wave field [u, v], we reverse the normalization, apply the inverse transform T −1 , and subsequently apply M −1/2 .We use a QC simulator which is run on a classical computer to investigate the correctness of our algorithm.There are two types of simulators: one where artificial noise is added to simulate the real behavior of QCs and one which is noise-free.We opt for the latter to test the theory.Therefore, any divergence from the ground truth is attributed to finite sampling noise in the QST procedure.Only in the limit of infinitely many samples, QST reproduces the actual quantum state.As the numerical ground truth, we use an ordinary differential equation (ODE) solver based on an explicit Runge-Kutta method (e.g., Dormand and Prince, 1980;Shampine, 1986) to directly integrate equation ( 34).We use Dirichlet boundary conditions, a uniform space discretization dx of 1 m, and monotonically increasing density and elastic modulus distributions.Since the quantum formalism is isomorphic to the discretized classical one, the HS should reproduce this solution, including finite-difference errors such as e.g.numerical dispersion.Figure 1 displays the simulation results of the classical ground truth (black circles) and of quantum simulators using 20 (orange) and 1000 (red) samples for the QST, respectively.We use a relative L2 error (RL2) of the wave field amplitudes u as a comparison metric calculated as with u g the ground truth solution obtained with the ODE solver.For the time steps shown in Figure 1, the RL2 error is roughly 20 % and 4 % when 20 and 1000 samples are used, respectively.Including even more samples would reduce the error further with a scaling of 1/ √ S, as discussed above.Because we use a matrix exponential scheme that evolves the state with a single matrix application (see Section 4), the error does not accumulate over time but exclusively results from the QST and is determined by the statistics of the probability mass function of the state at the particular time point.However, for more elaborate time evolution algorithms that can be implemented on future QCs, such as the one discussed in Appendix B, the error will increase over time.
Finally, we solve the problem on the real QC IBM Brisbane, requiring us to scale the problem down to 7 grid points, as discussed shortly.The results are displayed in Figure 2. We observe qualitative agreement between the numerical ground truth and the simulation on IBM Brisbane, however with RL2 errors of up to 60 %.These errors can be attributed to the high decoherence and noise rates on real QCs, as quantum states are extremely fragile and hard to isolate from the environment.We measure each quantum circuit with 1000 samples, reducing the RL2 contribution of the QST to approximately 3 %, which is far below the hardware noise level that dominates the error.Again, the error does not accumulate over time, as discussed above, but rather depends on the quantum circuit for a particular time step and its hardware implementation.
When a quantum algorithm is executed, a sequence of quantum gates is applied to each qubit as discussed in section 2. The number of gates applied to the qubits is known as the circuit depth.Each gate application has a certain success rate, the probability that the gate operation is executed correctly.IBM Brisbane's gates have an error per layered gate for a 100-qubit chain (EPLG) of 1.9 %, equaling an approximate average 5-qubit layer success rate of 99.9 % (Gambetta (2020); Collins and Easterly (2021)).For the 7 grid point problem, our implementation yields an average circuit depth of 870 quantum gates with 5 qubits, which results roughly in a probability of 0.999 870 ≈ 0.42 for a successful execution of the whole circuit.This indeed reflects approximately the error level of 60 % in our simulation.Increasing the number of grid points will increase the circuit depth and hence the error beyond the noise level, currently limiting the problem size to 7 grid points.This is not a limit imposed by the presented algorithm, but a limitation of currently available hardware.
The error characteristics of QCs vary substantially per hardware and problem.Different quantum processing units have different interaction topographies, implemented gate sets, and gate error rates of individual qubits and require continuous, often daily, re-calibration which changes the error characteristics.Because the problem specifics such as the number of grid points, evolution time, material parameters, initial  conditions, or boundary conditions alter the quantum circuit, they influence the error.However, since currently the error mostly depends on the hardware which is calibrated and updated rapidly, a detailed problem-specific error analysis is not meaningful at this stage.
While currently HS on QCs is still in the proof-of-principle state, improving the reliability of quantum chips, in combination with more advanced error mitigation methods, is a vivid research field.It is expected that substantial improvements will be achieved in the coming years and decades (Gambetta (2020); Kim et al. (2023)).Additionally, it has recently been demonstrated that HS is tractable on noisy QCs (Kikuchi et al. (2023)), which may open avenues for applications of HS in the short term.

DISCUSSION
In the following paragraphs, we provide additional details concerning potential computational speedups, the utilization of the final state, a continuous version of the above-described transformations in the context of alternative numerical methods, and future extensions of the 1-D elastic concept.
1-D elastic wave propagation on quantum computers 9
The first promise of quantum computing concerns the storage size.For classical computational methods, the number of bits required to store the entire wave field scales linearly with the number of grid points.In contrast, for quantum computational methods like the one discussed here, the number of required qubits scales only logarithmically because, reciprocally, the quantum state space grows exponentially with the number of qubits.For instance, assume a discretized wave field with two real-valued parameters per grid point, such as displacement and velocity, and Earth volume 1.083 × 10 12 km 3 .Representing a full wave field in mm 3 resolution requires merely 101 qubits.
The second promise is runtime speedup.As discussed in Section 4, the time evolution method implemented in our QC experiment has been chosen because it has a small circuit overhead, allowing for implementation on current hardware.This formulation does not provide any speedup over classical computers as it scales exponentially in the number of qubits.However, as detailed in Appendix B, the best quantum algorithms for HS have a runtime complexity that is polynomial in n and poly-logarithmic in N (Low and Chuang (2019); Babbush et al. (2023)).We analyze those algorithms only theoretically as they are not implementable on current hardware due to the large circuit overhead.
In contrast to the best HS algorithms, any classical algorithm necessarily has at least exponential complexity in n and polynomial in N (Babbush et al., 2023).In the foreseeable future, the specific N for which the quantum approach becomes faster might be above the limit given by the available hardware.However, the promise is highly enticing in the long term.

The read-out problem
As discussed in section 4, reading out the full state with QST scales polynomially in N and exponentially in n.Hence, the read-out would nullify the exponential speedup.Babbush et al. (2023) show that it is possible to read out certain properties, such as kinetic energy or potential energy of a subdomain of the model, in polynomial time, suggesting that the exponential speedup can be preserved in such cases.A similar logic will be required to preserve the exponential speedup for seismological purposes, e.g., by computing a misfit function on the fly instead of reading out the complete time evolution of the wave field.This topic will be the subject of a follow-up publication.

Continuous transformation and the use of other numerical schemes
In section 3, we performed the transformation of the elastic wave equation to the Schrödinger equation in the discrete domain in the interest of conciseness.A conceptually similar transformation can be designed in the continuous domain, thereby enabling the use of alternative numerical methods, including finite-element or spectral-element schemes.Furthermore, the continuous approach illuminates how transformations must be designed to achieve self-adjointness with respect to the natural inner product in quantum physics, loosely referred to as the quantum inner product in the following.To understand this, we consider the wave operator for the 1-D elastic wave equation We omit the dependence on space, x for brevity.Following the approach in section 3, we define a transformation utilising the square root of density ρ 1/2 , resulting in the density-transformed wave field ũ = ρ 1/2 u and a density-transformed wave operator We then transform the continuous wave equation to a first-order system where ṽ = ∂t ũ = ρ 1/2 ∂tu ∈ R are the density-transformed velocities.We now aim to construct an invertible transformation T that maps equation (34) to a continuous Schrödinger equation.Mathematically, this means that H = iT QT −1 has to be self-adjoint under the quantum inner product In appendix A, we show that there exists such a transformation that is invertible.The resulting Hamiltonian is given by and the transformed wave field vector is identified as |ϕ⟩ = T [ũ, ṽ] T .This leads to a continuous Schrödinger equation subject to the transformed initial conditions It is now possible to discretize the continuous Schrödinger equation directly by discretizing µ 1/2 ∂xρ −1/2 with any numerical method, provided that the boundary conditions can be satisfied.

Extensions and limitations
The work presented here is merely a first step toward QC wave simulations in geophysical applications.Generalizing to higher space dimensions will require a similar procedure as above to map the augmented wave operator Q or its discretized version Q to the Hamiltonian operator, ensuring self-adjointness under the quantum inner product.A bigger challenge may be the incorporation of visco-elastic attenuation because HS is inherently unitary and hence can only simulate systems that conserve energy.Another open problem is the incorporation of sources, which will be the topic of a follow-up study.

Relation to other quantum algorithms with potential for seismological research
In this article, the QC was used as a semi-analog simulator through HS.There exist other families of algorithms that have the potential to be harnessed in seismology.Among the most prominent are variational quantum algorithms (Cerezo et al. (2021)) which use a hybrid quantum-classical optimization procedure to find eigenvalues and eigenvectors (Abrams and Lloyd (1999); Kirby and Love ( 2021)) or solve linear systems (Pellow-Jarman et al. ( 2021)).The former may have applications in solving the wave equation in the frequency domain.The latter has been used to solve linear partial differential equations, including the wave equation (Trahan et al. (2023)), by successively solving the linear system of equations that arises in implicit time stepping schemes.While their computational benefits in comparison to classical solutions and HS still need to be worked out, they bring the advantage of being executable on today's and near-term quantum chips of the NISQ era (Cerezo et al. (2021)).
Again looking ahead into the potential future of error-corrected quantum chips, Harrow et al. ( 2009) have proposed a quantum algorithm, known as the HHL algorithm, that is capable of solving sparse linear systems with an exponential speedup over classical computers.Such algorithms can eventually be harnessed in a discrete-time integration scheme similar to Trahan et al. (2023) or in linearized inverse problems (Wiebe et al. (2012); Schuld et al. (2016)).Furthermore, despite concerns over vanishing gradients issues (McClean et al. (2018); Wang et al. (2021)), quantum machine learning (Biamonte et al. (2017); Melnikov et al. (2023)) is an emerging field with potential application in seismological research.
Finally, there are quantum annealers (QAs) or adiabatic quantum computers (Albash and Lidar (2018a)), which differ from gate-based QCs.QAs are capable of solving quadratic optimization problems and their application in geoscience has been investigated in a handful of studies (O'Malley (2018); Golden and O'Malley (2021); Souza et al. (2022); Dukalski et al. (2023)).It should be noted that there is controversy about whether or not QA can achieve the same speedup as gate-based quantum computation in real-world contexts (Van Dam et al. (2001); Aharonov et al. (2008); Albash and Lidar (2018b); Villanueva et al. (2023)).

CONCLUSIONS AND OUTLOOK
We presented a quantum computing concept for 1-D elastic wave propagation through heterogeneous media.Our approach involves (1) a finite-difference approximation of the wave equation, (2) a sparsity-preserving transformation to a Schrödinger equation, and (3) a direct simulation of the Schrödinger equation on a gate-based quantum computer.
The implementation of our approach on an error-free quantum simulator produces nearly exact results and serves as the foundation of wave field simulations with small numbers of grid points on the real quantum computer IBM Brisbane.While the simulation results qualitatively agree with the error-free version, they are contaminated by errors related to quantum decoherence and noise, which are typical issues for today's quantum computers.
A continuous version of the discrete transformation to the Schrödinger equation enables the modification of our method based on finite differences, including the use of other spatial discretization schemes, such as the spectral-element method.Finally, assuming that error-corrected quantum chips will become available, we discuss that the best quantum algorithms can solve the elastic wave equation with exponential speedup over classical solvers.
With the FWI as an end goal, our ongoing research focuses on incorporating sources, measuring misfit functions, extending the formalism to higher space dimensions, including visco-elastic attenuation, and computing gradients.) the application of the observable rotating the state into the computational basis, and 4.) a full state measurement of all qubits in the computational basis.

Figure 1 .
Figure 1.Comparison of elastic wave field simulations using a classical ODE solver (black) and two ideal quantum simulators with 20 samples (orange) and 1000 samples (red) used in the respective QST for a 128 grid-point problem.(a) Heterogeneous distributions of density and elastic modulus.(b-e) wave field at selected times.

Figure 2 .
Figure 2. Comparison of elastic wave field simulations using a classical ODE solver (black), an ideal quantum simulator (red), and the quantum computer IBM Brisbane (blue) with 1000 samples used in the QST for a 7 grid-point problem.(a) Heterogeneous distributions of density and elastic modulus.(b-e) wave field at selected times.
Figure A1.Circuit Representation of the wave simulation algorithm for a 7 grid point problem subject to the XZXZ Pauli observable.The circuits show four distinct algorithmic stages in differing levels of detail: 1.)The state preparation of |ψ(0)⟩, 2.) the time evolution from |ψ(0)⟩ → |ψ(t)⟩, 3.) the application of the observable rotating the state into the computational basis, and 4.) a full state measurement of all qubits in the computational basis.