Nonlinear dynamics as a ground-state solution on quantum computers

For the solution of time-dependent nonlinear differential equations, we present variational quantum algorithms (VQAs) that encode both space and time in qubit registers. The spacetime encoding enables us to obtain the entire time evolution from a single ground-state computation. We describe a general procedure to construct efficient quantum circuits for the cost function evaluation required by VQAs. To mitigate the barren plateau problem during the optimization, we propose an adaptive multigrid strategy. The approach is illustrated for the nonlinear Burgers equation. We classically optimize quantum circuits to represent the desired ground-state solutions, run them on IBM Q System One and Quantinuum System Model H1, and demonstrate that current quantum computers are capable of accurately reproducing the exact results.


I. INTRODUCTION
Methods to solve partial differential equations (PDEs), particularly those which are nonlinear, have long been of central importance in fields such as aerospace engineering [1] and energy science [2].Prominent applications include solving the Navier-Stokes equation in computational fluid dynamics and the numerical integration of continuum models in battery research [3,4].Despite their efficacy, conventional numerical methods encounter substantial constraints when tackling large-scale threedimensional models [5].
On NISQ devices, variational quantum algorithms (VQAs) have become a popular approach to solving PDEs [39].One famous example is the Variational Quantum Linear Solver, which is a generic solver for linear systems of equations and can be applied to linear PDEs [40][41][42].A standard benchmark is the Poisson equation [43][44][45][46][47].For nonlinear equations, different algorithmic primitives have been proposed [48][49][50][51][52]. * michael.lubasch@quantinuum.com, birger.horstmann@dlr.de One way of solving time-dependent PDEs using VQAs proceeds one time step at a time, analogous to established classical numerical solvers [53].In this time-stepping approach, however, a large number of time steps can lead to an accumulation of errors.Also note that time stepping leads to a computational complexity that is at least linear in the number of time steps.
In this article, our goal is to develop alternative VQAs that do not rely on time stepping.To that end, we use the formalism of Feynman [54,55] and Kitaev [56], in which time is encoded in a clock qubit register and the solution to a time-dependent problem at all points in time is contained in the ground state of a Hamiltonian.A proof-of-principle demonstration of this spacetime approach was introduced by McClean et al. for quantum chemical problems in 2013 [57], which was extended to open quantum systems by Tempel et al. in 2014 [58].A recent implementation of this method has been put forth by Barison et al. [59] for the dynamics of quantum systems, providing numerical evidence for a favorable scaling.
We present methods that extend the Feynman-Kitaev formalism to incorporate nonlinear dynamics [60].We show how to evaluate the Feynman-Kitaev Hamiltonian using quantum nonlinear processing units (QNPUs) [49,61,62].In order to mitigate barren plateaus, which are a serious problem for variational algorithms [63], and scale the variational algorithm to a larger number of qubits, we propose a multigrid optimization strategy with a customized ansatz structure [64].In this approach, a converged ansatz is passed from a coarser to a finer grid, which allows us to speed up the optimization by limiting the number of variational parameters.We apply the approach to the Burgers equation as a paradigmatic nonlinear equation.We show that the number of circuits required to evaluate the Feynman-Kitaev Hamiltonian only depends on the number of terms in the equation and the order of approximation in the time step, not on the number of time steps or spatial points, and that these circuits have a depth that scales at most linearly with the number of qubits.A schematic of the approach The nonlinear Hamiltonian is evaluated using QNPUs.The solution to the PDE at all times corresponds to the ground state of the nonlinear Hamiltonian, which can be calculated via the variational optimization of an ansatz U ( ⃗ θ) with n x + n t qubits and variational parameters ⃗ θ.(c) 3D plot of the solution corresponding to the result from the Quantinuum H1-1 computer.See Fig. 3a for more details. is depicted in Fig. 1.
This article has the following structure.In Sec.II, we first present the modified Feynman-Kitaev cost function and its implementation in terms of QNPUs.We then show how to evaluate this cost function on NISQ devices and how to optimize it using VQAs.In Sec.III, we present results from the IBMQ and Quantinuum devices.We then show the cost achieved for different variational ansatz structures and discuss the scalability of the algorithm.We provide our conclusions in Sec.IV and technical details in the Appendices.

II. METHODS
Here we describe how we use the Feynman-Kitaev formalism to solve nonlinear PDEs on NISQ devices.In Sec.II A, we derive the Feynman-Kitaev Hamiltonian for nonlinear systems.In Sec.II B, we show how to construct QNPU circuits to evaluate it.In Sec.II C, we describe the variational procedure which we use to find the ground state of the Hamiltonian.Finally, we discuss in Sec.II D what is necessary to implement this Hamiltonian on a NISQ device.

A. Feynman-Kitaev Hamiltonian for nonlinear systems
We want to apply the Feynman-Kitaev formalism to time-dependent differential equations of the form where L is a (linear or nonlinear) differential operator acting on the function f (x, t).To implement nonlinear equations, the operator L may have a functional dependence on f (x, t).Such an equation can be solved classically by defining a propagator [53] T (dt) = e dt L . ( The exponential of the operator can then be approximated using the Taylor series where the time step is assumed to be small, i.e., dt ≪ 1.Although this is formulated for a time-independent problem, the generalization to time-dependent problems using a time-ordered propagator is straightforward [57].
In contrast to quantum dynamics [59], the propagator T is not unitary for classical systems in general.
To encode the function f (x, t) in a quantum state |Ψ⟩, we define a spatial register of n x qubits and a time register of n t qubits.We then discretize the function f (x, t) on a grid of 2 nx = N x +1 points in space and 2 nt = N t +1 in time, and use the amplitude encoding [23,59] given by where binary(i) is the binary encoding of the number i as a bitstring on the space or time qubits.This quantum state is normalized as i,j |ψ i,j | 2 = 1.We now define the Feynman-Kitaev Hamiltonian, where the operator Ĉ0 , which enforces the initial condi- The constant c 0 can be tuned to speed up optimization; the ground state of the Hamiltonian, for which ⟨ Ĉ0 ⟩ = 0, does not depend on it.We use c 0 = 2 in all cases presented in this work.The normalization constant N = ∥f (x, t 0 )∥ 2 is used to rescale the results such that the norm of the first time point matches N and the initial condition is obeyed.The operator X, which describes the time evolution, is given by FIG. 2: Parametrizable quantum circuits employed throughout this work.(a) One layer of a brickwall ansatz, which can be repeated to increase the depth of the circuit.Each unitary U (red) represents a generic two-qubit gate.(b) Generic quantum MPS with bond dimension χ = 4 (Sec.II C) where each unitary U (blue) represents a generic three-qubit gate.(c) Sequential (blue) and reversed space (red) ordering of space and time qubits in the register.The arrows point from the most to the least significant qubit for time and space.(d-e) Example QNPU circuits to calculate the cost function, using the adder circuit Â [49] for the time shift operator, where Û ( ⃗ θ) is the ansatz circuit and H the Hadamard gate.Of these, (d) is used for the Ĉ1 term and (e) for the Ĉ2 term.In these two figures, the notation QNPU( ⃗ θ) is used as a short-hand for the QNPU itself and any duplicates of the ansatz, which depends on ⃗ θ, encoded on ancilla qubits.See Appendix B for details about the implementation of these circuits.This X ≈ dt(∂ t − L) describes an implicit (backward Euler) time integration scheme, known to be more stable for stiff PDEs than explicit schemes [53].The product X † X in Eq. ( 5) can be further decomposed into Ĉ1 − Ĉ2 with We now integrate the PDE of Eq. ( 1) for all time steps by finding the zero-energy ground state of Ĥ.This ground state is given by the history state We show numerically in Appendix A that the history state is indeed the ground state and that it is nondegenerate.Further details on the derivation of the cost function can also be found there.

B. Evaluation of the Hamiltonian
We evaluate the Feynman-Kitaev Hamiltonian ⟨ Ĥ⟩ described in the previous section using the quantum nonlinear processing unit (QNPU) formalism introduced by Lubasch et al. [49].A QNPU implements nonlinear terms by pointwise multiplication.Lubasch et al. present an efficient implementation of spatial derivatives and periodic boundary conditions using a quantum adder circuit [49].Nonlinear terms such as f (x, t) ∂f (x,t) ∂x in the Burgers equation (12) are represented by combining these two concepts with a duplicate of the ansatz.
We go beyond the approach of Lubasch et al. [49] by including the time register.Because our problem is not periodic in time, we add a time qubit to the register before applying the adder circuit Â for the time derivative, so the time does not wrap around from |N t ⟩ to |0⟩ (see Fig. 2e).Using QNPUs, we can estimate the expectation value ⟨ Ĥ⟩ in a number of circuits that does not depend on the number of discretization points, with a depth linear in the number of space and time qubits.In this way, we ensure the scalability of the QNPU implementation.The number of circuits to realize T and its product T † T as they appear in Eqs.(8) and (9) depends on the number of terms L = A + B + . . . in the PDE and the powers of L in the Taylor expansion of the time evolution operator T in Eq. (3).Every product of the terms A, B, etc. in the powers of L gives rise to a small number of QNPUs; see for example Eq. (B3).The number of circuits to implement the Hamiltonian for the Burgers equation depends on the order of the expansion in dt of T in Ĉ2 and T † T in Ĉ1 .Using T up to first order, which means T † T contains up to dt 2 , one can implement the Hamiltonian in 18 circuits, as shown in Appendix B. It is possible to reduce this to 11 circuits by cutting off 7 circuits encoding terms of order dt 2 in T † T .Notably, the number of circuits does not depend on the number of time or space qubits.
The Ĉ0 term enforcing the initial condition can be measured in a single circuit with a swap test [65] if one is able to prepare |ψ 0 ⟩ as a quantum state.For this swap test, efficient destructive implementations have recently been proposed [66].While preparing a function in amplitude encoding can require deep circuits in general, most engineering problems have a simple initial state that can be encoded in a shallow circuit on a quantum computer.Alternatively, one could generate this initial state using dissipative engineering [67,68].

C. Variational ansatz circuits
Throughout this study we parametrize the solution using a variational ansatz circuit where Û ( ⃗ θ) is a parametrizable quantum circuit of the forms shown in Fig. 2. We consider two different parametrizable quantum circuits as variational ansatz to find the ground state of the Hamiltonian Ĥ.First, we consider a short and dense brickwall ansatz of various depths, of which a single layer can be seen in Fig. 2a.The brickwall ansatz, also referred to as the checkerboard ansatz [46], is especially suited for NISQ devices [69].It has been used before in the context of the Variational Quantum Eigensolver (VQE) [70] as well as simulations of quantum dynamics [71].Second, we employ a tensor-inspired quantum matrix product state (MPS) ansatz, as seen in Fig. 2b.In classical simulations of dynamical systems, matrix product states represent a promising technique to improve the algorithmic efficiency thanks to their efficient representation of multidimensional problems [72,73].Quantum MPS have been employed successfully for the numerical simulation of one-dimensional systems with relatively small correlations [74].Our ansatz consists of successive three-qubit unitaries, which corresponds to a sparse representation of a classical MPS with bond dimension χ = 4.By varying the depth of the individual unitaries, we obtain an ansatz with 13, 26 or 39 CNOT gates.This is described in more detail in Appendix C.
In order to identify an optimal ordering of spacetime qubits that may take advantage of the underlying entanglement structure, we investigate two different ways to order the qubits in the register: a sequential structure with separate time and space registers (blue arrow in Fig. 2c), and one where the space register is reversed to have the qubits representing the coarsest structures in space and time next to each other (red arrow in Fig. 2c).This ordering is particularly suitable for the multigrid optimization method discussed below and provides a lowrank tensor representation [64].
In order to accelerate the minimization procedure and avoid barren plateaus, we initially use small values for coefficients such as the D or β in Eq. (12).We use the Adam optimizer [75] to find the ground state of Ĥ for this simpler problem.Then, we improve the solution in steps with the L-BFGS-B optimizer [76], as described in Appendix D. In Sec.III C, we discuss a multigrid approach to further improve the scalability, limiting the number of variational parameters.

D. Variational optimization on NISQ devices
We want to implement the nonlinear Hamiltonian (5) on NISQ devices such as IBMQ Ehningen and Quantinuum H1-1.IBMQ Ehningen is a superconducting quantum computer with 27 qubits, based on the IBM Quantum Falcon processor which allows for two-qubit gate infidelities below 10 −2 and coherence times around 100 µs [77].It features fast computation times, suitable for optimization problems with many iterations.Quantinuum H1-1 is an ion trap quantum computer with 20 qubits.It features all-to-all connectivity, a two-qubit gate infidelity around 2×10 −3 [78] and coherence times in the order of seconds [79].Therefore, complicated circuits are evaluated at high accuracy.
After performing the VQE optimization on the Qiskit statevector simulator [80], we sample the resulting state on both quantum computers.On IBMQ Ehningen, we use readout error mitigation from the M3 library [81], and transpile the circuits using the highest optimization level.This includes dynamical decoupling [82] to avoid decoherence during the idle time of the qubits.As a result, the two-qubit gate infidelity becomes the main source of errors.On Quantinuum H1-1, we run our circuit without error mitigation.The main sources of errors on this machine are qubit depolarization and the two-qubit gate infidelity [79].
We have validated the QNPU implementation for a 2 + 2-qubit problem on a noiseless simulator.For further benchmarking of our nonlinear Hamiltonian Ĥ, we decompose it into Pauli strings because this is more feasible for a low number of qubits.Nonlinear terms such as f (x, t) ∂f (x,t) ∂x in Eq. ( 12) depend on the ansatz state itself.Therefore, we perform a Pauli decomposition of a linearized Hamiltonian at each iteration of the optimization.While this is possible for small systems, the decomposition of an observable into Pauli strings scales exponentially with the number of qubits [83].

III. RESULTS AND DISCUSSION
In Sec.III A, we perform simulations on the Quantinuum and IBMQ quantum computers to show that our ansatz circuits are viable on NISQ devices.In Sec.III B, we investigate the performance of the different ansatz and entanglement structures in the minimization of the Hamiltonian Ĥ and numerically demonstrate the scaling.Finally, we discuss promising optimization strategies for our hybrid algorithm in Sec.III C and the applicability of our methods beyond variational algorithms in Sec.III D.

A. IBMQ and Quantinuum
We want to solve the Burgers equation, a nonlinear PDE that describes the convection-diffusion of particles in a viscous medium.Here, D is the diffusion coefficient and β determines the strength of the nonlinearity.If one sets β = 0, the nonlinear term vanishes and the diffusion or heat equation is obtained [84].
First, we study this equation for 3 + 3 qubits, which corresponds to a spacetime grid of 2 3 = 8 points in time and 2 3 = 8 points in space.We use β = 1 and D = 0.05 with periodic boundary conditions, with a time step of ∆t = 0.05.A detailed discussion of the discretization error in spacetime is given in Appendix E. To emphasize the nonlinearity in the form of a shock wave, we choose a Gaussian profile f (x, t 0 ) = exp(−(2πx−π) 2 ) for x ∈ [0, 1] as initial state.On a noiseless statevector simulator, we optimize for the VQE solution |ψ VQE ⟩ using a 4-layer brickwall ansatz (circuit depth 8) with reversed space entanglement structure (Fig. 2c).The cost ⟨ Ĥ⟩ of this state is low, 2.7×10 −4 .The infidelity 1−⟨ψ num |ψ VQE ⟩ = 3.3 × 10 −4 , where |ψ num ⟩ is the solution from the ODE solver from the SciPy library [85], proves that our VQE optimization has converged.
We sample this spacetime state |ψ VQE ⟩ on the IBMQ Ehningen and Quantinuum H1-1 quantum computers (see Figs. 3a and 3b).On Quantinuum H1-1, we take 2.5 × 10 4 shots, whereas on IBMQ Ehningen, we take 2 × 10 6 shots (in 20 series of 10 5 shots), which has a similar execution time.We observe that the Quantinuum computer is able to reproduce the rightwards shift of the wave, characteristic for the Burgers equation, whereas the IBM result is more affected by noise, especially in the yellow curve representing the last time point.The analogous simulation with a quantum MPS ansatz is plotted in Figs.15a and 15b, where the Quantinuum result stems from an emulator which models the noise of the Quantinuum H1-1 machine.
Next, we discard the nonlinear term from Eq. ( 12) by setting β = 0 and D = 1 with a time step of ∆t = 0.00625, which gives us the one-dimensional diffusion or heat equation.As initial condition at t 0 = 0, we consider a sinusoidal profile f (x, t 0 ) = 2+sin(2πx) for x ∈ [0, 1], with periodic boundary conditions.We have chosen the constant 2 instead of 1 for the experiments on both quantum computers, because we find that, in the presence of noise, we measure some counts for every bitstring, including the ones that correspond to values close to zero.This makes it difficult to distinguish values close to zero in amplitude encoding.This shift does not change the dynamics of the diffusion equation.For comparison, we show the unshifted simulations in Fig. 16 in Appendix E.
We calculate a VQE solution on the Qiskit statevector simulator using a 3-layer brickwall ansatz (circuit depth 6).The cost ⟨ Ĥ⟩ of this state is 4.7×10 −13 .To verify that this state with a low cost indeed matches the desired solution, we calculate the infidelity with regards to a solution from numerical integration 1−⟨ψ num |ψ VQE ⟩ = 3.2×10 −7 .These cost and infidelity values indicate that the upwards shift of the initial condition does not prevent convergence of the VQE.In Figs.3c and 3d, we sample it on the IBMQ Ehningen and Quantinuum H1-1 quantum computers.The analogous simulation using a quantum MPS ansatz is depicted in Figs.15c and 15d.

B. Circuit depth requirements
We show the achieved cost of the VQE solutions for all ansatz and entanglement structures in Fig. 4. Note that the initial condition for the diffusion equation is f (x, t 0 ) = 1 + sin(2πx) here.The circuit depth in this figure is defined by counting the consecutive layers of CNOT gates; for the brickwall ansatz with 6 qubits, this is the number of CNOT gates divided by 2.5 because one brickwall layer contains 5 blocks, but has a depth of 2. For the quantum MPS ansatz, it is equal to the number of CNOT gates.We find in Fig. 4c that our diffusion problem is very well described using a shallow reversed space ansatz with 3 or 4 brickwall layers.The low cost which we observe for these depths, and the subse- quent increase in cost for deeper circuits, indicates that these ansatzes can be optimized in less iterations than the deeper ones.The sequential ordering (see Fig. 2c) requires a few more layers and is more sensitive to random initializations.For the nonlinear problem in Fig. 4a, we see that the deeper ansatzes provide an advantage in expressibility, with only small differences between the entanglement structures.The quantum MPS ansatz requires a slightly higher number of CNOT gates, but is less sensitive to the ordering of the qubits in the register, as seen in Figs.4b and 4d.We conclude that there are particularly efficient (shallow) representations of our diffusion problem using the reversed space entanglement structure; however, the cost for these shallow circuits does show a significant spread depending on the random initial condition, indicating there are local minima to the cost function.This is consistent with literature about the trainability of shallow circuits [86,87].Our classical MPS calculations indicate that the entanglement entropy in the ansatz and the necessary ansatz depth depend significantly on the smoothness of the initial condition and the spectral content of its evolution.
We study the scaling of our algorithm.To this aim, we analyze the required circuit depth increasing the system size from 6 to 10 qubits (n x + n t = 3 + 3, 4 + 4 and 5 + 5).This gives us a measure for the necessary coherence time required to run any of these circuits on a real We have increased the number of iterations from 2, 500 to 5, 000 for every value of D for 8 qubits, and to 10, 000 for every value of D for 10 qubits.
quantum computer.The results for a 4-layer brickwall ansatz with 4 + 4 qubits and a 6-layer brickwall ansatz with 5+5 qubits are shown in Fig. 5.The costs for brickwall and quantum MPS ansatzes of different depths are shown in Fig. 6.We normalize them by dividing by the energy E 1 of the first excited state of the Hamiltonian for the respective number of qubits, such that we can use the distance to the first excited state as a measure of the quality of the solution, irrespective of the number of qubits.We observe in this figure that the low cost solutions for 3 + 3 qubits require an ansatz with more parameters than the dimension of the Hilbert space, whereas the solutions from Fig. 5 don't.From this analysis, we conclude that a crossover exists between coarse grids requiring an "overparametrized" ansatz for convergence (6 qubits), and finer grids (8 qubits or more) that can be efficiently encoded with less variational parameters than the dimension of the corresponding Hilbert space.

C. Barren plateaus and multigrid optimization
Variational cost functions commonly exhibit barren plateaus: gradients of the cost function exponentially vanish as the number of qubits increases when optimizing using random initial parameters [63].Fig. 4c indicates that it can be beneficial to optimize for a smaller depth first and then increase the depth if necessary.This motivates our approach of first solving a simpler problem with a small D and/or β and then ramping up these physical constants (see Appendix D).We calculate the variance of the gradient for different random initializations in Appendix F to evaluate this strategy.We do not observe any barren plateaus for the system sizes studied here.Now, we analyze a multigrid approach [64] for further improvements: first, the solution is found for a coarse grid using a combination of the Adam and L-BFGS optimizers, as described in Appendix D; second, optimization on a finer grid is started with the optimal parameters of the coarse grid.We use the reversed space entanglement structure.In this way, it becomes possible to add new time or space qubits while maintaining the ansatz structure for the coarser grid, as shown in Fig. 7a.This strategy is inspired by the fact that the entanglement entropy is greatest for the tensors encoding the coarsest partitions of the grid in a matrix product state representation [64].
In order to stay close to the ground state and avoid detrimental jumps in the cost when transitioning from a coarse to a fine grid, we want to implement a 'step' profile as shown in Fig. 7b.We achieve this by slightly modifying the brickwall ansatz so that every 2-qubit unitary contains two CNOT gates instead of one, and six rotations after every CNOT, duplicating the unit seen in Fig. 11d in Appendix C. There, we also show that this modification of the ansatz does not negatively affect the required circuit depth, because less layers are required if every layer contains twice the number of gates.By initializing all of these angles to zero, except for the very last rotations on the new qubits, the new qubits representing fine structures in space and time are not immediately entangled with the existing coarse ones.Importantly, the rotations at the end of the ansatz are initialized such that they represent a Hadamard gate.This results in a 'step' profile after initialization of the fine grid ansatz, as shown in Fig. 7b.We slightly shift the initial condition along the x-axis to make the generated step profile a better match.In Fig. 8a, we depict the cost after transitioning from a coarse to a fine grid.Our 'step' profile algorithm results in a smaller cost function than when initializing the new qubits randomly or in a zero state.It may be possible to engineer other initial states, such that they come even closer to the desired solution.
After adding the new qubits, we optimize for the new parameters and the blocks affecting the qubits x 2 and t 2 directly adjacent to the new ones, as shown in Fig. 7a.This guides the optimizer towards the solution of the refined problem.After this optimization round, we add the Comparison of the number of iterations versus cost for a 4 + 4-qubit problem using different optimization procedures, showing that the multigrid ansatz has a performance comparable to the Adam + L-BFGS procedure.For all three lines, the best result out of 20 random initializations was chosen.The green line describes an Adam + L-BFGS optimization for 3 + 3 qubits, followed by a transition to 4 + 4 qubits marked by the dashed black line, and subsequent L-BFGS optimization in red.This illustrates our multigrid strategy.The green and red lines for the multigrid strategy are shorter because the 3 + 3-qubit simulation considered only took 7,857 iterations instead of the total 10,000.
blocks affecting x 1 and t 1 and optimize again.Finally, we optimize all blocks.We compare the cost as a function of the number of iterations of this method with a direct optimization for 4 + 4 qubits in Fig. 8b, either using Adam only, without ramping up D on the way, or using Adam and L-BFGS.This figure shows that the multigrid strategy, with a transition from 3 + 3 to 4 + 4 qubits, needs a similar number of iterations as a direct optimization for 4 + 4 qubits with Adam and L-BFGS, whereby the multigrid method is faster because it optimizes less parameters.It is also clear that either method provides a significant advantage compared to solving with Adam only, which ends up in a local minimum almost two or-ders of magnitude higher than the other methods in this example.We conclude that the multigrid strategy described here is a promising method to solve problems on fine grids.

D. Beyond variational algorithms
The nonlinear Feynman-Kitaev Hamiltonian is applicable beyond the variational methods discussed in this article.Besides the variational quantum eigensolver (VQE), on which our methods are based, there are many other variational and non-variational techniques for find-ing the ground state.Imaginary time evolution, which we used in Appendix A 2 to study the stability of the solution, has a quantum analog that is applicable on NISQ devices [88].This quantum imaginary time evolution, which has previously been used to solve differential equations [89], is still variational.A non-variational strategy can be found in adiabatic quantum computing [90], which has similarities with the optimization strategy we used of increasing D and/or β in steps.An implementation on NISQ devices has been recently shown [91] and promising proposals exist to further compress the corresponding quantum circuits [92].As such, we expect that this can also be used to solve our ground-state problem.

IV. CONCLUSIONS
In this work, we present an extension of the Feynman-Kitaev formalism, originally developed for quantum dynamics, that is tailored to the integration of arbitrary PDEs with nonlinearities.We provide proof-of-principle calculations that demonstrate that nonlinear dissipative processes are well reproduced within this framework.Here, the full spacetime solution of the system can be retrieved in a single optimization routine of an appropriate cost function, which prevents the accumulation of errors of iterative time marching schemes.We see this as a necessary step to solve PDEs on quantum computers, as stability conditions (see Appendix A 2 for details) require that an increase in spatial resolution goes together with an increase in time resolution.Our results indicate that solutions to nonlinear PDEs can be obtained with shallow circuits using either a brickwall ansatz with a few layers or a sparse formulation of matrix product states (MPS) as a quantum circuit.We adopt a multigrid strategy where the ansatz is adaptively optimized for finer and finer solutions of the discretized PDE, thereby making the spacetime quantum solver potentially scalable in the number of qubits.We show how to adapt the ansatz so it can be expanded with additional qubits while keeping the coarser result.We expect this multigrid strategy to become more powerful as the resolution, and thereby smoothness of the function, increases.We adapt the Quantum Nonlinear Processing Unit (QNPU) [49] to evaluate the cost function.To implement time evolution, we show how to implement derivatives without periodic boundary conditions.Notably, the number of circuits to be measured is only dependent on the number of (non)linear terms in the PDE and the desired accuracy of the time evolution operator and not on the number of qubits, and the depth of the circuits is linear in the number of qubits.Depending on the nature of the problem and the available quantum resources, it can thus be favorable to use a first or second order propagator and a fine grid in time, instead of a higher order integration scheme on a coarser time grid, which is commonly done in classical simulations [53].Alternatively, one could start with a low order solution and use that as input for a higher order optimization with more circuits.
Current classical methods for nonlinear PDEs in fluid dynamics allow around 10 6 spatial points [1,72].This amounts to about 10 2 points in every dimension for a 3D grid, and is equivalent to 20 spatial qubits.Hence we believe that only around 100 qubits in total (including time qubits and duplicates of the ansatz for nonlinearities) would be necessary to go beyond current classical simulations.Our methods are also relevant in classical computing, as recent developments in tensor network representations of nonlinear problems have shown that it is possible to devise efficient classical algorithms inspired by quantum algorithms for PDEs [72,73,93].Finally, we remark that our methodology does not require a variational approach for finding the ground state of the Hamiltonian, and alternative techniques can be readily employed.

ACKNOWLEDGMENTS
This project was made possible through the support from the Competence Center Quantum Computing Baden-Württemberg in the context of the QuESt and QuESt+ projects (funded by the Ministry of Economic Affairs, Labour and Tourism Baden-Württemberg).The classical computations were performed on the JUSTUS 2 cluster, supported by the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grant no.INST 40/575-1 FUGG.We are grateful to Eric Brunner and Chris N. Self for proofreading the manuscript and valuable discussions.

Appendix A: Derivation of the cost function
We start with the general form of the differential equation as given by Eq. ( 1), where the differential operator L (linear or non-linear) is the generator of the underlying dynamics.For simplicity, let us assume that L is timeindependent.We encode the values of the function in a quantum state using the amplitude encoding of Eq. ( 4).This way we can define a quantum version of Eq. ( 1) as where the amplitudes ψ(x, t) = N f (x, t) of the state |Ψ⟩ represent the function values up to a certain normalization constant N .A general solution to this equation can be written in terms of the time evolution operator (propagator), where ⃗ ψ(t) is the state in space for a specific time t.We use this propagator to define an implicit integration scheme

Note on implicit integration
Using the implicit integration scheme of Eq. (A3), the Ĉ1 term (A7b) contains a term with T † (−dt) T (−dt).If one performs the derivation with an explicit integration scheme, there is a term T † (dt) T (dt) instead.For quantum dynamics [57,59], the propagator is unitary which means that both terms become an identity.In this case, there is no difference between an implicit and an explicit integration scheme.For our nonunitary T , we have chosen to use an implicit scheme.

Uniqueness and stability
In order to produce stable solutions in a finitedifference integration scheme, the distance travelled by the solution in one time step must be less than the distance between two points in the grid [53].For the Burgers equation (12), this stability condition translates into the requirement that D dt/dx 2 < 1/2.With our model parameters (D = 0.05, β = 1), we expect stable solutions for dt < ∼ 0.15 for 3 space qubits (dx = 1/8).To obtain stable solutions from the ground-state problem given by the Feynman-Kitaev Hamiltonian Ĥ, it is necessary that the ground state of Ĥ is unique.This is ensured if there is an energy gap between the ground state and the first excited state.In this section, we show numerical simulations with imaginary time evolution (ITE) using the mpnum library [94] to evaluate the existence of this gap under different model parameters.
The ground and first excited states of Ĥ must satisfy the conditions Due to the non-linear nature of the Hamiltonian, the states |Ψ 0 ⟩, |Ψ 1 ⟩ may not be orthogonal to each other.Furthermore, the lowest eigenvector |E ′ 0 ⟩ of Ĥ(|ϕ⟩) for an arbitrary state |ϕ⟩ may not coincide with |Ψ 0 ⟩.
The ITE method is based on the fact that the operator e −t Ĥ becomes a projector onto the ground state of Ĥ in the limit t → ∞.For an arbitrary, initial state |ϕ⟩ with a non-negligible overlap |⟨Ψ 0 |ϕ⟩| > 0 with the ground state, the ITE operator e −t Ĥ acting on |ϕ⟩ converges to the ground state Starting from a random state |ϕ⟩, we choose an Euler integration scheme for ITE by iterating We normalize the state after each time step, because the norm ⟨ϕ(t)|ϕ(t)⟩ diverges exponentially with t.
We turn the calculation of the first excited state |Ψ 1 ⟩ (A12) into another ground-state problem.This can be achieved by performing ITE with a modified Hamiltonian Ĥ′ , whose ground state energy is raised above the energy of its first excited state.One may be tempted to use the "shifted" Hamiltonian where c is a positive constant and |Ψ 0 ⟩ is the ground state calculated with ITE as described above.However, for non-linear Hamiltonians, the ground and excited states as defined in Eqs.(A11), (A12) may not be orthogonal (⟨Ψ 0 |Ψ 1 ⟩ ̸ = 0), and the shift c |Ψ 0 ⟩ ⟨Ψ 0 | in Eq. (A14) not only modifies the lowest energy subspace of Ĥ(|ϕ⟩), but it also perturbs the rest of the spectrum.To avoid this, we perform ITE simulations with the modified Hamiltonian where the second term raises the energy of the lowest eigenstate |E ′ 0 ⟩ of Ĥ(|ϕ(t)⟩), which changes with every ITE iteration.This way the ITE is directed towards a state that solves Eq. (A12).The only requisite for the constant c is that it should be larger than the second eigenvalue E ′ 1 of Ĥ(|ϕ⟩) for any |ϕ⟩.We use this method to calculate the energy gap of the nonlinear Hamiltonian Ĥ associated with the Burgers equation (12).
In Fig. 9a, we calculate the ground state for the Burgers equation ( 12) using the same parameters as in the main text (D = 0.05, β = 1), for a spacetime grid with 3 + 3 qubits (n x = n t = 3) and three different time steps (dt = 0.05, 0.1, 0.2).Our results demonstrate a stable solution for the cases with dt = 0.05 and dt = 0.1.However, for dt = 0.2, the ITE trajectories converge to an unphysical state, which is in accordance with the stability condition.In this scenario, the coefficients ψ ij (4) corresponding to later times become overrepresented.In Fig. 9b, we show E 1 as a function of the number of time qubits, for a fixed time step of dt = 0.05 and a spatial register of n x = 3 qubits.We see that for this time step, there is a gap for any number of time qubits, and that it decreases with the number of qubits for a constant dt as the number of eigenstates of the system increases.Our simulations with the modified Hamiltonian (A15) for dt = 0.2 show a constantly decreasing energy, without converging to a state that satisfies the eigenvalue condition (A12).For the states |ξ⟩ which we found on the way, Ĥ(|ξ⟩) presents a degenerate ground state.This suggests that there is indeed a closing gap for our modified Feynman-Kitaev Hamiltonian once we violate the stability condition.

Appendix B: Evaluation of the cost function
In order to measure the expectation value ⟨ Ĥ⟩ on a quantum computer, it needs to be expressed using unitary circuits.We have produced the numerical results in this work by decomposing the cost function into Pauli strings.This is not a scalable approach because the number of these strings, in general, scales exponentially with the number of qubits [83].When using a Pauli string decomposition for nonlinear problems, it is also necessary to decompose the cost function repeatedly during the optimization process to update the linearized function.In this section, we show that there is a decomposition into a constant number of circuits using quantum nonlinear processing units (QNPUs) [49], if one can efficiently create the initial state |ψ 0 ⟩ on a quantum computer.We show that the number of two-or three-qubit gates in these circuits scales linearly with the number of qubits.This is important for the scalability of the algorithm.
The QNPU approach allows us to separate space and time in the cost function, as the propagator T is independent of time.The time part of the Ĉ1 term (A7b) is then given by the sums over |i⟩⟨i| and |i + 1⟩⟨i + 1|.In general, the sum over all projectors |i⟩⟨i| is an identity.However, as these sums run until N t − 1, they represent an identity with one missing entry, given by We can then rewrite the Ĉ1 term of the cost function as This reduces the complexity of the time part of Ĉ1 to that of a single projector for both terms.The first term is determined in one circuit by sampling the ansatz and counting the results where the time is not |N t ⟩.The second term can be implemented efficiently by using an ancilla qubit to indicate whether the time was |0⟩, before performing a Hadamard test to measure T † T , then postselecting the results for times other than |0⟩.In general, the probability of this approaches 1 as n t is increased.This circuit is shown in Fig. 2d in the main text.We propose that one can combine it with a measurement of Ĉ0 (described below) as shown in Fig. 10a.
The Ĉ2 term (A7c) concerns the shift operator |i+1⟩⟨i| and its adjoint.This term can be encoded efficiently using an adder circuit, such as the one by Lubasch et al. [49], with a depth that scales linearly with the number of qubits.For small numbers of qubits one could use the one by Sato et al. [44] instead.The shift operator in Ĉ2 is not periodic, thus it is necessary to filter out the |N t ⟩⟨0| term.We achieve this by adding a qubit to our time register which we initialize to zero; instead of |N t ⟩⟨0|, there will be an |N t ⟩⟨N t + 1| term, for which the Hadamard test will measure an overlap of zero.This can be seen in Fig. 10b.
The spatial part of the Ĉ1 and Ĉ2 terms is given by the propagator T , as well as T † T .For the second-order Taylor approximation used in our main results, T contains terms up to the square of the Laplace operator, and T † T up to its fourth power.In general, one can formulate a quantum nonlinear processing unit (QNPU) as given by Lubasch et al. [49] to calculate any of the terms in the equation.In the propagator T to second order, one has the Laplacian ∂ 2 ∂x 2 , the nonlinear term f ∂ ∂x , (a) FIG. 10: (a) Circuit to evaluate Ĉ0 and the T † T term in Ĉ1 simultaneously, using an ancilla qubit to determine whether to apply a swap test for Ĉ0 if the time is |0⟩, or a QNPU for Ĉ1 otherwise.(b) Circuit to evaluate Ĉ2 , with an additional qubit for the time shift operator to make it non-periodic.This qubit represents the most significant bit of the adder.In these two figures, the notation QNPU( ⃗ θ) is used as a short-hand for the QNPU itself and any duplicates of the ansatz, which depends on ⃗ θ, encoded on ancilla qubits.(c) Building blocks of QNPU circuits, with the Â operator red and the pointwise multiplication ĉ in blue.The multi-controlled gates need to be understood as to be applied pairwise: qubit 0 on one time or space register controls qubit 0 on the other register, and so forth for every time and space qubit.(d) Circuit calculating ℜ( i ψ 2 i )/ √ N t N x with ℜ the real value.(e) Adder circuit for 4 qubits [49].
as well as the squares of both and the product of the two with each other.For T † T , one also needs to implement the third and fourth powers.The first derivative ∂ ∂x and the Laplacian ∂ 2 ∂x 2 can be produced by applying an adder circuit in space, implementing the finite difference approximations of either derivative as ( Â − Î)/∆x and ( Â − 2 Î + Â † )/(∆x) 2 , respectively.For the nonlinearity, one can create a duplicate of the ansatz and do a pointwise multiplication.The differential operator L from Eq. (1) can then be calculated using where Â represents the adder circuit and ĉ the multiplication with a duplicate of the ansatz.One can see this as a diagonal operator with the coefficients of the state on the diagonal.M is a constant with which the nonlinear term needs to be multiplied, to turn the amplitudes of the quantum state into the function values f (x, t).The value of this can be determined from the norm N of the initial state divided by the probability to measure |0⟩ time.
In this QNPU formulation, a complex-valued equation (A1) is being solved, instead of the real-valued Burgers equation, because the pointwise multiplication ĉ multiplies with the complex amplitudes of the state.To reproduce the real-valued equation, the state |Ψ⟩ has to be real-valued.This can be ensured using a real-valued ansatz [62], or by adding a penalty term to the cost function, such as Ĉ3 = c 3 [1 − ℜ( i ψ 2 i )], with ℜ the real value, which is zero if all amplitudes are real.A circuit that can be used to measure ℜ( i ψ 2 i ) up to a factor √ N t N x is shown in Fig. 10d.
In the representation of Eq. (B3), the propagator T can be expressed up to first order in 4 circuits, because one can take Â and its adjoint together.The operator T † T from the Ĉ1 can be expressed in 11 circuits up to first order.The term Î ⊗ ( Î − |N t ⟩⟨N t |) from Eq. (B2) can be measured in a single circuit, and the same holds for ψ 2 from Ĉ3 when using a complex-valued ansatz.One additional circuit is necessary to determine the probability to measure the zero time state, which is needed for the normalization constant M.This brings us to a total of 18 circuits for the Burgers equation with a first-order propagator.
The coefficients in front of all of these circuits depend on the Taylor series of the exponential of ∆t L from Eq. (3), with L from Eq. (B3).For the diffusion term, this means that the coefficients scale with powers of ∆t/(∆x) 2 , while for the nonlinear term, they scale with ∆t/∆x.It is important to make sure that these coefficients do not increase exponentially with the number of space or time qubits.A good scaling can be achieved when ∆t ∝ ∆x 2 , i.e., adding two time qubits to the system for every added space qubit.Interestingly, this coincides with the stability condition for a PDE that is first order in time and second order in space, which requires that ∆t/(∆x) 2 is bounded [53].

Initial state
The time part of Ĉ0 is the projector |0⟩⟨0|, and the space part is an identity minus a projector Î − |ψ 0 ⟩⟨ψ 0 |.If one can prepare the initial state |ψ 0 ⟩ as a quantum circuit, it is possible to do a swap test to measure the overlap between the ansatz and this initial state.This results in a linear number of controlled swap gates.One can make this more hardware-efficient using a destructive swap test [66], which only requires CNOT gates instead of controlled swap gates.However, it is necessary to sample the time register and filter out the time states that are not |0⟩.Because Ĉ0 only measures time |0⟩, and the T † T term in Ĉ1 only measures times other than |0⟩, we propose a circuit where one conditionally measures Ĉ0 or Ĉ1 , based on the state of the time qubits.This is shown in Fig. 10a.
This measurement can be simplified if the initial state is a basis vector of the computational basis, or if a basis transformation can be applied such that the initial state can be identified by a single bitstring.In this case, one can sample all qubits and determine the probability of finding time |0⟩ and a spatial state that is not |ψ 0 ⟩.In particular, the Ĉ0 term for a sine wave as initial condition could be measured using a quantum Fourier transform on the space qubits instead of a swap test.The brickwall ansatz (Fig. 11a) consists of layers of two-qubit unitaries.For the results in the main text, we used a two-qubit unitary consisting of six rotations and one CNOT gate, as given in Fig. 11d.We have also run the simulations with two CNOTs per block by duplicating this two-qubit unitary, thus creating layers that are twice as deep.The achieved cost for both diffusion and the Burgers equation can be seen in Fig. 12.We see here, that there is no notable disadvantage in most cases by using less, but deeper, layers.

MPS ansatz
We consider an MPS ansatz of bond dimension χ = 4 as given in Fig. 11b.In general, each of the threequbit unitaries in this ansatz circuit can be decomposed into six two-qubit unitaries [95].Of these six two-qubit unitaries, the last one overlaps with the first one of the next three-qubit unitary so it is left out, except at the end of the circuit.This means that a structure with five twoqubit unitaries is obtained.To limit the circuit depth, we use a 'sparse' MPS with only three two-qubit unitaries, as shown in Fig. 11c.In this circuit, we also vary the complexity of the two-qubit unitary.A generic two-qubit unitary can be composed of 15 rotations and 3 CNOT gates [66,96].We use three different depths of this circuit as shown in Fig. 11e.In Fig. 13, we present the cost as a function of circuit depth and compare with the 'full' χ = 4 quantum MPS (with five two-qubit unitaries of varying depths), as well as a χ = 2 quantum MPS, to verify that the sparse χ = 4 quantum MPS provides a good balance between circuit depth and expressibility.

Appendix D: Optimization protocol
For both the diffusion and Burgers equation, we use an optimization procedure where we start with a configuration with a small D and β = 0, and then ramp up either D or β.This is equivalent to Barison et al. [59], who instead increased the time step in steps to avoid barren plateaus.
We optimize this initial configuration with 2,500 steps of the Adam optimizer [75], and all following configurations with at most 2,500 steps L-BFGS-B optimizer [76].For the diffusion equation, we start with a diffusion coefficient of 1/8.After optimizing with Adam, we increase D in steps to 1/4, 1/2, and finally 1, optimizing with L-BFGS-B in-between.For the Burgers equation ( 12), we start with a diffusion coefficient of D = 0.05 and a nonlinearity coefficient β = 0, so the optimization with Adam can be performed on a linear equation, which is less computationally expensive.We then increase β in steps to 1/8, 1/4, 1/2, and finally 1, optimizing with L-BFGS-B in-between.This means that we end with β = 20D, which is sufficient to observe nonlinear effects.
We chose these optimizers because the Adam optimizer allows us to avoid local minima that can be encountered when starting from a random initial state, whereas the L-BFGS-B optimizer is more local, and thus well suited in cases where the initial state is already reasonably close to the solution.For either optimizer, we use automatic differentiation [97] to determine the gradients of the cost function with respect to all ansatz parameters.
For the L-BFGS-B optimization, we consider the L-BFGS-B optimization converged if the difference in cost between two successive iterations is less than 10 times machine precision.We choose a very low value here to be able to benchmark the cost function, which means that the maximum of 2,500 iterations per value of D is exhausted for most runs of the linear problem.
We run the above protocols for both types of ansatz, with different depths for the brickwall ansatz and different bond dimensions for the quantum MPS.We run each simulation 20 times with random initial conditions.We plot the cost for each random initialization as a function of the number of L-BFGS-B iterations in Fig. 14.We can see here that the simulations for shallow circuits converge around a cost of 10 −2 , while the deeper circuits either exhaust the maximum number of iterations or converge at a cost below 10 −10 , which is close to machine precision.Especially for the reversed space structure, there are many runs that converge at a low number of iterations, for all numbers of layers ≥ 3.
For shallow circuits, we find that the cost landscape is irregular, and thus a measurement protocol using different random initializations, then taking the one that leads to the lowest cost, is necessary.With deep circuits, we find that the optimization will eventually converge if the number of iterations is increased sufficiently, regardless of the initial condition.As we are focusing on shallow circuits for NISQ devices, we decided to present results from random initializations in this work.In Fig. 18, we show the results for two different stages of the optimization: after initializing randomly and after ramping up the diffusion coefficient from D = 0.5 to D = 1 (which is the last time that D is increased in the optimization procedure).Instead of taking the standard deviation of one component of the gradient, we have taken the standard deviation of the gradient norm divided by the square root of the number of parameters to be more rigorous.In this figure, there is no indication of an exponential decrease with the number of qubits.

FIG. 1 :
FIG. 1: Schematic summarizing the proposed approach.(a) The quantum state |ψ 0 ⟩ encodes the initial condition at t = 0. (b) The nonlinear Hamiltonian is evaluated using QNPUs.The solution to the PDE at all times corresponds to the ground state of the nonlinear Hamiltonian, which can be calculated via the variational optimization of an ansatz U ( ⃗ θ) with n x + n t qubits and variational parameters ⃗ θ.(c) 3D plot of the solution corresponding to the result from the Quantinuum H1-1 computer.See Fig. 3a for more details.

1 FIG. 4 :
FIG. 4: Cost ⟨ Ĥ⟩ with Ĥ from Eq. (5) as a function of the circuit depth on a noiseless simulator for 3+3 qubits.(a) Burgers equation, brickwall; (b) Burgers equation, quantum MPS; (c) diffusion, brickwall; (d) diffusion, quantum MPS.The data points in (a) and (c) represent 2 to 8 layers, those in (b) and (d) three quantum MPS ansatzes with 13 to 39 CNOT gates as detailed in Appendix C. The top and bottom axes indicate that the CNOT count and circuit depth are equivalent for the quantum MPS ansatz but different for brickwall circuits.The lines represent the median of the cost resulting from 20 different random initializations, which are plotted as scattered points.The maximum number of iterations is 4 × 2, 500 as described in Appendix D.

1 FIG. 6 :
FIG.6: Cost ⟨ Ĥ⟩ with Ĥ from Eq. (5) achieved for the diffusion equation with different numbers of qubits, divided by the energy E 1 of the first excited state, plotted as a function of number of parameters divided by Hilbert space dimension (64 for 6 qubits, 256 for 8 qubits, 1024 for 10 qubits.)The shaded area indicates where the number of parameters exceeds the dimension of the Hilbert space.The dots in the background represent the 20 individual runs for which the line shows the median.(a) Brickwall ansatz of 2 to 8 layers.The two states from Fig.5are marked with red circles.(b) Three quantum MPS circuits, as detailed in Appendix C. A reversed space entanglement structure is used in both cases.We have increased the number of iterations from 2, 500 to 5, 000 for every value of D for 8 qubits, and to 10, 000 for every value of D for 10 qubits.

1 FIG. 7 : 1 FIG. 8 :
FIG.7: Multigrid optimization strategy.(a) The reversed space ansatz of two brickwall layers is expanded from 6 to 8 qubits.The red blocks correspond to the new spacetime qubits (x 3 , t 3 ) and the optimization is initially only performed for the angles within the red and blue blocks, while the unitaries of the qubits encoding the coarser structures (x 0 , t 0 ) are fixed to the parameters of the converged circuit with 6 qubits.(b) Converged 3 + 3-qubit result expanded to 4 + 4 qubits, showing the "step" profile that is created by initializing the new qubits with a Hadamard gate.The dashed lines correspond to the solution obtained via numerical integration.

nt 5 × 10 − 3 1 × 10 − 2 2 × 10 − 2 E1 1 FIG. 9 :
FIG. 9: (a) Ground-state profiles of Ĥ for the Burgers equation (D = 0.05, β = 1) calculated with imaginary time evolution (ITE).Converged solutions can be achieved for dt = 0.05 and dt = 0.1, but the breakdown of stable solutions to the discretized PDE with dt = 0.2 results in unphysical states.(b) Energies of the first excited state E 1 (A12) of Ĥ as a function of the number of time qubits n t for a fixed time step of dt = 0.05 and a spatial register of n x = 3 qubits.

1 FIG. 12 :
FIG. 12: Comparison of the cost ⟨ Ĥ⟩ using the brickwall ansatz with r = 1 or 2 CNOTs per block, with a reversed space entanglement structure, and a brickwall ansatz with 2 to 8 layers for r = 1 and 1 to 4 layers for r = 2.The lines represent the median of the cost resulting from 20 different random initializations, which are plotted as scattered points.(a) Diffusion.(b) Burgers equation.

1 FIG. 13 : 1 FIG. 14 :
FIG.13: Comparison of the cost ⟨ Ĥ⟩ achieved using different MPS ansatz circuits with χ=2 as well as χ = 4, with three versions of the two-qubit unitary (three data points on every line) and two versions of the 3-qubit unitary (orange and green lines), with a reversed space entanglement structure.The lines represent the median of the cost resulting from 20 different random initializations, which are plotted as scattered points.(a) Diffusion.(b) Burgers equation.

1 FIG. 19 :
FIG. 19: Result of noisy optimization of the diffusion equation on 2 + 2 qubits.The optimization was done using the IBMQ Montreal noise model and the result has been plotted without noise.The dashed lines correspond to the solution obtained via numerical integration.
d) FIG.11: Schematics of a brickwall ansatz with a single layer and a generic MPS ansatz of bond dimension χ = 4. (a) One layer of the brickwall ansatz.(b) Generic quantum MPS with χ = 4. (c) Three-qubit unitary, with the part taken out for the sparse χ = 4 quantum MPS shown in red.The dashed unitary is left out if it overlaps with one of the next block.(d) Two-qubit unitary with one CNOT for the brickwall ansatz.(e) Two-qubit unitary for the MPS ansatz, for a number of CNOTs per block of r = 1...3.The dashed gate is left out if it overlaps with an R z gate of the previous block.